fix: build skeleton only once and pass through pipeline

Author: 404Simon

README.md

@@ -26,9 +26,9 @@ Inside the **Analyze Vessels** widget, tune extraction settings and use **Save C
 Use the same JSON preset exported from napari to batch-process images.
 
 ```sh
-uv run vesskel config-init --out config.json
-uv run vesskel validate-config --config config.json
-uv run vesskel run --input HRF/manual1 --config config.json --out outputs
+vesskel config-init --out config.json
+vesskel validate-config --config config.json
+vesskel run --input HRF/manual1 --config config.json --out outputs
 ```
 
 CLI outputs:

tests/test_2d_thinning_regression.py

@@ -10,8 +10,9 @@
 
 import numpy as np
 import pytest
+from skan import summarize
 
-from vesskel.features import extract_vessel_features
+from vesskel.features import build_vessel_graph, extract_vessel_features
 from vesskel.hrf import HRFDataset, preprocess_segmentation
 from vesskel.thin import lee94_thin
 
@@ -50,7 +51,9 @@ def test_skeleton_matches_baseline(self, dataset, index, request):
         info = dataset.image_list[index]
         name = info["name"]
         skeleton = _compute_skeleton(dataset, index)
-        features = extract_vessel_features(skeleton)
+        graph = build_vessel_graph(skeleton)
+        branch_data = summarize(graph, separator="-")
+        features = extract_vessel_features(skeleton, graph, branch_data)
         baseline_file = skeleton_path(name)
         feature_file = feature_path(name)
 

tests/test_3d_thinning_regression.py

@@ -2,9 +2,10 @@
 
 import numpy as np
 import pytest
+from skan import summarize
 from skimage import data
 
