Skip to content

analysis

occulus.analysis

Point cloud analysis — volume computation and cross-section extraction.

Available functions
  • :func:compute_volume — cut/fill volume between two surfaces
  • :func:extract_cross_section — profile extraction along a polyline
  • :func:extract_profiles — perpendicular profiles at regular intervals
Data classes
  • :class:VolumeResult — cut/fill/net volume and area statistics
  • :class:CrossSection — station-elevation profile with source points

All implementations use pure NumPy. No optional dependencies required.

CrossSection dataclass

A cross-section profile extracted from a point cloud.

Attributes:

Name Type Description
station NDArray[float64]

Station distances along the profile, shape (M,).
elevation NDArray[float64]

Elevation values at each station, shape (M,).
points NDArray[float64]

The 3D points contributing to this cross-section, shape (K, 3).
width float

Corridor half-width used to select points.
polyline NDArray[float64]

The polyline vertices defining this cross-section, shape (P, 2) or (P, 3).
Source code in src/occulus/analysis/cross_section.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
@dataclass
class CrossSection:
    """A cross-section profile extracted from a point cloud.

    Attributes
    ----------
    station : NDArray[np.float64]
        Station distances along the profile, shape (M,).
    elevation : NDArray[np.float64]
        Elevation values at each station, shape (M,).
    points : NDArray[np.float64]
        The 3D points contributing to this cross-section, shape (K, 3).
    width : float
        Corridor half-width used to select points.
    polyline : NDArray[np.float64]
        The polyline vertices defining this cross-section, shape (P, 2) or (P, 3).
    """

    station: NDArray[np.float64]
    elevation: NDArray[np.float64]
    points: NDArray[np.float64]
    width: float
    polyline: NDArray[np.float64]

VolumeResult dataclass

Result of a cut/fill volume computation.

Attributes:

Name Type Description
cut_volume float

Volume of material removed (surface above reference), in cubic units.
fill_volume float

Volume of material added (surface below reference), in cubic units.
net_volume float

Net volume change (cut - fill). Positive means net removal.
resolution float

Grid cell size used for the computation.
area float

Total planimetric area covered by the computation grid.
cut_area float

Planimetric area where the surface is above the reference.
fill_area float

Planimetric area where the surface is below the reference.
Source code in src/occulus/analysis/volume.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
@dataclass
class VolumeResult:
    """Result of a cut/fill volume computation.

    Attributes
    ----------
    cut_volume : float
        Volume of material removed (surface above reference), in cubic units.
    fill_volume : float
        Volume of material added (surface below reference), in cubic units.
    net_volume : float
        Net volume change (cut - fill). Positive means net removal.
    resolution : float
        Grid cell size used for the computation.
    area : float
        Total planimetric area covered by the computation grid.
    cut_area : float
        Planimetric area where the surface is above the reference.
    fill_area : float
        Planimetric area where the surface is below the reference.
    """

    cut_volume: float
    fill_volume: float
    net_volume: float
    resolution: float
    area: float
    cut_area: float
    fill_area: float

compute_volume(surface, reference, resolution=1.0, method='grid')

Compute cut/fill volumes between two point cloud surfaces.

Both surfaces are rasterized onto a common 2D grid. For each cell the mean elevation of each surface is computed, and the per-cell difference (surface_z - reference_z) is multiplied by the cell area to obtain volume.

Parameters:

Name Type Description Default
surface PointCloud

The measured (as-built or post-event) surface.

 required
reference PointCloud

The reference (design or pre-event) surface.

 required
resolution float

Grid cell size in coordinate units, by default 1.0.

 1.0
method str

Volume computation method. Currently only "grid" is supported.

 'grid'

Returns:

Type Description
VolumeResult

Dataclass containing cut, fill, and net volumes together with area statistics.

Raises:

Type Description
OcculusValidationError

