feat(preprocessing): add binary closing and hole filling option

tests/test_preprocessing.py

@@ -0,0 +1,164 @@
+import numpy as np
+
+from vesskel.config import ExtractionConfig, OutputConfig, PipelineConfig
+from vesskel.pipeline import analyze_binary_image
+
+
+class TestClosing:
+    def test_closing_bridges_1px_gap_in_thick_line(self):
+        gap = np.zeros((16, 16), dtype=np.uint8)
+        gap[7:10, 4:7] = 1
+        gap[7:10, 8:12] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(closing_iterations=1),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(gap, "gap", config)
+        # closing bridges the 1px gap at column 7 in the 3px-thick line
+        assert result.skeleton[7:10, 7].any()
+
+    def test_closing_bridges_gap_in_thick_region(self):
+        block = np.zeros((16, 16), dtype=np.uint8)
+        block[6:10, 4:7] = 1
+        block[6:10, 8:12] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(closing_iterations=1),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(block, "block", config)
+        # the 1px gap at column 7 should be bridged before thinning
+        assert result.skeleton[6:10, 7].any()
+
+    def test_closing_off_does_not_modify(self):
+        gap = np.zeros((16, 16), dtype=np.uint8)
+        gap[7:10, 4:7] = 1
+        gap[7:10, 8:12] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(closing_iterations=0),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(gap, "gap", config)
+        # without closing the gap remains
+        assert not result.skeleton[7:10, 7].any()
+
+
+class TestFillHoles:
+    def test_fill_holes_fills_enclosed_void(self):
+        ring = np.ones((16, 16), dtype=np.uint8)
+        ring[6:10, 6:10] = 0
+        ring[6:10, 6] = 1
+        ring[6:10, 9] = 1
+        ring[6, 6:10] = 1
+        ring[9, 6:10] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(fill_holes=True),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(ring, "ring", config)
+        # preprocessed_binary should have the hole filled
+        assert result.preprocessed_binary is not None
+        assert result.preprocessed_binary[7:9, 7:9].all()
+
+
+class TestMaxHoleSize:
+    def test_max_hole_size_skips_large_holes(self):
+        ring = np.ones((32, 32), dtype=np.uint8)
+        ring[8:24, 8:24] = 0
+        ring[8:24, 8] = 1
+        ring[8:24, 23] = 1
+        ring[8, 8:24] = 1
+        ring[23, 8:24] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(fill_holes=True, max_hole_size=50),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(ring, "ring", config)
+        # 15x15 = 225 pixels, above threshold, so preprocessed binary keeps it as a hole
+        assert result.preprocessed_binary is not None
+        assert not result.preprocessed_binary[9:23, 9:23].any()
+
+    def test_max_hole_size_fills_small_holes(self):
+        ring = np.ones((16, 16), dtype=np.uint8)
+        ring[7:9, 7:9] = 0
+        ring[7:9, 7] = 1
+        ring[7:9, 8] = 1
+        ring[7, 7:9] = 1
+        ring[8, 7:9] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(fill_holes=True, max_hole_size=50),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(ring, "ring", config)
+        # 2x2 = 4 pixels, under threshold, so filled in preprocessed
+        assert result.preprocessed_binary is not None
+        assert result.preprocessed_binary[7:9, 7:9].all()
+
+    def test_max_hole_size_zero_fills_all(self):
+        ring = np.ones((16, 16), dtype=np.uint8)
+        ring[7:11, 7:11] = 0
+        ring[7:11, 7] = 1
+        ring[7:11, 10] = 1
+        ring[7, 7:11] = 1
+        ring[10, 7:11] = 1
+
+        config = PipelineConfig(
+            extraction=ExtractionConfig(fill_holes=True, max_hole_size=0),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(ring, "ring", config)
+        # zero = unlimited, so 4x4 hole gets filled
+        assert result.preprocessed_binary is not None
+        assert result.preprocessed_binary[8:10, 8:10].all()
+
+
+class TestPreprocessedBinary:
+    def test_set_when_closing_active(self):
+        img = np.zeros((16, 16), dtype=np.uint8)
+        img[8, 4:12] = 1
+        config = PipelineConfig(
+            extraction=ExtractionConfig(closing_iterations=1),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(img, "test", config)
+        assert result.preprocessed_binary is not None
+        assert result.preprocessed_binary.shape == img.shape
+
+    def test_set_when_fill_holes_active(self):
+        ring = np.ones((16, 16), dtype=np.uint8)
+        ring[7:9, 7:9] = 0
+        ring[7:9, 7] = 1
+        ring[7:9, 8] = 1
+        ring[7, 7:9] = 1
+        ring[8, 7:9] = 1
+        config = PipelineConfig(
+            extraction=ExtractionConfig(fill_holes=True),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(ring, "test", config)
+        assert result.preprocessed_binary is not None
+
+    def test_none_when_preprocessing_disabled(self):
+        img = np.zeros((16, 16), dtype=np.uint8)
+        img[8, 4:12] = 1
+        config = PipelineConfig(
+            extraction=ExtractionConfig(),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(img, "test", config)
+        assert result.preprocessed_binary is None
+
+    def test_none_when_show_preprocessed_checked_but_no_preprocessing(self):
+        img = np.zeros((16, 16), dtype=np.uint8)
+        img[8, 4:12] = 1
+        config = PipelineConfig(
+            extraction=ExtractionConfig(show_preprocessed=True),
+            output=OutputConfig(),
+        )
+        result = analyze_binary_image(img, "test", config)
+        assert result.preprocessed_binary is None

vesskel/_napari.py

@@ -48,6 +48,10 @@ def _extraction_params(
             include_vessel_radius: bool = False,
             junction_cleanup: bool = False,
             cleanup_threshold_factor: float = 2.5,
+            fill_holes: bool = False,
+            closing_iterations: int = 0,
+            max_hole_size: int = 0,
+            show_preprocessed: bool = False,
         ) -> None:
             return None
 
@@ -67,6 +71,24 @@ def _extraction_params(
                 "max": 10.0,
                 "step": 0.1,
             },
+            fill_holes={"annotation": bool, "value": False},
+            closing_iterations={
+                "annotation": int,
+                "value": 0,
+                "widget_type": "SpinBox",
+                "min": 0,
+                "max": 10,
+                "step": 1,
+            },
+            max_hole_size={
+                "annotation": int,
+                "value": 0,
+                "widget_type": "SpinBox",
+                "min": 0,
+                "max": 100000,
+                "step": 100,
+            },
+            show_preprocessed={"annotation": bool, "value": False},
         )
 
         # ---------- output parameters (magicgui) ----------
@@ -129,6 +151,23 @@ def _output_params(
         advanced_group.append(self.junction_cleanup_widget)
         advanced_group.append(self.cleanup_threshold_widget)
 
+        self.fill_holes_widget = extraction_gui.fill_holes
+        self.fill_holes_widget.label = "Fill holes in segmentation"
+
+        self.closing_iterations_widget = extraction_gui.closing_iterations
+        self.closing_iterations_widget.label = "Closing iterations"
+
+        self.max_hole_size_widget = extraction_gui.max_hole_size
+        self.max_hole_size_widget.label = "Max hole size (pixels)"
+
+        self.show_preprocessed_widget = extraction_gui.show_preprocessed
+        self.show_preprocessed_widget.label = "Show preprocessed binary layer"
+
+        advanced_group.append(self.fill_holes_widget)
+        advanced_group.append(self.closing_iterations_widget)
+        advanced_group.append(self.max_hole_size_widget)
+        advanced_group.append(self.show_preprocessed_widget)
+
         # ============================================================
         # Output Settings (CLI file export options)
         # ============================================================
@@ -202,6 +241,10 @@ def _get_current_pipeline_config(self) -> PipelineConfig:
                 vessel_radius=self.include_vessel_radius_widget.value,
                 junction_cleanup=junction_cleanup,
                 cleanup_threshold_factor=self.cleanup_threshold_widget.value,
+                fill_holes=self.fill_holes_widget.value,
+                closing_iterations=self.closing_iterations_widget.value,
+                max_hole_size=self.max_hole_size_widget.value,
+                show_preprocessed=self.show_preprocessed_widget.value,
             ),
             output=OutputConfig(
                 write_skeleton_npy=self.write_skeleton_npy_widget.value,
@@ -221,6 +264,10 @@ def _set_pipeline_config(self, config: PipelineConfig) -> None:
         self.include_vessel_radius_widget.value = e.vessel_radius
         self.junction_cleanup_widget.value = e.junction_cleanup
         self.cleanup_threshold_widget.value = e.cleanup_threshold_factor
+        self.fill_holes_widget.value = e.fill_holes
+        self.closing_iterations_widget.value = e.closing_iterations
+        self.max_hole_size_widget.value = e.max_hole_size
+        self.show_preprocessed_widget.value = e.show_preprocessed
 
         o = config.output
         self.write_skeleton_npy_widget.value = o.write_skeleton_npy
@@ -299,6 +346,19 @@ def _on_analyze(self) -> None:
                 )
             )
 
+            # -- optional: show preprocessed binary layer ----------------
+            if (
+                pipeline_config.extraction.show_preprocessed
+                and result.preprocessed_binary is not None
+            ):
+                self.viewer.add_layer(
+                    Layer.create(
+                        result.preprocessed_binary,
+                        {"name": f"{img.name}_preprocessed"},
+                        "labels",
+                    )
+                )
+
             for data, meta, layer_type in result.layers:
                 try:
                     layer = Layer.create(data, meta, layer_type)

vesskel/config.py

@@ -31,6 +31,10 @@ class ExtractionConfig:
     vessel_radius: bool = False
     junction_cleanup: bool = False
     cleanup_threshold_factor: float = 2.5
+    closing_iterations: int = 0
+    fill_holes: bool = False
+    max_hole_size: int = 0
+    show_preprocessed: bool = False
 
     def to_dict(self) -> dict[str, Any]:
         return {
@@ -41,6 +45,10 @@ def to_dict(self) -> dict[str, Any]:
             "vessel_radius": self.vessel_radius,
             "junction_cleanup": self.junction_cleanup,
             "cleanup_threshold_factor": self.cleanup_threshold_factor,
+            "closing_iterations": self.closing_iterations,
+            "fill_holes": self.fill_holes,
+            "max_hole_size": self.max_hole_size,
+            "show_preprocessed": self.show_preprocessed,
         }
 
     @classmethod
@@ -54,6 +62,10 @@ def from_dict(cls, data: dict[str, Any]) -> ExtractionConfig:
             vessel_radius=data.get("vessel_radius", False),
             junction_cleanup=data.get("junction_cleanup", False),
             cleanup_threshold_factor=data.get("cleanup_threshold_factor", 2.5),
+            closing_iterations=data.get("closing_iterations", 0),
+            fill_holes=data.get("fill_holes", False),
+            max_hole_size=data.get("max_hole_size", 0),
+            show_preprocessed=data.get("show_preprocessed", False),
         )
 
 

vesskel/pipeline.py

@@ -6,6 +6,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+from scipy import ndimage as ndi
 from skan import summarize
 
 from vesskel._utils import to_binary
@@ -34,6 +35,7 @@ class AnalysisResult:
     summary_features: dict[str, float]
     branch_records: list[dict[str, object]]
     radius_matrix: np.ndarray | None = None
+    preprocessed_binary: np.ndarray | None = None
 
 
 def analyze_binary_image(
@@ -53,6 +55,33 @@ def analyze_binary_image(
         Full pipeline configuration.
     """
     binary = to_binary(image)
+
+    # -- optional: morphological preprocessing -------------------------
+    preprocessed_binary: np.ndarray | None = None
+    if config.extraction.closing_iterations > 0:
+        structure = ndi.generate_binary_structure(binary.ndim, 1)
+        binary = ndi.binary_closing(
+            binary, structure=structure, iterations=config.extraction.closing_iterations
+        ).astype(binary.dtype)
+
+    if config.extraction.fill_holes:
+        before_fill = binary.copy() if config.extraction.max_hole_size > 0 else None
+        binary = ndi.binary_fill_holes(binary).astype(binary.dtype)
+
+        if config.extraction.max_hole_size > 0 and before_fill is not None:
+            diff = binary.astype(np.int8) - before_fill.astype(np.int8)
+            filled = diff > 0
+            if filled.any():
+                labels, _ = ndi.label(filled)
+                sizes = np.bincount(labels.ravel())
+                big = sizes > config.extraction.max_hole_size
+                big[0] = False
+                revert = big[labels]
+                binary[revert] = 0
+
+    if config.extraction.closing_iterations > 0 or config.extraction.fill_holes:
+        preprocessed_binary = binary
+
     skeleton = lee94_thin(binary)
 
     if not skeleton.any():
@@ -129,4 +158,5 @@ def analyze_binary_image(
         summary_features=summary_features,
         branch_records=branch_records,
         radius_matrix=radius_matrix,
+        preprocessed_binary=preprocessed_binary,
     )