[codex] Prepare v0.4.1 release

README.md

@@ -14,11 +14,11 @@
 
 ## Overview
 
-PyForestScan is a Python library designed for analyzing and visualizing forest structure using airborne 3D point cloud data. The library helps derive important forest metrics such as Canopy Height, Plant Area Index (PAI), Canopy Cover, Plant Area Density (PAD), and Foliage Height Diversity (FHD).
+PyForestScan is a Python library designed for analyzing and visualizing forest structure using airborne 3D point cloud data. The library helps derive important forest metrics such as Canopy Height, Plant Area Index (PAI), Canopy Cover, Plant Area Density (PAD), Foliage Height Diversity (FHD), and Rumple Index.
 
 ## Features
 
-- **Forest Metrics**: Calculate and visualize key metrics like Canopy Height, PAI, PAD, and FHD.
+- **Forest Metrics**: Calculate and visualize key metrics like Canopy Height, PAI, PAD, FHD, Canopy Cover, and Rumple Index.
 - **Large Point Cloud Support**: Utilizes efficient data formats such as EPT for large point cloud processing.
 - **Visualization**: Create 2D and 3D visualizations of forest structure and structural metrics
 - **Extensibility**: Easily add custom filters and visualization techniques to suit your needs.

docs/usage/forest-structure/intro.md

@@ -10,6 +10,7 @@ Most calculations for the metrics follows the method outlined in Kamoske et al.
 ## Metrics
 
 * [Canopy Height Models (CHM)](chm.md)
+* [Rumple Index](rumple.md)
 * [Plant Area Density (PAD)](pad.md)
 * [Plant Area Index (PAI)](pai.md)
 * [Foliage Height Diversity (FHD)](fhd.md)

docs/usage/forest-structure/rumple.md