If resolution is not positive, either cloud is empty, or an unsupported method is specified.
Source code in src/occulus/analysis/volume.py 
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def compute_volume(
    surface: PointCloud,
    reference: PointCloud,
    resolution: float = 1.0,
    method: str = "grid",
) -> VolumeResult:
    """Compute cut/fill volumes between two point cloud surfaces.

    Both surfaces are rasterized onto a common 2D grid. For each cell the
    mean elevation of each surface is computed, and the per-cell difference
    ``(surface_z - reference_z)`` is multiplied by the cell area to obtain
    volume.

    Parameters
    ----------
    surface : PointCloud
        The measured (as-built or post-event) surface.
    reference : PointCloud
        The reference (design or pre-event) surface.
    resolution : float, optional
        Grid cell size in coordinate units, by default 1.0.
    method : str, optional
        Volume computation method. Currently only ``"grid"`` is supported.

    Returns
    -------
    VolumeResult
        Dataclass containing cut, fill, and net volumes together with area
        statistics.

    Raises
    ------
    OcculusValidationError
        If ``resolution`` is not positive, either cloud is empty, or an
        unsupported ``method`` is specified.
    """
    if resolution <= 0:
        raise OcculusValidationError(f"resolution must be positive, got {resolution}")
    if surface.n_points == 0:
        raise OcculusValidationError("surface point cloud is empty")
    if reference.n_points == 0:
        raise OcculusValidationError("reference point cloud is empty")
    if method != "grid":
        raise OcculusValidationError(
            f"Unsupported method '{method}'. Only 'grid' is currently supported."
        )

    # Determine common bounding box
    all_xy = np.vstack((surface.xyz[:, :2], reference.xyz[:, :2]))
    x_min, y_min = all_xy.min(axis=0)
    x_max, y_max = all_xy.max(axis=0)

    x_edges = np.arange(x_min, x_max + resolution, resolution)
    y_edges = np.arange(y_min, y_max + resolution, resolution)
    nx = len(x_edges) - 1
    ny = len(y_edges) - 1

    if nx < 1 or ny < 1:
        raise OcculusValidationError(
            "Point clouds are too small to form a grid at the given resolution. "
            f"Extents: x=[{x_min}, {x_max}], y=[{y_min}, {y_max}], resolution={resolution}"
        )

    surface_grid = _rasterize_mean_z(surface.xyz, x_min, y_min, resolution, nx, ny)
    reference_grid = _rasterize_mean_z(reference.xyz, x_min, y_min, resolution, nx, ny)

    # Only compute where both surfaces have data
    valid = np.isfinite(surface_grid) & np.isfinite(reference_grid)
    diff = np.full_like(surface_grid, np.nan)
    diff[valid] = surface_grid[valid] - reference_grid[valid]

    cell_area = resolution * resolution

    cut_mask = valid & (diff > 0)
    fill_mask = valid & (diff < 0)

    cut_volume = float(np.nansum(diff[cut_mask]) * cell_area)
    fill_volume = float(np.nansum(np.abs(diff[fill_mask])) * cell_area)
    net_volume = cut_volume - fill_volume

    n_valid = int(valid.sum())
    n_cut = int(cut_mask.sum())
    n_fill = int(fill_mask.sum())

    logger.info(
        "compute_volume: cut=%.2f, fill=%.2f, net=%.2f (grid %dx%d, res=%.2f)",
        cut_volume,
        fill_volume,
        net_volume,
        nx,
        ny,
        resolution,
    )

    return VolumeResult(
        cut_volume=cut_volume,
        fill_volume=fill_volume,
        net_volume=net_volume,
        resolution=resolution,
        area=float(n_valid * cell_area),
        cut_area=float(n_cut * cell_area),
        fill_area=float(n_fill * cell_area),
    )

extract_cross_section(cloud, polyline, width=1.0, resolution=0.1)

Extract a cross-section profile along a polyline.

Points within width distance of the polyline are projected onto the polyline centreline. The resulting station-elevation profile is binned at the specified resolution.

