Merge pull request #47 from iosefa/feature/add-voxeldownsample

process: add voxel-grid downsampling with mode selection

Author: iosefa

pyforestscan/filters.py

@@ -1,8 +1,9 @@
 from typing import List
+import math
 
 from pyforestscan.handlers import _build_pdal_pipeline
 from pyforestscan.pipeline import _filter_hag, _filter_ground, _filter_statistical_outlier, _filter_smrf, \
-    _filter_radius, _select_ground
+    _filter_radius, _select_ground, _filter_voxeldownsize
 
 
 def filter_hag(arrays, lower_limit=0, upper_limit=None) -> List:
@@ -184,3 +185,49 @@ def downsample_poisson(arrays, thin_radius) -> List:
     """
     pipeline = _build_pdal_pipeline(arrays, [_filter_radius(thin_radius)])
     return pipeline.arrays
+
+
+def downsample_voxel(arrays, cell, mode) -> List:
+    """
+    Downsample point cloud arrays using PDAL's voxel-based thinning.
+
+    Uses PDAL ``filters.voxeldownsize`` to partition space into cubic voxels of
+    edge length ``cell`` and keep one representative point per occupied voxel.
+    The representative is chosen by ``mode``:
+
+      - "center": keeps the first point encountered in each voxel and overwrites
+        its X/Y/Z with the voxel center coordinates (shifts positions).
+      - "first": keeps the first point encountered in each voxel unchanged.
+
+    Args:
+        arrays (list): Point cloud arrays to be downsampled (as expected by your
+            pipeline builder).
+        cell (float): Voxel edge length in the same units as the coordinates.
+            Must be a positive, finite number (> 0).
+        mode (str): One of {"center", "first"} (case-insensitive).
+
+    Returns:
+        list: Downsampled point cloud arrays from the executed pipeline.
+
+    Raises:
+        ValueError: If ``cell`` is not positive/finite or if ``mode`` is invalid.
+        TypeError: If ``mode`` is not a string.
+
+    Notes:
+        - "center" modifies point coordinates to voxel centers; use "first" if
+          you must preserve original XY/Z.
+    """
+    if not isinstance(cell, (int, float)) or not math.isfinite(cell) or cell <= 0:
+        raise ValueError(f"'cell' must be a positive, finite number (got {cell!r}).")
+
+    if not isinstance(mode, str):
+        raise TypeError("'mode' must be a string: 'center' or 'first'.")
+
+    mode_lc = mode.lower()
+    allowed_modes = {"center", "first"}
+    if mode_lc not in allowed_modes:
+        raise ValueError(f"'mode' must be one of {sorted(allowed_modes)} (got {mode!r}).")
+
+    pipeline = _build_pdal_pipeline(arrays, [_filter_voxeldownsize(float(cell), mode_lc)])
+    return pipeline.arrays
+

pyforestscan/pipeline.py

@@ -76,14 +76,19 @@ def _filter_hag(lower_limit=0, upper_limit=None):
 
 def _filter_radius(radius):
     """
-    todo: update: this should be described as random down sample. Uses Poisson sampling.
-    Generates a filter configuration with a specified radius.
+    Build a PDAL Poisson-disk sampler configuration.
+
+    Creates a PDAL filter config that performs fixed-radius thinning in the
+    XY plane: within any circle of radius ``radius``, at most one point is
+    retained. This enforces a minimum spacing of ``radius`` between kept
+    points and helps homogenize point density prior to downstream processing.
 
     :param radius:
-        The radius to be used in the filter configuration.
-    :type radius: int
+        Sampling radius in the same units as the point coordinates
+        (typically meters). Must be > 0 for thinning to take effect.
+    :type radius: float
     :return:
-        A dictionary containing the filter type and radius.
+        A PDAL JSON fragment for ``filters.sample``.
     :rtype: dict
     """
     return {
@@ -92,6 +97,35 @@ def _filter_radius(radius):
     }
 
 
+def _filter_voxeldownsize(cell, mode):
+    """
+    Build a PDAL voxel-grid downsampling configuration.
+
+    Partitions space into cubic voxels of edge length ``cell`` and reduces
+    each occupied voxel to a single representative point. The selection
+    strategy is controlled by ``mode`` and is passed through to PDAL.
+
+    :param cell:
+        Edge length of the cubic voxel (same units as coordinates,
+        typically meters). Must be > 0 for downsampling to take effect.
+    :type cell: float
+    :param mode:
+        Strategy PDAL uses to choose the representative point per voxel.
+        Accepted values depend on your PDAL version; this string is forwarded
+        unchanged (e.g., selecting the first point encountered or a summary
+        statistic like a centroid/mean, if supported).
+    :type mode: str
+    :return:
+        A PDAL JSON fragment for ``filters.voxeldownsize``.
+    :rtype: dict
+    """
+    return {
+        "type": "filters.voxeldownsize",
+        "cell": cell,
+        "mode": mode
+    }
+
+
 def _filter_ground():
     """
     Generate a PDAL classification filter to remove ground points (classification 2).

pyforestscan/process.py

@@ -6,7 +6,7 @@
 from tqdm import tqdm
 
 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
+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.utils import get_bounds_from_ept, get_srs_from_ept
@@ -53,7 +53,9 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                        hag_dtm=False, dtm=None, bounds=None, interpolation=None, remove_outliers=False,
                        cover_min_height: float = 2.0, cover_k: float = 0.5,
                        skip_existing: bool = False, verbose: bool = False,
-                       thin_radius: float | None = None) -> None:
+                       thin_radius: float | None = None,
+                       voxelgrid_cell: float | None = None,
+                       voxelgrid_mode: str = "first") -> 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.
@@ -81,6 +83,10 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
         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.
             Units are in the same CRS as the data (e.g., meters). Defaults to None.
+        voxelgrid_cell (float or None, optional): If provided (> 0), apply PDAL voxel-grid downsampling per tile before metrics
+            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".
 
     Returns:
         None
@@ -189,7 +195,7 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                 else:
                     tile_points = arrays[0]
 
-                # Optional radius-based thinning before metrics
+                # Optional thinning before metrics
                 if thin_radius is not None and thin_radius > 0:
                     thinned = downsample_poisson([tile_points], thin_radius=thin_radius)
                     tile_points = thinned[0] if thinned else tile_points
@@ -200,6 +206,22 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                         pbar.update(1)
                         continue
 
+                if voxelgrid_cell is not None and voxelgrid_cell > 0:
+                    try:
+                        vthinned = downsample_voxel([tile_points], cell=float(voxelgrid_cell), mode=voxelgrid_mode)
+                    except Exception as e:
+                        # Fail soft per tile to keep processing going
+                        if verbose:
+                            print(f"Warning: Voxel-grid downsampling failed for tile ({i}, {j}): {e}. Proceeding without.")
+                        vthinned = None
+                    if vthinned:
+                        tile_points = vthinned[0]
+                    if tile_points.size == 0:
+                        if verbose:
+                            print(f"Warning: Tile ({i}, {j}) empty after voxel-grid thinning. Skipping.")
+                        pbar.update(1)
+                        continue
+
                 buffer_pixels_x = int(np.ceil(buffer_x / voxel_size[0]))
                 buffer_pixels_y = int(np.ceil(buffer_y / voxel_size[1]))
 

tests/test_process.py

@@ -158,6 +158,102 @@ def _fake_downsample(arrays, thin_radius):
     assert len(created_tifs) >= 1, "Expected at least one output tile for PAI."
 
 
