Implement multi-component skeletonization

Author: constantinpape

MIGRATION_GUIDE.md

@@ -2502,8 +2502,9 @@ Important differences:
 ## kimimaro / TEASAR skeletonization
 
 Kimimaro skeletonizes densely labelled volumes and splits labels into connected
-components internally. The first `bioimage-cpp` TEASAR entry point is narrower:
-it accepts one binary 3D object and returns its graph directly.
+components internally. `bioimage-cpp` exposes the same two useful input
+semantics explicitly: `teasar` treats all nonzero values as one binary class,
+while `teasar_labels` preserves integer semantic-label identities.
 
 Kimimaro:
 
@@ -2536,17 +2537,33 @@ vertices, edges, radii = bic.skeleton.teasar(
     pdrf_exponent=4,
     number_of_threads=4,
 )
+
+skeletons = bic.skeleton.teasar_labels(
+    labels,
+    background=0,
+    spacing=(2.0, 1.0, 1.0),
+    scale=1.5,
+    constant=0,
+    pdrf_scale=100000,
+    pdrf_exponent=4,
+    number_of_threads=4,
+)
+# {original_label: (vertices, edges, radii), ...}
 ```
 
 Important differences and current scope:
 
-- `mask` is binary (nonzero means foreground), must be three-dimensional, and
-  a nonempty mask must contain exactly one 26-connected component. Component
-  splitting and multi-label dispatch are intentionally deferred.
+- `teasar(mask)` treats every nonzero value as one foreground class. It
+  skeletonizes every 26-connected component independently and returns their
+  deterministic concatenation as one forest tuple.
+- `teasar_labels(labels)` accepts native-endian `uint8`, `uint16`, `uint32`,
+  `uint64`, `int32`, or `int64` arrays. It returns one dictionary entry per
+  original non-background label. Disconnected occurrences of one label form
+  one forest, while touching distinct labels remain separate.
 - Coordinates follow NumPy axis order. `vertices` is `float64` physical
   `(z, y, x)` coordinates, `edges` is `(E, 2)` `uint64`, and `radii` is
-  per-vertex physical distance-to-boundary in `float32`. The graph is a
-  deterministic tree; an empty mask returns typed empty arrays.
+  per-vertex physical distance-to-boundary in `float32`. Empty binary masks
+  return typed empty arrays; all-background labeled inputs return `{}`.
 - The implementation follows the core TEASAR loop: a padded exact Euclidean
   distance-to-boundary field, a deterministic two-sweep root, a physical
   Dijkstra distance-from-root field, penalized repeated Dijkstra paths, and
@@ -2557,21 +2574,23 @@ Important differences and current scope:
   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.
-- The C++ core remains dependency-free. `number_of_threads=1` is the default;
-  `0` uses hardware concurrency. The thread budget controls the exact distance
-  transform. Compact-foreground root sweeps and rails remain ordered on the
-  optimized heap: threshold benchmarks showed that staged delta stepping nearly
-  doubled runtime for large solid objects. Compact IDs avoid full-volume
-  shortest-path fields. The public Dijkstra functions remain dense and exact
-  FP64, with delta stepping reserved for broad physical multi-source fields.
-
-Correctness tests are under `tests/skeleton/test_teasar.py`. The independent
+- 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
+  component fan-out and each component's exact distance transform without
+  nested oversubscription.
+- Compact-foreground root sweeps and rails remain ordered on the optimized
+  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
 Dijkstra benchmark is `development/distance/benchmark_dijkstra.py`; the
-end-to-end synthetic branching-tube benchmark is
-`development/skeleton/benchmark_teasar.py` and can optionally add kimimaro with
-`--kimimaro` when it is installed. Pass `--sequential-backends` to reproduce
-the dense/compact FP64 design matrix; extended density, spacing, and PDRF
-regimes are in `development/skeleton/benchmark_teasar_sequential.py`.
+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`.
 
 ## I/O and Build Dependencies
 

README.md

@@ -10,7 +10,8 @@ Image processing and segmentation functionality in C++ with light-weight python
 The package includes dependency-free triangle-mesh extraction from 3D volumes
 and segmentation masks, plus Laplacian mesh smoothing, under `bioimage_cpp.mesh`.
 It also includes exact masked-grid Dijkstra paths under `bioimage_cpp.distance`