Parameters:

Name Type Description Default
cloud PointCloud

Input point cloud.

 required
polyline NDArray[float64]

Polyline vertices as (P, 2) or (P, 3) array. Only XY coordinates are used for the corridor selection; Z is ignored if present.

 required
width float

Half-width of the corridor around the polyline, by default 1.0. Points within this distance of the polyline centreline are included.

 1.0
resolution float

Station spacing for the output profile, by default 0.1.

 0.1

Returns:

Type Description
CrossSection

Extracted profile with station distances and elevations.

Raises:

Type Description
OcculusValidationError

If inputs are invalid (empty cloud, degenerate polyline, non-positive width or resolution).
Source code in src/occulus/analysis/cross_section.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def extract_cross_section(
    cloud: PointCloud,
    polyline: NDArray[np.float64],
    width: float = 1.0,
    resolution: float = 0.1,
) -> CrossSection:
    """Extract a cross-section profile along a polyline.

    Points within ``width`` distance of the polyline are projected onto the
    polyline centreline. The resulting station-elevation profile is binned at
    the specified ``resolution``.

    Parameters
    ----------
    cloud : PointCloud
        Input point cloud.
    polyline : NDArray[np.float64]
        Polyline vertices as (P, 2) or (P, 3) array. Only XY coordinates are
        used for the corridor selection; Z is ignored if present.
    width : float, optional
        Half-width of the corridor around the polyline, by default 1.0.
        Points within this distance of the polyline centreline are included.
    resolution : float, optional
        Station spacing for the output profile, by default 0.1.

    Returns
    -------
    CrossSection
        Extracted profile with station distances and elevations.

    Raises
    ------
    OcculusValidationError
        If inputs are invalid (empty cloud, degenerate polyline, non-positive
        width or resolution).
    """
    _validate_inputs(cloud, polyline, width, resolution)

    polyline_2d = np.asarray(polyline[:, :2], dtype=np.float64)
    xy = cloud.xyz[:, :2]
    z = cloud.xyz[:, 2]

    # Select points within the corridor
    stations, distances = _project_onto_polyline(xy, polyline_2d)
    corridor_mask = distances <= width

    if not corridor_mask.any():
        logger.warning(
            "extract_cross_section: no points found within width=%.2f of polyline",
            width,
        )
        return CrossSection(
            station=np.array([], dtype=np.float64),
            elevation=np.array([], dtype=np.float64),
            points=np.zeros((0, 3), dtype=np.float64),
            width=width,
            polyline=polyline,
        )

    sel_stations = stations[corridor_mask]
    sel_z = z[corridor_mask]
    sel_points = cloud.xyz[corridor_mask]

    # Bin into regular station intervals
    total_length = _polyline_length(polyline_2d)
    bin_edges = np.arange(0.0, total_length + resolution, resolution)

    if len(bin_edges) < 2:
        logger.warning(
            "extract_cross_section: polyline too short for resolution=%.3f",
            resolution,
        )
        return CrossSection(
            station=np.array([], dtype=np.float64),
            elevation=np.array([], dtype=np.float64),
            points=sel_points,
            width=width,
            polyline=polyline,
        )

    bin_indices = np.clip(np.digitize(sel_stations, bin_edges) - 1, 0, len(bin_edges) - 2)

    n_bins = len(bin_edges) - 1
    station_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0
    mean_z = np.full(n_bins, np.nan, dtype=np.float64)

    for i in range(n_bins):
        mask = bin_indices == i
        if mask.any():
            mean_z[i] = float(sel_z[mask].mean())

    # Remove empty bins
    valid = np.isfinite(mean_z)
    out_station = station_centres[valid]
    out_elevation = mean_z[valid]

    logger.debug(
        "extract_cross_section: %d bins with data out of %d (width=%.2f, res=%.3f)",
        int(valid.sum()),
        n_bins,
        width,
        resolution,
    )

    return CrossSection(
        station=out_station,
        elevation=out_elevation,
        points=sel_points,
        width=width,
        polyline=polyline,
    )

extract_profiles(cloud, polyline, interval=10.0, width=5.0, resolution=0.1)

Extract perpendicular profiles at regular intervals along a polyline.

At each station along the centreline, a perpendicular cross-section of length 2 * width is extracted and sampled at resolution.

Parameters:

Name Type Description Default
cloud PointCloud

Input point cloud.

 required
polyline NDArray[float64]

Centreline polyline vertices as (P, 2) or (P, 3) array.

 required
interval float

Spacing between profile stations along the centreline, by default 10.0.

 10.0
width float

Half-width of each perpendicular profile, by default 5.0.

 5.0
resolution float

Station spacing within each profile, by default 0.1.

 0.1

Returns:

Type Description
list[CrossSection]

One CrossSection per profile station. Empty profiles (no points in corridor) are included with zero-length arrays.

Raises:

Type Description
OcculusValidationError

If inputs are invalid (empty cloud, degenerate polyline, non-positive interval, width, or resolution).
Source code in src/occulus/analysis/cross_section.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
def extract_profiles(
    cloud: PointCloud,
    polyline: NDArray[np.float64],
    interval: float = 10.0,
    width: float = 5.0,
    resolution: float = 0.1,
) -> list[CrossSection]:
    """Extract perpendicular profiles at regular intervals along a polyline.

    At each station along the centreline, a perpendicular cross-section of
    length ``2 * width`` is extracted and sampled at ``resolution``.

    Parameters
    ----------
    cloud : PointCloud
        Input point cloud.
    polyline : NDArray[np.float64]
        Centreline polyline vertices as (P, 2) or (P, 3) array.
    interval : float, optional
        Spacing between profile stations along the centreline, by default 10.0.
    width : float, optional
        Half-width of each perpendicular profile, by default 5.0.
    resolution : float, optional
        Station spacing within each profile, by default 0.1.

    Returns
    -------
    list[CrossSection]
        One ``CrossSection`` per profile station. Empty profiles (no points
        in corridor) are included with zero-length arrays.

    Raises
    ------
    OcculusValidationError
        If inputs are invalid (empty cloud, degenerate polyline, non-positive
        interval, width, or resolution).
    """
    if cloud.n_points == 0:
        raise OcculusValidationError("Cannot extract profiles from an empty cloud")
    polyline = np.asarray(polyline, dtype=np.float64)
    if polyline.ndim != 2 or polyline.shape[0] < 2 or polyline.shape[1] < 2:
        raise OcculusValidationError(
            f"polyline must be (P, 2) or (P, 3) with P >= 2, got shape {polyline.shape}"
        )
    if interval <= 0:
        raise OcculusValidationError(f"interval must be positive, got {interval}")
    if width <= 0:
        raise OcculusValidationError(f"width must be positive, got {width}")
    if resolution <= 0:
        raise OcculusValidationError(f"resolution must be positive, got {resolution}")

    polyline_2d = polyline[:, :2]
    total_length = _polyline_length(polyline_2d)

    # Generate station distances along the centreline
    station_distances = np.arange(0.0, total_length + interval * 0.5, interval)

    profiles: list[CrossSection] = []
    cloud.xyz[:, :2]

    for dist in station_distances:
        # Find the point and perpendicular direction at this station
        centre, perp = _point_and_perp_at_station(polyline_2d, dist)

        # Build a 2-point perpendicular polyline
        p1 = centre - perp * width
        p2 = centre + perp * width
        perp_line = np.array([p1, p2], dtype=np.float64)

        profile = extract_cross_section(cloud, perp_line, width=width, resolution=resolution)
        profiles.append(profile)

    logger.info(
        "extract_profiles: %d profiles at interval=%.1f along %.1fm centreline",
        len(profiles),
        interval,
        total_length,
    )

    return profiles