+@patch("pyforestscan.process.downsample_voxel")
+@patch("pyforestscan.process.pdal.Pipeline")
+def test_process_with_tiles_voxelgrid_applied(mock_pipeline_cls, mock_downsample_voxel, tmp_path):
+    """
+    When voxelgrid_cell is provided, ensure downsample_voxel is called and output is produced.
+    """
+    dtype = [("X", "f8"), ("Y", "f8"), ("HeightAboveGround", "f8")]
+    pts = np.zeros(80, dtype=dtype)
+    pts["X"] = np.random.uniform(0, 10, size=80)
+    pts["Y"] = np.random.uniform(0, 10, size=80)
+    pts["HeightAboveGround"] = np.random.uniform(0, 5, size=80)
+
+    mock_pipeline = MagicMock()
+    mock_pipeline.execute.return_value = True
+    mock_pipeline.arrays = [pts]
+    mock_pipeline_cls.return_value = mock_pipeline
+
+    def _fake_voxel(arrays, cell, mode):
+        arr = arrays[0]
+        # Keep every 3rd point to mimic downsampling
+        return [arr[::3]]
+
+    mock_downsample_voxel.side_effect = _fake_voxel
+
+    out_dir = tmp_path / "test_pai_voxelgrid"
+    out_dir.mkdir()
+
+    process_with_tiles(
+        ept_file="fake_ept_path",
+        tile_size=(20, 20),
+        output_path=str(out_dir),
+        metric="pai",
+        voxel_size=(2, 2, 1),
+        voxel_height=1.0,
+        buffer_size=0.0,
+        srs="EPSG:32610",
+        hag=False,
+        hag_dtm=False,
+        dtm=None,
+        bounds=([0, 20], [0, 20], [0, 10]),
+        voxelgrid_cell=1.5,
+        voxelgrid_mode="first",
+    )
+
+    assert mock_downsample_voxel.called, "downsample_voxel should be called when voxelgrid_cell is set"
+    created_tifs = list(out_dir.glob("tile_*_pai.tif"))
+    assert len(created_tifs) >= 1, "Expected at least one output tile for PAI with voxel-grid downsampling."
+
+
+@patch("pyforestscan.process.downsample_voxel")
+@patch("pyforestscan.process.pdal.Pipeline")
+def test_process_with_tiles_voxelgrid_empty_skips(mock_pipeline_cls, mock_downsample_voxel, tmp_path):
+    """
+    If voxel-grid thinning removes all points, the tile should be skipped gracefully.
+    """
+    dtype = [("X", "f8"), ("Y", "f8"), ("HeightAboveGround", "f8")]
+    pts = np.zeros(30, dtype=dtype)
+    pts["X"] = np.random.uniform(0, 10, size=30)
+    pts["Y"] = np.random.uniform(0, 10, size=30)
+    pts["HeightAboveGround"] = np.random.uniform(0, 1, size=30)
+
+    mock_pipeline = MagicMock()
+    mock_pipeline.execute.return_value = True
+    mock_pipeline.arrays = [pts]
+    mock_pipeline_cls.return_value = mock_pipeline
+
+    def _empty_voxel(arrays, cell, mode):
+        # Simulate all points removed by voxel downsampling
+        empty = np.array([], dtype=dtype)
+        return [empty]
+
+    mock_downsample_voxel.side_effect = _empty_voxel
+
+    out_dir = tmp_path / "test_voxelgrid_empty"
+    out_dir.mkdir()
+
+    process_with_tiles(
+        ept_file="fake_ept_path",
+        tile_size=(20, 20),
+        output_path=str(out_dir),
+        metric="fhd",
+        voxel_size=(2, 2, 1),
+        buffer_size=0.0,
+        srs="EPSG:32610",
+        hag=False,
+        hag_dtm=False,
+        dtm=None,
+        bounds=([0, 20], [0, 20], [0, 10]),
+        voxelgrid_cell=1.0,
+        voxelgrid_mode="first",
+        verbose=True,
+    )
+
+    created_tifs = list(out_dir.glob("*.tif"))
+    assert len(created_tifs) == 0, "No output should be produced when voxel-grid thinning empties the tile."
+
 @patch("pyforestscan.process.pdal.Pipeline")
 def test_process_with_tiles_pai_handles_low_top_height(mock_pipeline_cls, tmp_path):
     """