@@ -0,0 +1,62 @@
+# Rumple Index
+
+## Theory
+
+Rumple is a measure of canopy surface complexity. It is defined as the ratio
+of canopy surface area to projected ground area:
+
+$$
+\text{Rumple} = \frac{A_{\text{surface}}}{A_{\text{planar}}}
+$$
+
+Where:
+
+- \( A_{\text{surface}} \) is the area of the canopy surface.
+- \( A_{\text{planar}} \) is the projected ground area beneath that surface.
+
+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.
+
+## Calculating Rumple
+
+To calculate rumple:
+
+```python
+from pyforestscan.handlers import read_lidar
+from pyforestscan.calculate import calculate_chm, calculate_rumple
+
+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)
+
+print(f"Rumple: {rumple:.3f}")
+```
+
+## Notes
+
+- `calculate_rumple` returns a single scalar value, not a raster.
+- `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.
+
+## References
+
+McElhinny, Chris, Phillip Gibbons, Cris Brack, and Juergen Bauhus. 2005.
+"Forest and woodland stand structural complexity: Its definition and
+measurement." Forest Ecology and Management 218 (1-3): 1-24.
+<https://doi.org/10.1016/j.foreco.2005.08.034>.
+
+Kane, Van R., Jonathan D. Bakker, Robert J. McGaughey, James A. Lutz,
+Rolf F. Gersonde, and Jerry F. Franklin. 2010. "Examining conifer canopy
+structural complexity across forest ages and elevations with LiDAR data."
+Canadian Journal of Forest Research 40 (4): 774-787.
+<https://doi.org/10.1139/X10-064>.

mkdocs.yml

@@ -42,6 +42,7 @@ nav:
         - Forest Structure:
             - usage/forest-structure/intro.md
             - usage/forest-structure/chm.md
+            - usage/forest-structure/rumple.md
             - usage/forest-structure/pad.md
             - usage/forest-structure/pai.md
             - usage/forest-structure/fhd.md
@@ -66,4 +67,4 @@ extra_css:
     - overrides/custom.css
 
 extra_javascript:
- - https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.0/es5/tex-mml-chtml.js
+ - https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.0/es5/tex-mml-chtml.js

pyforestscan/__init__.py

@@ -4,6 +4,7 @@
     calculate_pai,
     calculate_fhd,
     calculate_chm,
+    calculate_rumple,
     calculate_point_density,
     calculate_voxel_stat,
     generate_dtm,
@@ -16,6 +17,7 @@
     "calculate_pai",
     "calculate_fhd",
     "calculate_chm",
+    "calculate_rumple",
     "calculate_point_density",
     "calculate_voxel_stat",
     "generate_dtm",

pyforestscan/calculate.py

@@ -584,6 +584,76 @@ 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:
+    """
+    Calculate the canopy rumple index from a Canopy Height Model (CHM).
+
+    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.
+
+    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.
+
+    Returns:
+        float: Rumple index (>= 1 for valid surfaces) or NaN if no valid 2x2
+            surface patches remain after masking.
+
+    Raises:
+        ValueError: If the CHM is not 2D, if cell_resolution is invalid, or if
+            dx/dy are not positive.
+    """
+    chm = np.asarray(chm, dtype=float)
+    if chm.ndim != 2:
+        raise ValueError(f"chm must be a 2D array (got shape {chm.shape})")
+
+    if len(cell_resolution) != 2:
+        raise ValueError("cell_resolution must be a (dx, dy) tuple")
+
+    dx, dy = map(float, cell_resolution)
+    if dx <= 0 or dy <= 0:
+        raise ValueError("cell_resolution components must be > 0")
+
+    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):
+        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
+    )
+    tri2 = 0.5 * np.sqrt(
+        (dy * (z01 - z11)) ** 2 +
+        (dx * (z11 - z10)) ** 2 +
+        (dx * dy) ** 2
+    )
+
+    surface_area = np.sum((tri1 + tri2)[valid], dtype=float)
+    planar_area = float(np.count_nonzero(valid)) * dx * dy
+    return surface_area / planar_area
+
+
 def _calc_valid_region_mask(arr):
     """Create valid region mask using morphological operations."""
     kernel = ndimage.generate_binary_structure(2, 1)

pyforestscan/process.py

@@ -8,7 +8,7 @@
 from pyforestscan.calculate import calculate_fhd, calculate_pad, calculate_pai, assign_voxels, calculate_chm, calculate_canopy_cover
 from pyforestscan.filters import remove_outliers_and_clean, downsample_poisson, downsample_voxel
 from pyforestscan.handlers import create_geotiff
-from pyforestscan.pipeline import _hag_raster, _hag_delaunay
+from pyforestscan.pipeline import _filter_expression, _filter_statistical_outlier, _hag_raster, _hag_delaunay
 from pyforestscan.utils import get_bounds_from_ept, get_srs_from_ept
 
 
@@ -51,13 +51,16 @@ def _crop_dtm(dtm_path, tile_min_x, tile_min_y, tile_max_x, tile_max_y):
 def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                        voxel_height=1, buffer_size=0.1, srs=None, hag=False,
                        hag_dtm=False, dtm=None, bounds=None, interpolation=None, remove_outliers=False,
+                       outlier_mean_k: int = 8, outlier_multiplier: float = 3.0,
                        cover_min_height: float = 2.0, cover_k: float = 0.5,
                        pai_min_height: float = 1.0,
                        fhd_min_height: float = 0.0,
                        skip_existing: bool = False, verbose: bool = False,
                        thin_radius: float | None = None,
                        voxelgrid_cell: float | None = None,
-                       voxelgrid_mode: str = "first") -> None:
+                       voxelgrid_mode: str = "first",
+                       tile_indices: set[tuple[int, int]] | None = None,
+                       outliers_before_hag: bool = False) -> None:
     """
     Process a large EPT point cloud by tiling, compute CHM or other metrics for each tile,
     and write the results to the specified output directory.
@@ -79,6 +82,10 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
             If None, tiling is done over the entire dataset.
         interpolation (str or None, optional): Interpolation method for CHM calculation ("linear", "cubic", "nearest", or None).
         remove_outliers (bool, optional): Whether to remove statistical outliers before calculating metrics. Defaults to False.
+        outlier_mean_k (int, optional): Number of nearest neighbors used by the statistical outlier filter.
+            Used only when remove_outliers is True. Defaults to 8.
+        outlier_multiplier (float, optional): Standard deviation multiplier used by the statistical outlier filter.
+            Used only when remove_outliers is True. Defaults to 3.0.
         cover_min_height (float, optional): Height threshold (in meters) for canopy cover (used when metric == "cover"). Defaults to 2.0.
         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.
@@ -91,6 +98,11 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
             using cell edge length of ``voxelgrid_cell``. Defaults to None.
         voxelgrid_mode (str, optional): Representative selection for voxel-grid downsampling. Common values are
             "first" (keep first point unchanged) or "center" (snap kept point to voxel center). Defaults to "first".
+        tile_indices (set[tuple[int, int]] or None, optional): If provided, process only these
+            zero-based tile indices while preserving the full tiling grid and buffer behavior.
+        outliers_before_hag (bool, optional): If True with remove_outliers=True, apply PDAL's
+            statistical outlier filter and remove classification 7 points before HAG. This can
+            avoid Delaunay failures caused by bad points. Defaults to False.
 
     Returns:
         None
@@ -124,6 +136,10 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
     with tqdm(total=total_tiles, desc="Processing tiles") as pbar:
         for i in range(num_tiles_x):
             for j in range(num_tiles_y):
+                if tile_indices is not None and (i, j) not in tile_indices:
+                    pbar.update(1)
+                    continue
+
                 # Apply buffer+crop for CHM and for PAI/COVER to avoid seam artifacts.
                 if metric in ["chm", "pai", "cover"]:
                     current_buffer_size = buffer_size
@@ -163,6 +179,15 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                     tile_bounds = ([tile_min_x, tile_max_x], [tile_min_y, tile_max_y])
                 tile_pipeline_stages = []
 
+                if remove_outliers and outliers_before_hag:
+                    tile_pipeline_stages.extend([
+                        _filter_statistical_outlier(
+                            mean_k=outlier_mean_k,
+                            multiplier=outlier_multiplier,
+                        ),
+                        _filter_expression("!(Classification == 7)"),
+                    ])
+
                 if hag:
                     tile_pipeline_stages.append(_hag_delaunay())
                 elif hag_dtm:
@@ -194,8 +219,19 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                     pbar.update(1)
                     continue
 
-                if remove_outliers:
-                    tile_points = remove_outliers_and_clean(arrays)[0]
+                if remove_outliers and not outliers_before_hag:
+                    cleaned_arrays = remove_outliers_and_clean(
+                        arrays,
+                        mean_k=outlier_mean_k,
+                        multiplier=outlier_multiplier,
+                        remove=True,
+                    )
+                    tile_points = cleaned_arrays[0]
+                    if tile_points.size == 0:
+                        if verbose:
+                            print(f"Warning: Tile ({i}, {j}) empty after outlier removal. Skipping.")
+                        pbar.update(1)
+                        continue
                 else:
                     tile_points = arrays[0]
 

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="pyforestscan",
-    version="0.3.8",
+    version="0.4.1",
     author="Joseph Emile Honour Percival",
     author_email="ipercival@gmail.com",
     description="Analyzing forest structure using aerial LiDAR data",

tests/test_calculate.py

@@ -8,7 +8,8 @@
     calculate_pai,
     calculate_fhd,
     calculate_chm,
-    calculate_canopy_cover
+    calculate_canopy_cover,
+    calculate_rumple
 )
 import math
 
@@ -628,6 +629,72 @@ def test_calculate_chm_large_heights():
     assert np.all(chm >= 1000)
 
 
+# ----------------------------
+# Tests for calculate_rumple
+# ----------------------------
+
+def test_calculate_rumple_flat_surface():
+    chm = np.full((4, 4), 12.0)
+    rumple = calculate_rumple(chm, (2.0, 3.0))
+    assert np.isclose(rumple, 1.0)
+
+
+def test_calculate_rumple_planar_slope_matches_analytical_ratio():
+    dx, dy = 2.0, 3.0
+    a, b = 0.5, 0.25
+    x = np.arange(5) * dx
+    y = np.arange(4) * dy
+    chm = a * x[:, None] + b * y[None, :]
+
+    rumple = calculate_rumple(chm, (dx, dy))
+    expected = math.sqrt(1.0 + a ** 2 + b ** 2)
+    assert np.isclose(rumple, expected)
+
+
+def test_calculate_rumple_rough_surface_exceeds_flat_surface():
+    flat = np.full((3, 3), 10.0)
+    rough = np.array([
+        [10.0, 10.0, 10.0],
+        [10.0, 15.0, 10.0],
+        [10.0, 10.0, 10.0],
+    ])
+
+    flat_rumple = calculate_rumple(flat, (1.0, 1.0))
+    rough_rumple = calculate_rumple(rough, (1.0, 1.0))
+    assert np.isclose(flat_rumple, 1.0)
+    assert rough_rumple > flat_rumple
+
+
+def test_calculate_rumple_min_height_can_remove_all_valid_patches():
+    chm = np.array([
+        [3.0, 3.0, 3.0],
+        [3.0, 1.0, 3.0],
+        [3.0, 3.0, 3.0],
+    ])
+    rumple = calculate_rumple(chm, (1.0, 1.0), min_height=2.0)
+    assert np.isnan(rumple)
+
+
+def test_calculate_rumple_returns_nan_when_no_valid_surface_exists():
+    chm = np.array([
+        [np.nan, np.nan],
+        [np.nan, np.nan],
+    ])
+    rumple = calculate_rumple(chm, (1.0, 1.0))
+    assert np.isnan(rumple)
+
+
+def test_calculate_rumple_invalid_inputs():
+    with pytest.raises(ValueError):
+        calculate_rumple(np.ones((2, 2, 2)), (1.0, 1.0))
+
+    with pytest.raises(ValueError):
+        calculate_rumple(np.ones((2, 2)), (0.0, 1.0))
+
+    with pytest.raises(ValueError):
+        calculate_rumple(np.ones((2, 2)), (1.0,))
+
+
 # ----------------------------
 # Tests for calculate_canopy_cover
 # ----------------------------