Implement primitives for blockwise skeletonization

CMakeLists.txt

@@ -23,6 +23,7 @@ nanobind_add_module(_core
     src/bindings/mesh.cxx
     src/bindings/segmentation.cxx
     src/bindings/skeleton.cxx
+    src/bindings/skeleton_distributed.cxx
     src/bindings/transformation.cxx
     src/bindings/util.cxx
     src/bindings/utils.cxx

MIGRATION_GUIDE.md

@@ -2589,12 +2589,14 @@ Important differences and current scope:
   distance-to-boundary field, a deterministic two-sweep root, a physical
   Dijkstra distance-from-root field, penalized repeated Dijkstra paths, and
   rolling invalidation with radius `scale * radius + constant`.
-- This first correctness-oriented version uses a physical axis-aligned
-  invalidation cube. It does not implement kimimaro's soma handling, hole
-  filling, border stitching, manual targets, component dust filtering,
-  cross-section metadata, or postprocessing heuristics. These differences can
-  change branch positions and vertex counts, so output is not expected to be
-  vertex-for-vertex identical to kimimaro.
+- This correctness-oriented implementation uses a physical axis-aligned
+  invalidation cube. The ordinary in-core entry points do not automatically
+  perform block stitching or accept manual targets; the dedicated block APIs
+  below provide mandatory interface targets and exact consolidation. Soma
+  handling, hole filling, component dust filtering, cross-section metadata,
+  and kimimaro's other postprocessing heuristics are not implemented. These
+  differences can change branch positions and vertex counts, so output is not
+  expected to be vertex-for-vertex identical to kimimaro.
 - The C++ core remains dependency-free. Component discovery uses x-runs and a
   union-find rather than dense component-label images. `number_of_threads=1`
   is the default and `0` uses hardware concurrency; one shared budget covers
@@ -2604,14 +2606,76 @@ Important differences and current scope:
   heap. Compact IDs avoid full-volume shortest-path fields; the public Dijkstra
   functions remain dense and exact FP64.
 
-Correctness tests are under `tests/skeleton/test_teasar.py` and
-`tests/skeleton/test_teasar_labels.py`. The independent
+### Blockwise skeletonization and exact stitching
+
+`bioimage_cpp.skeleton.distributed` provides the dependency-light primitives
+needed by an external block scheduler. It does not perform I/O, enumerate
+blocks, launch tasks, or serialize artifacts. The implemented layout uses one
+shared voxel plane: a core block with interval `[a, b)` is processed through
+`b + 1` when it has a high-side neighbor, while the neighbor begins at `b`.
+
+```python
+dist = bic.skeleton.distributed
+
+# `left` includes its high-side overlap plane and starts at `left_origin`.
+targets = dist.block_border_targets(
+    left,
+    faces=[(2, "high")],
+    origin=left_origin,
+    spacing=spacing,
+)
+left_fragment = dist.block_teasar(
+    left,
+    open_faces=[(2, "high")],
+    origin=left_origin,
+    required_targets=targets,
+    spacing=spacing,
+)
+
+# Repeat independently for every processing block. Neighboring calls select
+# identical global targets on their shared plane.
+merged = dist.merge_block_skeletons(block_fragments)
+forest = dist.minimum_spanning_forest(merged, spacing=spacing)  # optional
+vertices, edges, radii = dist.lattice_to_physical(forest, spacing=spacing)
+```
+
+Block fragments use global `int64` lattice vertices, `uint64` edges, and
+`float32` physical radii. Integer coordinates are the exact merge identity;
+physical proximity is never used to invent a connection. Duplicate radii are
+reduced with `maximum`, and canonical sort-and-unique output makes
+`merge_block_skeletons` associative, commutative, idempotent, and suitable for
+hierarchical reduction. `merge_block_skeleton_maps` applies the same operation
+per exact semantic label, including negative and large `uint64` labels.
+
+`open_faces` is mandatory even when empty. It identifies artificial cuts where
+the distance transform may look through the processing-block boundary; paths
+remain restricted to real input foreground. This prevents interface radii from
+collapsing to one voxel. When present, the lexicographically greatest required
+target on an open face becomes the deterministic component root.
+
+The labeled counterparts are `block_border_targets_labels` and
+`block_teasar_labels`. Border targets are selected per 8-connected same-label
+face patch by maximum anisotropic 2D distance transform with deterministic,
+orientation-stable plateau resolution. Several valid targets can create cycles
+after graph union, so exact consolidation keeps cycles and
+`minimum_spanning_forest` removes them only when explicitly called. Blocked
+skeletons are connected through their anchors but are not expected to be
+vertex-for-vertex identical to whole-volume TEASAR because each block still has
+less path-selection context.
+
+Correctness tests are under `tests/skeleton/test_teasar.py`,
+`tests/skeleton/test_teasar_labels.py`, and
+`tests/skeleton/test_distributed.py`. The independent
 Dijkstra benchmark is `development/distance/benchmark_dijkstra.py`; the
 end-to-end benchmark is `development/skeleton/benchmark_teasar.py`. Use
 `--suite all --kimimaro` to compare paired binary and multi-label scenarios,
 `--memory` for fresh-process peak-RSS probes, or `--sequential-backends` for
 the dense/compact FP64 binary design matrix. Extended density, spacing, and
 PDRF regimes are in `development/skeleton/benchmark_teasar_sequential.py`.
+The development-only serial block harness is
+`development/skeleton/blockwise_stitching.py`; the rigorous one-thread layered
+comparison against whole and blocked Kimimaro is
+`development/skeleton/compare_blockwise_kimimaro.py`.
 
 ## I/O and Build Dependencies
 

development/skeleton/PERFORMANCE_NOTES.md

@@ -736,3 +736,58 @@ leave the initial linear mode. Skeleton output remained array-exact across
 backends and worker counts. Final verification after the optimization and
 benchmark addition: `1182 passed` with third-party pytest plugin autoload
 disabled.
+
+## Blockwise open-boundary correctness and Kimimaro comparison (2026-07-16)
+
+The first block prototype reused the ordinary per-component zero crop halo for
+its DBF. On an artificial cut this made every interface radius exactly one
+voxel. Thick branching tubes remained connected, but accumulated 2.1--3.1x the
+whole-volume Kimimaro cable length. A controlled Kimimaro run reproduced the
+failure when `black_border=True`, isolating the problem from target selection,
+exact consolidation, and cycle removal.
+
+The corrected block path keeps its real foreground path mask unchanged and
+uses a separate EDT-only mask with one ghost layer across explicitly declared
+`open_faces`. One deterministic face target roots each affected component.
+Ghost voxels can never become graph vertices. The complete reference workflow
+is:
+
+```bash
+python development/skeleton/compare_blockwise_kimimaro.py \
+    --warmup 1 --repeats 7 --check \
+    --json /tmp/blockwise-kimimaro.json
+```
+
+The run used bioimage-cpp 0.7.0, Kimimaro 5.8.1, EDT 3.1.1, NumPy 2.4.6,
+SciPy 1.17.1, isotropic spacing, and one thread for every backend. Kimimaro's
+blocked fragments used `fix_borders=True`; both block backends then used the
+same exact lattice merge and minimum-spanning forest.
+
+| case | bio block | Kimimaro whole | Kimimaro block | bio / Kimimaro-block cable | direct p95 / Hausdorff |
+| --- | ---: | ---: | ---: | ---: | ---: |
+| thin line, 16 blocks | 3.52 ms | 1.66 ms | 6.95 ms | 1.000 | 0.00 / 0.00 |
+| oblique tube `64^3`, anisotropic, 8 blocks | 4.41 ms | 11.44 ms | 17.78 ms | 1.000 | 0.00 / 0.00 |
+| tube `64^3`, 8 blocks | 6.73 ms | 18.34 ms | 46.97 ms | 0.988 | 0.20 / 1.41 |
+| tube `96^3`, 8 blocks | 13.04 ms | 48.95 ms | 72.01 ms | 0.916 | 1.41 / 3.74 |
+| tube `128^3`, seam-aligned, 8 blocks | 21.70 ms | 102.29 ms | 123.76 ms | 0.988 | 0.00 / 4.24 |
+| tube `128^3`, 27 blocks | 24.85 ms | 102.02 ms | 118.78 ms | 0.903 | 1.27 / 5.00 |
+
+The thin line remains graph-exact. The anisotropic oblique case has identical
+blocked vertex geometry and cable length. All 112 expected local interface
+anchors were present in both controlled implementations, and their radii
+matched Kimimaro exactly. Across all thick cases, final graphs were connected
+and acyclic, cable length stayed within 10% of Kimimaro's blocked result, and
+every quality gate passed. On the nontrivial timing cases, bioimage-cpp was
+2.6--4.7x faster than whole-volume one-thread Kimimaro and 4.0--7.0x faster than
+its serial block workflow.
+
+On the final `128^3` calls, the last measured phase breakdowns were:
+
+| layout | targets | local TEASAR | merge | forest |
+| --- | ---: | ---: | ---: | ---: |
+| 8 blocks | 2.39 ms | 17.91 ms | 0.25 ms | 0.09 ms |
+| 27 blocks | 4.68 ms | 15.89 ms | 0.35 ms | 0.10 ms |
+
+Merge and cycle removal remain negligible. More blocks mainly increase face
+analysis and Python orchestration. The comparison intentionally does not apply
+Kimimaro's heuristic nearest-component joining, dust removal, or tick pruning.

development/skeleton/blockwise_stitching.py

@@ -0,0 +1,173 @@
+"""Minimal serial harness for blockwise TEASAR and exact stitching.
+
+This is deliberately development-only orchestration. It performs no I/O,
+parallel scheduling, retries, or artifact serialization; those responsibilities
+belong to downstream block-processing systems.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+
+import numpy as np
+
+import bioimage_cpp as bic
+
+
+dist = bic.skeleton.distributed
+
+
+def _processing_blocks(volume: np.ndarray, block_shape: Sequence[int]):
+    blocking = bic.utils.Blocking([0, 0, 0], list(volume.shape), list(block_shape))
+    for block_id in range(blocking.number_of_blocks):
+        core = blocking.get_block(block_id)
+        begin = [int(value) for value in core.begin]
+        end = [int(value) for value in core.end]
+        faces = []
+        for axis in range(3):
+            if blocking.get_neighbor_id(block_id, axis, lower=True) >= 0:
+                faces.append((axis, "low"))
+            if blocking.get_neighbor_id(block_id, axis, lower=False) >= 0:
+                faces.append((axis, "high"))
+                end[axis] += 1
+        slices = tuple(slice(begin[axis], end[axis]) for axis in range(3))
+        yield blocking, block_id, np.ascontiguousarray(volume[slices]), begin, faces
+
+
+def _unique_coordinates(parts: list[np.ndarray]) -> np.ndarray:
+    if not parts:
+        return np.empty((0, 3), dtype=np.int64)
+    return np.unique(np.concatenate(parts, axis=0), axis=0)
+
+
+def _unique_target_map(parts: list[dict[int, np.ndarray]]) -> dict[int, np.ndarray]:
+    labels = sorted({label for part in parts for label in part})
+    return {
+        label: _unique_coordinates([part[label] for part in parts if label in part])
+        for label in labels
+    }
+
+
+def _assert_matching_interfaces(blocking, per_face) -> None:
+    for block_id in range(blocking.number_of_blocks):
+        for axis in range(3):
+            neighbor = blocking.get_neighbor_id(block_id, axis, lower=False)
+            if neighbor < 0:
+                continue
+            left = per_face[(block_id, axis, "high")]
+            right = per_face[(neighbor, axis, "low")]
+            if isinstance(left, dict):
+                if left.keys() != right.keys():
+                    raise AssertionError(
+                        f"target labels disagree across blocks {block_id}/{neighbor}"
+                    )
+                for label in left:
+                    np.testing.assert_array_equal(left[label], right[label])
+            else:
+                np.testing.assert_array_equal(left, right)
+
+
+def run_blockwise_binary(
+    mask: np.ndarray,
+    block_shape: Sequence[int],
+    *,
+    spacing=(1.0, 1.0, 1.0),
+    remove_cycles: bool = False,
+):
+    """Run the complete binary block pipeline serially and return a lattice graph."""
+    fragments = []
+    per_face = {}
+    blocking = None
+    for blocking, block_id, block, origin, faces in _processing_blocks(
+        np.asarray(mask), block_shape
+    ):
+        face_targets = []
+        for axis, side in faces:
+            targets = dist.block_border_targets(
+                block,
+                [(axis, side)],
+                origin=origin,
+                spacing=spacing,
+                number_of_threads=1,
+            )
+            per_face[(block_id, axis, side)] = targets
+            face_targets.append(targets)
+        targets = _unique_coordinates(face_targets)
+        fragments.append(
+            dist.block_teasar(
+                block,
+                open_faces=faces,
+                origin=origin,
+                required_targets=targets,
+                spacing=spacing,
+                number_of_threads=1,
+            )
+        )
+    if blocking is None:
+        raise ValueError("mask must be a three-dimensional array")
+    _assert_matching_interfaces(blocking, per_face)
+    merged = dist.merge_block_skeletons(fragments)
+    if remove_cycles:
+        merged = dist.minimum_spanning_forest(merged, spacing=spacing)
+    return merged
+
+
+def run_blockwise_labels(
+    labels: np.ndarray,
+    block_shape: Sequence[int],
+    *,
+    background: int = 0,
+    spacing=(1.0, 1.0, 1.0),
+    remove_cycles: bool = False,
+):
+    """Run the labeled block pipeline serially and return lattice graphs by label."""
+    fragment_maps = []
+    per_face = {}
+    blocking = None
+    for blocking, block_id, block, origin, faces in _processing_blocks(
+        np.asarray(labels), block_shape
+    ):
+        face_targets = []
+        for axis, side in faces:
+            targets = dist.block_border_targets_labels(
+                block,
+                [(axis, side)],
+                origin=origin,
+                background=background,
+                spacing=spacing,
+                number_of_threads=1,
+            )
+            per_face[(block_id, axis, side)] = targets
+            face_targets.append(targets)
+        targets = _unique_target_map(face_targets)
+        fragment_maps.append(
+            dist.block_teasar_labels(
+                block,
+                open_faces=faces,
+                origin=origin,
+                required_targets=targets,
+                background=background,
+                spacing=spacing,
+                number_of_threads=1,
+            )
+        )
+    if blocking is None:
+        raise ValueError("labels must be a three-dimensional array")
+    _assert_matching_interfaces(blocking, per_face)
+    merged = dist.merge_block_skeleton_maps(fragment_maps)
+    if remove_cycles:
+        merged = {
+            label: dist.minimum_spanning_forest(fragment, spacing=spacing)
+            for label, fragment in merged.items()
+        }
+    return merged
+
+
+if __name__ == "__main__":
+    example = np.zeros((17, 17, 25), dtype=np.uint8)
+    example[8, 8, 2:23] = 1
+    result = run_blockwise_binary(example, (9, 9, 8), remove_cycles=True)
+    print(
+        f"stitched {len(result[0])} lattice vertices and {len(result[1])} edges "
+        "from serial processing blocks"
+    )