Merge pull request #41 from iosefa/feat/add-thin-radius-to-process-tiles

process: add Poisson-based thinning and improve PAI tile handling

Author: iosefa

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
+from pyforestscan.filters import remove_outliers_and_clean, downsample_poisson
 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
@@ -52,7 +52,8 @@ 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,
                        cover_min_height: float = 2.0, cover_k: float = 0.5,
-                       skip_existing: bool = False, verbose: bool = False) -> None:
+                       skip_existing: bool = False, verbose: bool = False,
+                       thin_radius: float | None = None) -> 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.
@@ -78,6 +79,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.
         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.
+            Units are in the same CRS as the data (e.g., meters). Defaults to None.
 
     Returns:
         None
@@ -111,7 +114,8 @@ 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 metric == "chm":
+                # Apply buffer+crop for CHM and for PAI/COVER to avoid seam artifacts.
+                if metric in ["chm", "pai", "cover"]:
                     current_buffer_size = buffer_size
                 else:
                     current_buffer_size = 0.0
@@ -171,17 +175,31 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
 
                 tile_pipeline = pdal.Pipeline(json.dumps(tile_pipeline_json))
                 tile_pipeline.execute()
-                if remove_outliers:
-                    tile_points = remove_outliers_and_clean(tile_pipeline.arrays)[0]
-                else:
-                    tile_points = tile_pipeline.arrays[0]
 
-                if tile_points.size == 0:
+                # Extract points from pipeline output safely
+                arrays = tile_pipeline.arrays if hasattr(tile_pipeline, "arrays") else []
+                if not arrays or arrays[0].size == 0:
                     if verbose:
                         print(f"Warning: No data in tile ({i}, {j}). Skipping.")
                     pbar.update(1)
                     continue
 
+                if remove_outliers:
+                    tile_points = remove_outliers_and_clean(arrays)[0]
+                else:
+                    tile_points = arrays[0]
+
+                # Optional radius-based 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
+
+                    if tile_points.size == 0:
+                        if verbose:
+                            print(f"Warning: Tile ({i}, {j}) empty after 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]))
 
@@ -235,8 +253,13 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                         if np.all(pad == 0):
                             result = np.zeros((pad.shape[0], pad.shape[1]))
                         else:
-                            result = calculate_pai(pad, voxel_height)
-                        result = np.where(np.isfinite(result), result, 0)
+                            # Guard against empty integration range when top height < default min_height
+                            effective_max_height = pad.shape[2] * voxel_size[-1]
+                            default_min_height = 1.0
+                            if default_min_height >= effective_max_height:
+                                result = np.zeros((pad.shape[0], pad.shape[1]))
+                            else:
+                                result = calculate_pai(pad, voxel_height)
                     elif metric == "cover":
                         if not voxel_height:
                             raise ValueError(f"voxel_height is required for metric {metric}")
@@ -252,18 +275,25 @@ def process_with_tiles(ept_file, tile_size, output_path, metric, voxel_size,
                                 max_height=None,
                                 k=cover_k,
                             )
-                        result = np.where(np.isfinite(result), result, 0)
 
                     if current_buffer_size > 0:
                         if buffer_pixels_x * 2 >= result.shape[1] or buffer_pixels_y * 2 >= result.shape[0]:
-                            print(
-                                f"Warning: Buffer size exceeds {metric.upper()} dimensions for tile ({i}, {j}). "
-                                f"Adjusting buffer size."
-                            )
+                            if verbose:
+                                print(
+                                    f"Warning: Buffer size exceeds {metric.upper()} dimensions for tile ({i}, {j}). "
+                                    f"Adjusting buffer size."
+                                )
                             buffer_pixels_x = max(0, result.shape[1] // 2 - 1)
                             buffer_pixels_y = max(0, result.shape[0] // 2 - 1)
 
-                        result = result[buffer_pixels_y:-buffer_pixels_y, buffer_pixels_x:-buffer_pixels_x]
+                        # Safe crop (avoid -0 slicing)
+                        start_x = buffer_pixels_x
+                        end_x = result.shape[1] - buffer_pixels_x if buffer_pixels_x > 0 else result.shape[1]
+                        start_y = buffer_pixels_y
+                        end_y = result.shape[0] - buffer_pixels_y if buffer_pixels_y > 0 else result.shape[0]
+
+                        if end_x > start_x and end_y > start_y:
+                            result = result[start_y:end_y, start_x:end_x]
 
                     core_extent = (
                         tile_min_x + buffer_x,

tests/test_process.py

@@ -72,3 +72,127 @@ def test_process_with_tiles_no_data(mock_pipeline_cls, tmp_path):
 
     created_tifs = list(out_dir.glob("*.tif"))
     assert len(created_tifs) == 0, "No data => should produce no TIF output."
+
+
+@patch("pyforestscan.process.pdal.Pipeline")
+def test_process_with_tiles_no_data_with_thinning(mock_pipeline_cls, tmp_path):
+    """
+    If the pipeline returns no points and thinning is requested, ensure we still skip gracefully.
+    """
+    mock_pipeline = MagicMock()
+    mock_pipeline.execute.return_value = True
+    mock_pipeline.arrays = [np.array([], dtype=[('X', 'f8'), ('Y', 'f8'), ('Z', 'f8')])]
+    mock_pipeline_cls.return_value = mock_pipeline
+
+    out_dir = tmp_path / "test_fhd_thin"
+    out_dir.mkdir()
+
+    process_with_tiles(
+        ept_file="fake_ept_path",
+        tile_size=(50, 50),
+        output_path=str(out_dir),
+        metric="fhd",
+        voxel_size=(1, 1, 1),
+        buffer_size=0.0,
+        srs="EPSG:32610",
+        hag=False,
+        hag_dtm=False,
+        dtm=None,
+        bounds=([0, 50], [0, 50], [0, 50]),
+        thin_radius=1.0
+    )
+
+    created_tifs = list(out_dir.glob("*.tif"))
+    assert len(created_tifs) == 0, "No data => should produce no TIF output (even with thinning)."
+
+
+@patch("pyforestscan.process.downsample_poisson")
+@patch("pyforestscan.process.pdal.Pipeline")
+def test_process_with_tiles_thin_radius_applied(mock_pipeline_cls, mock_downsample, tmp_path):
+    """
+    When thin_radius is provided, ensure downsample_poisson is called and output is produced.
+    """
+    # Mock PDAL pipeline to return a small synthetic point set with required fields
+    dtype = [("X", "f8"), ("Y", "f8"), ("HeightAboveGround", "f8")]
+    pts = np.zeros(100, dtype=dtype)
+    pts["X"] = np.random.uniform(0, 10, size=100)
+    pts["Y"] = np.random.uniform(0, 10, size=100)
+    pts["HeightAboveGround"] = np.random.uniform(0, 5, size=100)
+
+    mock_pipeline = MagicMock()
+    mock_pipeline.execute.return_value = True
+    mock_pipeline.arrays = [pts]
+    mock_pipeline_cls.return_value = mock_pipeline
+
+    # Mock downsample to return half the points, and track calls
+    def _fake_downsample(arrays, thin_radius):
+        arr = arrays[0]
+        return [arr[::2]]
+
+    mock_downsample.side_effect = _fake_downsample
+
+    out_dir = tmp_path / "test_pai_thin"
+    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]),
+        thin_radius=1.0
+    )
+
+    # Should have called thinning at least once
+    assert mock_downsample.called, "downsample_poisson should be called when thin_radius is set"
+
+    # And produced at least one PAI output
+    created_tifs = list(out_dir.glob("tile_*_pai.tif"))
+    assert len(created_tifs) >= 1, "Expected at least one output tile for PAI."
+
+
+@patch("pyforestscan.process.pdal.Pipeline")
+def test_process_with_tiles_pai_handles_low_top_height(mock_pipeline_cls, tmp_path):
+    """
+    When the available height is below the default PAI min_height (1 m),
+    the process should not raise and should produce a tile (zeros allowed).
+    """
+    dtype = [("X", "f8"), ("Y", "f8"), ("HeightAboveGround", "f8")]
+    pts = np.zeros(50, dtype=dtype)
+    pts["X"] = np.random.uniform(0, 10, size=50)
+    pts["Y"] = np.random.uniform(0, 10, size=50)
+    # HAG strictly below 1 m to force a single Z layer when dz=1
+    pts["HeightAboveGround"] = np.random.uniform(0, 0.5, size=50)
+
+    mock_pipeline = MagicMock()
+    mock_pipeline.execute.return_value = True
+    mock_pipeline.arrays = [pts]
+    mock_pipeline_cls.return_value = mock_pipeline
+
+    out_dir = tmp_path / "test_pai_lowtop"
+    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),  # dz=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])
+    )
+
+    created_tifs = list(out_dir.glob("tile_*_pai.tif"))
+    assert len(created_tifs) >= 1, "Expected a PAI output tile even when top height < 1 m."