[codex] Make rumple a gridded metric

docs/usage/forest-structure/rumple.md

@@ -17,36 +17,46 @@ Where:
 A flat canopy has a rumple value of 1.0. More structurally complex or
 corrugated canopies have values greater than 1.0.
 
-PyForestScan computes rumple from a Canopy Height Model (CHM) by treating the
-CHM as a triangulated surface over the raster grid and summing surface area
-over valid 2x2 CHM patches.
+PyForestScan computes rumple as a gridded point-cloud metric. For each output
+cell, points inside the cell are treated as a local triangulated canopy surface.
+Rumple is then calculated as the 3D surface area of those triangles divided by
+their projected planar ground area.
 
 ## Calculating Rumple
 
 To calculate rumple:
 
 ```python
 from pyforestscan.handlers import read_lidar
-from pyforestscan.calculate import calculate_chm, calculate_rumple
+from pyforestscan.calculate import calculate_rumple
+from pyforestscan.visualize import plot_metric
 
 file_path = "../example_data/20191210_5QKB020880.laz"
 arrays = read_lidar(file_path, "EPSG:32605", hag=True)
 points = arrays[0]
 
-cell_resolution = (5.0, 5.0)
-chm, extent = calculate_chm(points, cell_resolution, interpolation="linear")
-rumple = calculate_rumple(chm, cell_resolution, min_height=2.0)
+voxel_resolution = (10.0, 10.0)
+rumple, extent = calculate_rumple(points, voxel_resolution, min_height=2.0)
 
-print(f"Rumple: {rumple:.3f}")
+plot_metric(
+    "Rumple Index",
+    rumple,
+    extent,
+    metric_name="Rumple",
+    cmap="viridis",
+)
 ```
 
 ## Notes
 
-- `calculate_rumple` returns a single scalar value, not a raster.
+- `calculate_rumple` returns a 2D raster and extent, matching the pattern used
+  by gridded metrics such as CHM.
 - `min_height` can be used to exclude low vegetation before calculating
   canopy surface complexity.
-- Interpolating the CHM before calculating rumple may fill gaps, but it can
-  also smooth the canopy surface and reduce rumple slightly.
+- Cells with fewer than three unique point locations cannot form a triangulated
+  surface and are returned as `NaN`.
+- The projected ground area is the planar area covered by valid triangles
+  inside each output cell, not the full rectangular cell footprint.
 
 ## References
 

pyforestscan/calculate.py

@@ -3,6 +3,7 @@
 from scipy.interpolate import griddata
 from scipy.stats import entropy
 from scipy import ndimage
+from scipy.spatial import Delaunay, QhullError
 from typing import List, Tuple, Optional
 
 