-and binary 3D TEASAR skeletonization under `bioimage_cpp.skeleton`.
+and binary-forest plus semantic multi-label 3D TEASAR skeletonization under
+`bioimage_cpp.skeleton`.
 
 The `bioimage_cpp` python library can be installed via pip:
 ```bash

development/skeleton/PERFORMANCE_NOTES.md

@@ -513,3 +513,95 @@ The following proposed optimizations were deliberately deferred:
   calculation was shared.
 - Candidate-target and state-reset changes are therefore benchmark follow-ups,
   not pending correctness fixes.
+
+## Multi-component and multi-label dispatch
+
+The run-based dispatcher and `teasar_labels` implementation were measured on
+2026-07-15 on the same 4-core / 8-hardware-thread Intel Core i7-1185G7 host.
+The normal optimized build used Python 3.13.13, NumPy 2.4.6, kimimaro 5.8.1,
+and EDT 3.1.1. Every row used one warmup, five measured calls, deterministic
+backend-order shuffling, spacing `(1.5, 1.0, 1.0)`, and the established TEASAR
+parameters. Kimimaro's soma, border, hole-filling, dust, and progress features
+were disabled as in the earlier comparisons.
+
+Reproduce the complete matrix and fresh-process memory probes with:
+
+```bash
+python development/skeleton/benchmark_teasar.py \
+    --large --suite all --kimimaro --threads 1 2 4 8 \
+    --repeats 5 --warmup 1 --memory --memory-threads 1 4 \
+    --json /tmp/teasar_multilabel_full.json
+```
+
+The table reports five-repeat medians. `loop` is the deliberately naive
+`np.unique` plus one full-volume `labels == L` and binary TEASAR call per label.
+
+| scenario | case | components / labels | bio t1 | bio t4 | loop t1 | kimi t1 | kimi t4 |
+| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: |
+| binary | branching tube `128³` | 1 / 1 | 35.45 ms | 23.52 ms | -- | 132.00 ms | 154.82 ms |
+| binary | branching tube `192³` | 1 / 1 | 109.90 ms | 71.65 ms | -- | 350.35 ms | 370.80 ms |
+| binary | branching tube `256³` | 1 / 1 | 272.59 ms | 179.95 ms | -- | 903.28 ms | 846.81 ms |
+| binary | packed tubes | 27 / 1 | 52.84 ms | 17.70 ms | -- | 203.60 ms | 177.64 ms |
+| labels | packed distinct | 27 / 27 | 53.20 ms | 18.47 ms | 264.04 ms | 206.36 ms | 199.80 ms |
+| labels | packed fragmented | 27 / 1 | 53.34 ms | 17.86 ms | 108.68 ms | 202.99 ms | 201.72 ms |
+| labels | imbalanced | 65 / 65 | 39.65 ms | 40.29 ms | 749.70 ms | 309.99 ms | 347.58 ms |
+| labels | dense balls | 8 / 8 | 49.01 ms | 14.17 ms | 69.15 ms | 166.97 ms | 130.41 ms |
+
+The packed binary and 27-label inputs contain exactly the same 33,426
+foreground voxels and produced array-exact forests after canonical grouping.
+Their t1 medians differ by less than 1%, so preserving semantic labels adds no
+material cost beyond equality-aware run scanning. The balanced label case
+scales 2.88x from one to four workers and the dense case scales 3.46x. The
+imbalanced case does not scale because one large component dominates the 64
+three-voxel tasks; the dynamic cursor avoids a regression through four workers,
+but eight workers add scheduling overhead.
+
+Native dispatch is 4.96x faster than the repeated-volume loop on the packed
+27-label case and 18.9x faster on the 65-label imbalanced case at one worker.
+The connected binary medians are also below the pre-dispatch measurements
+recorded above, clearing the no-regression gate. Bioimage-cpp is 3.4--7.8x
+faster than kimimaro at one worker across the headline multi-label cases and
+8.6--11.3x faster at four workers.
+
+### Worst-case runs and memory
+
+The explicit run-count stress case is a `48³` alternating volume: 110,592
+voxels, 110,592 x-runs, two semantic labels, and two 26-connected components.
+Three-repeat medians were 750.25 ms (t1) and 381.23 ms (t4) for bioimage-cpp,
+versus 1.673 s and 937.81 ms for kimimaro. It completed without a dense
+component-label fallback.
+
+Linux memory probes start each backend in a fresh worker and sample aggregate
+`/proc` RSS for the worker and all descendants; the worker's own `ru_maxrss`
+delta covers native calls shorter than the sampling interval. Incremental
+process-tree peaks were:
+
+| case | bio t1 / t4 | kimimaro t1 / t4 |
+| --- | ---: | ---: |
+| binary tube `256³` | 64.2 / 61.8 MiB | 182.0 / 751.7 MiB |
+| packed distinct labels | 0.0 / 2.7 MiB | 37.9 / 655.8 MiB |
+| packed fragmented label | 0.2 / 2.8 MiB | 37.2 / 654.7 MiB |
+| imbalanced labels | 5.6 / 5.6 MiB | 61.0 / 765.9 MiB |
+| alternating labels | 3.7 / 7.6 MiB | 28.4 / 429.5 MiB |
+
+A zero incremental value means the short native call stayed below the
+post-import/input high-water mark; it does not mean zero allocations. The
+important structural check is also enforced in code: component discovery owns
+only run metadata, row offsets, union-find state, and concurrently active tight
+crops, with no full-volume `uint32`/`uint64` component image.
+
+### Dispatcher profile
+
+A profile build on the 27-label packed case at one worker reported:
+
+| phase | time | share |
+| --- | ---: | ---: |
+| component scan | 4.6 ms | 8.7% |
+| run union | 0.2 ms | 0.4% |
+| descriptors | 0.2 ms | 0.4% |
+| component TEASAR | 48.0 ms | 90.4% |
+| forest assembly | <0.1 ms | <0.1% |
+
+The normal optimized build was restored after profiling. Skeleton outputs were
+exact across worker counts `1, 2, 4, 8`, and the benchmark validates `E = V-C`,
+semantic keys, paired binary/label parity, and raw samples before writing JSON.