-from vesskel.features import extract_vessel_features
+from vesskel.features import build_vessel_graph, extract_vessel_features
 from vesskel.thin import lee94_thin
 
 from ._helpers import (
@@ -31,7 +32,9 @@ def image(self):
 
     def test_skeleton_matches_baseline(self, image, request):
         skeleton = _compute_skeleton(image)
-        features = extract_vessel_features(skeleton)
+        graph = build_vessel_graph(skeleton)
+        branch_data = summarize(graph, separator="-")
+        features = extract_vessel_features(skeleton, graph, branch_data)
         name = "brain"
         baseline_file = skeleton_path(name)
         feature_file = feature_path(name)

vesskel/cli.py

@@ -164,10 +164,12 @@ def _run_batch(args: argparse.Namespace) -> int:
     out_dir.mkdir(parents=True, exist_ok=True)
 
     summary_rows: list[dict[str, object]] = []
-
     output_name_counts: dict[str, int] = {}
+    total = len(input_paths)
+
+    for idx, in_path in enumerate(input_paths, 1):
+        print(f"[{idx}/{total}] {in_path.name}", flush=True)
 
-    for in_path in input_paths:
         image = _load_image(in_path)
         result = analyze_binary_image(
             image=image, base_name=in_path.stem, config=config

vesskel/extraction.py

@@ -3,14 +3,10 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from skan import summarize
+from skan import Skeleton
 
 from vesskel.config import ExtractionConfig
-from vesskel.features import (
-    build_vessel_graph,
-    compute_tortuosity,
-    extract_vessel_features,
-)
+from vesskel.features import compute_tortuosity
 
 if TYPE_CHECKING:
     from napari.types import LayerDataTuple
@@ -19,6 +15,8 @@
 def extract_skeleton_layers(
     skeleton: np.ndarray,
     base_name: str,
+    graph: Skeleton,
+    branch_data,
     config: ExtractionConfig | None = None,
     features: dict[str, float] | None = None,
 ) -> list["napari.types.LayerDataTuple"]:
@@ -30,19 +28,22 @@ def extract_skeleton_layers(
         Binary 2D or 3D skeleton array.
     base_name : str
         Base name for layer naming.
+    graph : Skeleton
+        Pre-built skan Skeleton graph (e.g. from `build_vessel_graph`).
+    branch_data : DataFrame
+        Pre-computed branch summary (e.g. from `skan.summarize(graph, ...)`).
     config : ExtractionConfig, optional
         Configuration for what to extract. Defaults to all except fractal_dimension.
     features : dict, optional
-        Pre-computed feature dictionary to avoid recomputation when caller
-        already has it (e.g. from extract_vessel_features).
+        Pre-computed summary feature dictionary (e.g. from `extract_vessel_features`).
     """
     if config is None:
         config = ExtractionConfig()
 
     layers = []
 
     if config.branches:
-        branch_layer = _extract_branch_features_layer(skeleton, base_name)
+        branch_layer = _extract_branch_features_layer(base_name, graph, branch_data)
         if branch_layer is not None:
             layers.append(branch_layer)
 
@@ -54,7 +55,6 @@ def extract_skeleton_layers(
         summary_layer = _extract_summary_features_layer(
             skeleton,
             base_name,
-            include_fractal=config.fractal_dimension,
             features=features,
         )
         layers.append(summary_layer)
@@ -63,33 +63,40 @@ def extract_skeleton_layers(
 
 
 def _extract_branch_features_layer(
-    skeleton: np.ndarray,
     base_name: str,
+    graph: Skeleton,
+    branch_data,
 ) -> "napari.types.LayerDataTuple | None":
     """Extract branch features and generate paths layer.
 
-    Returns None if skeleton has no branches.
+    Parameters
+    ----------
+    base_name : str
+        Base name used for layer naming.
+    graph : Skeleton
+        Pre-built skan Skeleton graph.
+    branch_data : DataFrame
+        Pre-computed branch summary from `skan.summarize`.
+
+    Returns
+    -------
+    LayerDataTuple or None
+        Napari shapes layer for branch paths, or None if skeleton has no branches.
     """
-    graph = build_vessel_graph(skeleton)
-    branch_data = summarize(graph, separator="-")
-
     if branch_data.empty:
         return None
 
     branch_data = branch_data.reset_index(drop=True).copy()
     branch_data["branch_id"] = np.arange(len(branch_data), dtype=np.int64)
 
-    # Compute tortuosity
     euclidean = branch_data["euclidean-distance"].to_numpy(dtype=float)
     branch_len = branch_data["branch-distance"].to_numpy(dtype=float)
     tortuosity = compute_tortuosity(branch_len, euclidean)
     tortuosity = np.nan_to_num(tortuosity, nan=1.0)
     branch_data["tortuosity"] = tortuosity
 
-    # Get branch path coordinates
     path_data = [graph.path_coordinates(i) for i in range(len(branch_data))]
 
-    # Determine if tortuosity varies significantly
     finite_tortuosity = tortuosity[np.isfinite(tortuosity)]
     varied_tortuosity = finite_tortuosity.size > 0 and float(
         np.min(finite_tortuosity)
@@ -120,11 +127,9 @@ def _extract_branch_text_layer(
     branch_layer: "napari.types.LayerDataTuple",
     base_name: str,
 ) -> "napari.types.LayerDataTuple":
-    """Create text labels for branches."""
     path_data = branch_layer[0]
     branch_data = branch_layer[1]["properties"]
 
-    # Compute label positions as mean of each path
     label_points = []
     for coords in path_data:
         if len(coords) == 0:
@@ -154,24 +159,22 @@ def _extract_branch_text_layer(
 def _extract_summary_features_layer(
     skeleton: np.ndarray,
     base_name: str,
-    include_fractal: bool = False,
-    features: dict[str, float] | None = None,
+    features: dict[str, float],
 ) -> "napari.types.LayerDataTuple":
-    """Extract global skeleton features and create summary point layer.
+    """Create a summary point layer displaying global skeleton features.
 
     Parameters
     ----------
-    features : dict, optional
-        Pre-computed feature dictionary. If provided, skips computation.
+    skeleton : ndarray
+        Binary 2D or 3D skeleton array. Used to position the summary label.
+    base_name : str
+        Base name for layer naming.
+    features : dict[str, float]
+        Pre-computed summary feature dictionary
+        (e.g. from `extract_vessel_features`).
     """
-    if features is None:
-        features = extract_vessel_features(
-            skeleton,
-            include_fractal=include_fractal,
-        )
     meta_features = {k: [v] for k, v in features.items()}
 
-    # Find center of foreground
     fg = np.argwhere(skeleton > 0)
     if fg.size:
         center = fg.mean(axis=0, dtype=float)

vesskel/features.py

@@ -2,7 +2,7 @@
 from numba import njit, prange
 from scipy.sparse import coo_matrix
 from scipy.sparse.csgraph import connected_components
-from skan import Skeleton, summarize
+from skan import Skeleton
 
 from ._utils import to_binary
 
@@ -122,16 +122,6 @@ def build_vessel_graph(skeleton: np.ndarray) -> Skeleton:
     return Skeleton(to_binary(skeleton))
 
 
-def summarize_skeleton(skeleton: np.ndarray):
-    """Build graph representation and return summarized branch data."""
-    if not np.any(skeleton):
-        from pandas import DataFrame
-
-        return DataFrame()
-    graph = build_vessel_graph(skeleton)
-    return summarize(graph, separator="-")
-
-
 _EMPTY_FEATURES: dict[str, float] = {
     "num_nodes": 0.0,
     "num_edges": 0.0,
@@ -154,19 +144,25 @@ def summarize_skeleton(skeleton: np.ndarray):
 
 def extract_vessel_features(
     skeleton: np.ndarray,
+    graph: Skeleton,
+    branch_data,
+    *,
     include_fractal: bool = True,
 ) -> dict[str, float]:
     """Extract graph-topology and segment statistics from a vessel skeleton.
 
     Parameters
     ----------
     skeleton : ndarray
-        Binary 2D or 3D skeleton array.
+        Binary 2D or 3D skeleton array. Used only for fractal dimension
+        computation, not for graph topology.
+    graph : Skeleton
+        Pre-built skan Skeleton graph (e.g. from `build_vessel_graph`).
+    branch_data : DataFrame
+        Pre-computed branch summary (e.g. from `skan.summarize(graph, ...)`).
     include_fractal : bool, optional
         Whether to compute fractal dimension (expensive-ish). Default is True.
     """
-    graph = build_vessel_graph(skeleton)
-    branch_data = summarize(graph, separator="-")
     if branch_data.empty:
         return dict(_EMPTY_FEATURES)
 

vesskel/pipeline.py

@@ -6,11 +6,12 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+from skan import summarize
 
 from vesskel._utils import to_binary
 from vesskel.config import PipelineConfig
 from vesskel.extraction import extract_skeleton_layers
-from vesskel.features import extract_vessel_features, summarize_skeleton
+from vesskel.features import build_vessel_graph, extract_vessel_features
 from vesskel.thin import lee94_thin
 
 if TYPE_CHECKING:
@@ -54,26 +55,30 @@ def analyze_binary_image(
             branch_records=[],
         )
 
+    graph = build_vessel_graph(skeleton)
+    branch_data = summarize(graph, separator="-")
+
     summary_features: dict[str, float] = {}
     if config.extraction.summary:
         summary_features = extract_vessel_features(
             skeleton,
+            graph,
+            branch_data,
             include_fractal=config.extraction.fractal_dimension,
         )
 
     layers = extract_skeleton_layers(
         skeleton,
         base_name,
-        config.extraction,
+        graph=graph,
+        branch_data=branch_data,
+        config=config.extraction,
         features=summary_features if config.extraction.summary else None,
     )
 
     branch_records: list[dict[str, object]] = []
-    if config.extraction.branches:
-        branch_table = summarize_skeleton(skeleton)
-        if not branch_table.empty:
-            # Keep native column names for traceability with skan output.
-            branch_records = branch_table.to_dict(orient="records")
+    if config.extraction.branches and not branch_data.empty:
+        branch_records = branch_data.to_dict(orient="records")
 
     return AnalysisResult(
         skeleton=skeleton,