@@ -584,74 +585,161 @@ def calculate_chm(arr, voxel_resolution, interpolation="linear",
     return chm, extent
 
 
-def calculate_rumple(chm: np.ndarray,
-                     cell_resolution: Tuple[float, float],
-                     min_height: float | None = None) -> float:
+def calculate_rumple(arr,
+                     voxel_resolution,
+                     min_height: float | None = None) -> Tuple[np.ndarray, List]:
     """
-    Calculate the canopy rumple index from a Canopy Height Model (CHM).
+    Calculate gridded canopy rumple index from point-cloud data.
 
-    Rumple is defined here as the ratio of canopy surface area to planar
-    ground area. The CHM is treated as a triangulated surface over the raster
-    grid, and the surface area is summed over valid 2x2 CHM patches.
+    Rumple is calculated independently for each (X, Y) grid cell as the ratio
+    of canopy surface area to projected ground area. Within each cell, the
+    points are treated as a local triangulated surface using their ``X``,
+    ``Y``, and ``HeightAboveGround`` values.
 
     Args:
-        chm (np.ndarray): 2D array of canopy heights.
-        cell_resolution (tuple[float, float]): CHM cell size as (dx, dy).
-        min_height (float | None, optional): If provided, CHM cells below this
-            height threshold are masked before computing the rumple index.
-            Defaults to None.
+        arr (np.ndarray): Input structured numpy array containing point cloud
+            data with fields ``X``, ``Y``, and ``HeightAboveGround``.
+        voxel_resolution (tuple of float): The output grid resolution for the
+            X and Y dimensions, specified as ``(x_resolution, y_resolution)``.
+            A three-component voxel resolution is also accepted; the Z
+            component is ignored.
+        min_height (float | None, optional): If provided, points below this
+            height-above-ground threshold are excluded before calculating
+            rumple. Defaults to None.
 
     Returns:
-        float: Rumple index (>= 1 for valid surfaces) or NaN if no valid 2x2
-            surface patches remain after masking.
+        tuple of (np.ndarray, list): A 2D array of rumple values and the
+            spatial extent as ``[x_min, x_max, y_min, y_max]``. Cells with too
+            few valid points to form a surface are set to ``np.nan``.
 
     Raises:
-        ValueError: If the CHM is not 2D, if cell_resolution is invalid, or if
-            dx/dy are not positive.
+        ValueError: If required fields are missing, if the input is empty, or
+            if voxel_resolution is invalid.
     """
-    chm = np.asarray(chm, dtype=float)
-    if chm.ndim != 2:
-        raise ValueError(f"chm must be a 2D array (got shape {chm.shape})")
+    dtype_names = getattr(getattr(arr, "dtype", None), "names", None)
+    required_fields = {"X", "Y", "HeightAboveGround"}
+    if dtype_names is None:
+        raise ValueError("Input array must be a structured NumPy array.")
 
-    if len(cell_resolution) != 2:
-        raise ValueError("cell_resolution must be a (dx, dy) tuple")
+    missing_fields = sorted(required_fields.difference(dtype_names))
+    if missing_fields:
+        raise ValueError(
+            "Input array is missing required fields: "
+            f"{', '.join(missing_fields)}."
+        )
+
+    if len(voxel_resolution) < 2:
+        raise ValueError("voxel_resolution must contain at least X and Y resolutions")
+
+    x_resolution, y_resolution = map(float, voxel_resolution[:2])
+    if x_resolution <= 0 or y_resolution <= 0:
+        raise ValueError("voxel_resolution X and Y components must be > 0")
+
+    if len(arr) == 0:
+        raise ValueError("Input array must contain at least one point.")
 
-    dx, dy = map(float, cell_resolution)
-    if dx <= 0 or dy <= 0:
-        raise ValueError("cell_resolution components must be > 0")
+    x = np.asarray(arr["X"], dtype=float)
+    y = np.asarray(arr["Y"], dtype=float)
+    z = np.asarray(arr["HeightAboveGround"], dtype=float)
 
+    valid_points = np.isfinite(x) & np.isfinite(y) & np.isfinite(z)
     if min_height is not None:
-        chm = np.where(chm >= float(min_height), chm, np.nan)
-
-    z00 = chm[:-1, :-1]
-    z10 = chm[1:, :-1]
-    z01 = chm[:-1, 1:]
-    z11 = chm[1:, 1:]
-
-    valid = (
-        np.isfinite(z00) &
-        np.isfinite(z10) &
-        np.isfinite(z01) &
-        np.isfinite(z11)
-    )
-    if not np.any(valid):
+        valid_points &= z >= float(min_height)
+
+    x = x[valid_points]
+    y = y[valid_points]
+    z = z[valid_points]
+    if x.size == 0:
+        raise ValueError("No valid points available to calculate rumple.")
+
+    x_min, x_max = x.min(), x.max()
+    y_min, y_max = y.min(), y.max()
+
+    nx = int(np.ceil((x_max - x_min) / x_resolution))
+    ny = int(np.ceil((y_max - y_min) / y_resolution))
+    nx = max(nx, 1)
+    ny = max(ny, 1)
+
+    rumple = np.full((nx, ny), np.nan, dtype=float)
+
+    x_indices = np.floor((x - x_min) / x_resolution).astype(int)
+    y_indices = np.floor((y - y_min) / y_resolution).astype(int)
+
+    np.minimum(x_indices, nx - 1, out=x_indices)
+    np.minimum(y_indices, ny - 1, out=y_indices)
+
+    flat_indices = x_indices * ny + y_indices
+    order = np.argsort(flat_indices, kind="mergesort")
+    sorted_flat = flat_indices[order]
+    unique, first = np.unique(sorted_flat, return_index=True)
+    counts = np.diff(np.append(first, sorted_flat.size))
+
+    for flat_idx, start, count in zip(unique, first, counts):
+        point_idx = order[start:start + count]
+        xi = flat_idx // ny
+        yi = flat_idx % ny
+        rumple[xi, yi] = _calculate_tin_rumple(
+            x[point_idx],
+            y[point_idx],
+            z[point_idx],
+        )
+
+    rumple = np.flip(rumple, axis=1)
+    extent = [x_min, x_min + nx * x_resolution, y_min, y_min + ny * y_resolution]
+    return rumple, extent
+
+
+def _calculate_tin_rumple(x, y, z) -> float:
+    """
+    Calculate rumple for one point set using a local triangulated surface.
+    """
+    if len(x) < 3:
+        return np.nan
+
+    xy = np.column_stack((x, y))
+    unique_xy, inverse = np.unique(xy, axis=0, return_inverse=True)
+    if unique_xy.shape[0] < 3:
         return np.nan
 
-    # Approximate the CHM as a triangular mesh over each 2x2 raster patch.
-    tri1 = 0.5 * np.sqrt(
-        (dy * (z10 - z00)) ** 2 +
-        (dx * (z01 - z00)) ** 2 +
-        (dx * dy) ** 2
+    unique_z = np.full(unique_xy.shape[0], -np.inf, dtype=float)
+    np.maximum.at(unique_z, inverse, z)
+    valid = np.isfinite(unique_z)
+    unique_xy = unique_xy[valid]
+    unique_z = unique_z[valid]
+    if unique_xy.shape[0] < 3:
+        return np.nan
+
+    try:
+        triangulation = Delaunay(unique_xy)
+    except QhullError:
+        return np.nan
+
+    triangles = triangulation.simplices
+    p0_xy = unique_xy[triangles[:, 0]]
+    p1_xy = unique_xy[triangles[:, 1]]
+    p2_xy = unique_xy[triangles[:, 2]]
+
+    planar_area = 0.5 * np.abs(
+        (p1_xy[:, 0] - p0_xy[:, 0]) * (p2_xy[:, 1] - p0_xy[:, 1]) -
+        (p2_xy[:, 0] - p0_xy[:, 0]) * (p1_xy[:, 1] - p0_xy[:, 1])
     )
-    tri2 = 0.5 * np.sqrt(
-        (dy * (z01 - z11)) ** 2 +
-        (dx * (z11 - z10)) ** 2 +
-        (dx * dy) ** 2
+    valid_triangles = planar_area > 0
+    if not np.any(valid_triangles):
+        return np.nan
+
+    p0 = np.column_stack((p0_xy, unique_z[triangles[:, 0]]))
+    p1 = np.column_stack((p1_xy, unique_z[triangles[:, 1]]))
+    p2 = np.column_stack((p2_xy, unique_z[triangles[:, 2]]))
+    surface_area = 0.5 * np.linalg.norm(
+        np.cross(p1 - p0, p2 - p0),
+        axis=1,
     )
 
-    surface_area = np.sum((tri1 + tri2)[valid], dtype=float)
-    planar_area = float(np.count_nonzero(valid)) * dx * dy
-    return surface_area / planar_area
+    total_planar_area = np.sum(planar_area[valid_triangles], dtype=float)
+    total_surface_area = np.sum(surface_area[valid_triangles], dtype=float)
+    if total_planar_area <= 0:
+        return np.nan
+    return total_surface_area / total_planar_area
 
 
 def _calc_valid_region_mask(arr):

pyforestscan/process.py

@@ -5,7 +5,15 @@
 
 from tqdm import tqdm
 
-from pyforestscan.calculate import calculate_fhd, calculate_pad, calculate_pai, assign_voxels, calculate_chm, calculate_canopy_cover
+from pyforestscan.calculate import (
+    calculate_fhd,
+    calculate_pad,
+    calculate_pai,
+    assign_voxels,
+    calculate_chm,
+    calculate_canopy_cover,
+    calculate_rumple,
+)
 from pyforestscan.filters import remove_outliers_and_clean, downsample_poisson, downsample_voxel
 from pyforestscan.handlers import create_geotiff
 from pyforestscan.pipeline import _filter_expression, _filter_statistical_outlier, _hag_raster, _hag_delaunay
@@ -55,6 +63,7 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                        cover_min_height: float = 2.0, cover_k: float = 0.5,
                        pai_min_height: float = 1.0,
                        fhd_min_height: float = 0.0,
+                       rumple_min_height: float | None = None,
                        skip_existing: bool = False, verbose: bool = False,
                        thin_radius: float | None = None,
                        voxelgrid_cell: float | None = None,
@@ -69,7 +78,7 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
         ept_file (str): Path to the EPT file containing the point cloud data.
         tile_size (tuple): Size of each tile as (tile_width, tile_height).
         output_path (str): Directory where the output files will be saved.
-        metric (str): Metric to compute for each tile ("chm", "fhd", "pai", or "cover").
+        metric (str): Metric to compute for each tile ("chm", "fhd", "pai", "cover", or "rumple").
         voxel_size (tuple): Voxel resolution as (x_resolution, y_resolution, z_resolution).
         voxel_height (float, optional): Height of each voxel in meters. Required if metric is "fhd", "pai", or "cover".
         buffer_size (float, optional): Fractional buffer size relative to tile size (e.g., 0.1 for 10% buffer). Defaults to 0.1.
@@ -90,6 +99,8 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
         cover_k (float, optional): Beer–Lambert extinction coefficient for canopy cover. Defaults to 0.5.
         pai_min_height (float, optional): Minimum height (m) to integrate PAI. Defaults to 1.0.
         fhd_min_height (float, optional): Minimum height (m) to include in FHD entropy. Defaults to 0.0.
+        rumple_min_height (float or None, optional): Minimum height (m) to include in rumple triangulation.
+            If None, all finite HeightAboveGround values are included. Defaults to None.
         skip_existing (bool, optional): If True, skip tiles whose output file already exists. Defaults to False.
         verbose (bool, optional): If True, print warnings for empty/invalid tiles and buffer adjustments. Defaults to False.
         thin_radius (float or None, optional): If provided (> 0), apply Poisson radius-based thinning per tile before metrics.
@@ -111,7 +122,7 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
         ValueError: If an unsupported metric is requested, if buffer or voxel sizes are invalid, or required arguments are missing.
         FileNotFoundError: If the EPT or DTM file does not exist, or a required file for processing is missing.
     """
-    if metric not in ["chm", "fhd", "pai", "cover"]:
+    if metric not in ["chm", "fhd", "pai", "cover", "rumple"]:
         raise ValueError(f"Unsupported metric: {metric}")
 
     (min_z, max_z) = (None, None)
@@ -140,8 +151,8 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                     pbar.update(1)
                     continue
 
-                # Apply buffer+crop for CHM and for PAI/COVER to avoid seam artifacts.
-                if metric in ["chm", "pai", "cover"]:
+                # Apply buffer+crop for gridded metrics to avoid seam artifacts.
+                if metric in ["chm", "pai", "cover", "rumple"]:
                     current_buffer_size = buffer_size
                 else:
                     current_buffer_size = 0.0
@@ -304,8 +315,15 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                     )
 
                     create_geotiff(chm, result_file, srs, core_extent)
-                elif metric in ["fhd", "pai", "cover"]:
-                    voxels, spatial_extent = assign_voxels(tile_points, voxel_size)
+                elif metric in ["fhd", "pai", "cover", "rumple"]:
+                    if metric == "rumple":
+                        result, spatial_extent = calculate_rumple(
+                            tile_points,
+                            voxel_size,
+                            min_height=rumple_min_height,
+                        )
+                    else:
+                        voxels, spatial_extent = assign_voxels(tile_points, voxel_size)
 
                     if metric == "fhd":
                         if voxel_size[-1] <= 0: