Skip to content

Commit b142364

Browse files
constantinpapeconstantinpape
authored
Marching cubes codex sol (#60)
* Implement marching cubes * Improve implementation, tests, and benchmark
1 parent f1e0ae9 commit b142364

18 files changed

Lines changed: 3192 additions & 0 deletions

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ nanobind_add_module(_core
2020
src/bindings/graph.cxx
2121
src/bindings/ground_truth.cxx
2222
src/bindings/label_multiset.cxx
23+
src/bindings/mesh.cxx
2324
src/bindings/segmentation.cxx
2425
src/bindings/transformation.cxx
2526
src/bindings/util.cxx

‎LICENSE‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,3 +19,33 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
1919
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
2020
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
2121
SOFTWARE.
22+
23+
Third-party notice — scikit-image Marching Cubes tables
24+
--------------------------------------------------------
25+
26+
The source and binary distributions also include lookup-table data derived
27+
from scikit-image 0.26.0. Copyright (c) 2009-2022 the scikit-image team.
28+
All rights reserved.
29+
30+
Redistribution and use in source and binary forms, with or without
31+
modification, are permitted provided that the following conditions are met:
32+
33+
1. Redistributions of source code must retain the above copyright notice,
34+
this list of conditions and the following disclaimer.
35+
2. Redistributions in binary form must reproduce the above copyright notice,
36+
this list of conditions and the following disclaimer in the documentation
37+
and/or other materials provided with the distribution.
38+
3. Neither the name of the copyright holder nor the names of its contributors
39+
may be used to endorse or promote products derived from this software
40+
without specific prior written permission.
41+
42+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
43+
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
44+
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
45+
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
46+
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
47+
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
48+
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
49+
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
50+
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
51+
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

‎MIGRATION_GUIDE.md‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1872,6 +1872,60 @@ Notes:
18721872

18731873
## Skimage
18741874

1875+
### Marching Cubes
1876+
1877+
`bioimage-cpp` provides dependency-free isosurface extraction under
1878+
`bic.mesh`, including both the topology-resolving Lewiner/MC33 method used by
1879+
default in scikit-image and the classic Lorensen lookup-table variant.
1880+
1881+
```python
1882+
import bioimage_cpp as bic
1883+
1884+
# Extract one object from a label image. `pad=True` closes objects that touch
1885+
# the volume boundary by adding a temporary zero-valued halo.
1886+
vertices, faces, normals, values = bic.mesh.marching_cubes(
1887+
labels == label_id,
1888+
level=0.5,
1889+
spacing=(z_spacing, y_spacing, x_spacing),
1890+
method="lewiner",
1891+
pad=True,
1892+
)
1893+
```
1894+
1895+
The signature follows `skimage.measure.marching_cubes`: `level`, `spacing`,
1896+
`gradient_direction`, `step_size`, `allow_degenerate`, `method`, and an
1897+
optional boolean `mask` have the same purpose. Coordinates use NumPy
1898+
`(z, y, x)` order. Vertices are `float32` at unit spacing and `float64` after
1899+
non-unit spacing; faces are consistently `int32`, and normals/values are
1900+
`float32`.
1901+
1902+
Important details:
1903+
1904+
- Inputs are converted to contiguous `float32` before extraction. Any real
1905+
numeric or boolean input dtype is accepted; complex inputs are rejected.
1906+
- `method="lewiner"` resolves ambiguous cases and is the default;
1907+
`method="lorensen"` selects the original 256-case algorithm.
1908+
- `spacing` accepts either one positive finite scalar for isotropic data or a
1909+
length-three `(z, y, x)` sequence.
1910+
- Normals and local-range values follow scikit-image semantics. As in
1911+
scikit-image, `gradient_direction` reverses face winding without changing
1912+
normals, and anisotropic spacing scales vertices without transforming
1913+
normals.
1914+
- `pad=False` matches scikit-image's open-boundary behavior. The additional
1915+
`pad=True` option uses a zero-valued halo and is intended for
1916+
foreground-positive segmentation masks. The iso-level is determined from
1917+
the original unpadded volume.
1918+
- Spacing entries must be positive and finite, and faces remain `int32` when
1919+
degenerate faces are removed. Duplicate vertices in a collapsed face are
1920+
merged transitively with the first vertex as representative, and faces that
1921+
still collapse after remapping are discarded. This guarantees in-range face
1922+
indices and intentionally avoids a rare scikit-image negative-index
1923+
remapping quirk. These validation and cleanup choices are intentional
1924+
differences from scikit-image edge cases.
1925+
1926+
See `development/mesh/check_marching_cubes.py` for reference comparisons and
1927+
`development/mesh/benchmark_marching_cubes.py` for reproducible timings.
1928+
18751929
### Anti-Aliased Resampling
18761930

18771931
`affine_transform` itself never pre-smooths the input; downsampling without

‎README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@
77

88
Image processing and segmentation functionality in C++ with light-weight python bindings through nanobind and minimal dependencies to enable distribution via pip.
99

10+
The package includes dependency-free triangle-mesh extraction from 3D volumes
11+
and segmentation masks under `bioimage_cpp.mesh`.
12+
1013
The `bioimage_cpp` python library can be installed via pip:
1114
```bash
1215
pip install bioimage-cpp

‎THIRD_PARTY_NOTICES.md‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# Third-party notices
2+
3+
## scikit-image Marching Cubes tables
4+
5+
`include/bioimage_cpp/mesh/detail/mc33_luts.hxx` contains lookup-table data
6+
derived from scikit-image 0.26.0's Marching Cubes implementation. The
7+
corresponding MC33 control flow in `marching_cubes.hxx` follows that reference
8+
port of the algorithm by Lewiner et al. The relevant scikit-image material is
9+
licensed under the BSD 3-Clause License:
10+
11+
Copyright (c) 2009-2022 the scikit-image team. All rights reserved.
12+
13+
Redistribution and use in source and binary forms, with or without
14+
modification, are permitted provided that the following conditions are met:
15+
16+
1. Redistributions of source code must retain the above copyright notice,
17+
this list of conditions and the following disclaimer.
18+
2. Redistributions in binary form must reproduce the above copyright notice,
19+
this list of conditions and the following disclaimer in the documentation
20+
and/or other materials provided with the distribution.
21+
3. Neither the name of the copyright holder nor the names of its contributors
22+
may be used to endorse or promote products derived from this software
23+
without specific prior written permission.
24+
25+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26+
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27+
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28+
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29+
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30+
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31+
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32+
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33+
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34+
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

‎development/mesh/PERFORMANCE_NOTES.md‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
# Marching cubes performance
2+
3+
## Codex-Sol / Claude consolidation (2026-07-11)
4+
5+
The consolidation pass kept the Codex-Sol float32 MC33 kernel, exact
6+
scikit-image normal/value semantics, and zero-copy NumPy handoff. It moved
7+
axis conversion and winding into triangle emission, added robust transitive
8+
degenerate-vertex merging with the shared `UnionFind`, and broadened the
9+
correctness and benchmark suites. Scalar spacing is now accepted as an
10+
isotropic convenience.
11+
12+
The benchmark now has a backward-compatible single-size mode and a full
13+
scaling suite matching the independent implementation comparison:
14+
15+
```bash
16+
python development/mesh/benchmark_marching_cubes.py --size medium
17+
python development/mesh/benchmark_marching_cubes.py \
18+
--suite scaling --repeats 5 --batches 1 --warmup 1 --memory \
19+
--json /tmp/marching-cubes.json
20+
```
21+
22+
The scaling suite alternates paired bioimage-cpp/scikit-image calls and covers
23+
Lewiner spheres/scalar fields through 512³, dense 10%-foreground masks through
24+
256³, and all three Lorensen workloads at 128³. Fresh subprocesses measure
25+
peak RSS before the timing process allocates any large volumes.
26+
27+
Same-harness before/after results on the machine described below:
28+
29+
| case | before ms | after ms | change | scikit-image / after |
30+
|---|---:|---:|---:|---:|
31+
| sphere, Lewiner, 512³ | 1,898.68 | 1,786.24 | -5.9% | 1.06× |
32+
| dense mask, Lewiner, 256³ | 2,119.26 | 2,063.38 | -2.6% | 1.98× |
33+
| scalar field, Lewiner, 512³ | 2,268.58 | 2,295.50 | +1.2% | 1.30× |
34+
| sphere, Lorensen, 128³ | 30.31 | 28.39 | -6.3% | 1.31× |
35+
| dense mask, Lorensen, 128³ | 192.19 | 212.75 | +10.7% | 2.16× |
36+
| scalar field, Lorensen, 128³ | 48.72 | 49.80 | +2.2% | 1.84× |
37+
38+
Absolute timings moved with CPU state: across all 15 cases the geometric-mean
39+
wall time improved 2.4%, while normalization by each paired scikit-image time
40+
showed a 1.1% improvement. Targeted same-session A/B repeats against the exact
41+
pre-change commit resolved the apparent small-scalar and dense-Lorensen
42+
regressions: scalar 64³ retained a ~2.04–2.06× paired speedup, while dense
43+
Lorensen improved from ~2.01–2.02× to ~2.08×. No reproducible regression above
44+
3% remained.
45+
46+
Peak RSS was unchanged within measurement noise:
47+
48+
| case | before MiB | after MiB | scikit-image MiB |
49+
|---|---:|---:|---:|
50+
| sphere 512³ | 725.8 | 725.7 | 734.6 |
51+
| dense mask 256³ | 643.5 | 643.5 | 1,108.6 |
52+
| scalar field 512³ | 703.3 | 703.4 | 843.3 |
53+
54+
Correctness validation finished with 29 mesh tests, 1,024 full-suite tests,
55+
all 254 nontrivial binary cube configurations for both methods, 1,024 random
56+
scalar cubes for both methods, and 1,000 randomized multi-cube mask/stride/
57+
degeneracy cases. Twelve `allow_degenerate=False` cases intentionally differed
58+
from scikit-image's negative-index remapping quirk; all returned faces were
59+
valid and contained no collapsed-coordinate triangles.
60+
61+
## Two-slice face-cache optimization
62+
63+
`bic.mesh.marching_cubes` remains deterministic and single-threaded. The
64+
optimization pass replaced the volume-growing edge hash map with Lewiner's
65+
two-slice face cache. Several smaller candidates were measured first and
66+
rejected when their gains did not survive repeated benchmarks.
67+
68+
## Measurement setup
69+
70+
- CPU: Intel Core i7-1185G7, 4 cores / 8 threads; the kernel uses one thread.
71+
- OS: Linux 5.15 x86_64.
72+
- Compiler/build: conda-forge GCC 14.3, editable release build (`-O3`).
73+
- Python 3.13.13, NumPy 2.4.6, scikit-image 0.26.0.
74+
- Workloads: a binary sphere, a reproducible 10%-foreground random binary
75+
mask, and a deterministic smooth scalar field.
76+
- Medium timings are medians of three seven-call batch medians. Small timings
77+
use three batches of three calls; large timings use one batch of three calls.
78+
- Every timed case passes the geometry/topology reference comparison first.
79+
80+
Reproduce:
81+
82+
```bash
83+
python development/mesh/benchmark_marching_cubes.py --size small --repeats 3 --batches 3
84+
python development/mesh/benchmark_marching_cubes.py --size medium --repeats 7 --batches 3
85+
python development/mesh/benchmark_marching_cubes.py --size large --repeats 3 --batches 1
86+
```
87+
88+
## Candidate evaluation
89+
90+
Each candidate was rebuilt, checked against the mesh tests and scikit-image
91+
oracle, then timed twice at 96³. A candidate needed a repeatable improvement
92+
of at least 3% without a regression above 3%.
93+
94+
| candidate | observed result | decision |
95+
|---|---|---|
96+
| Slice-sized output reserves and scalar normal appends | First run helped spheres ~3%, repeat regressed spheres 3–4% and dense Lewiner 11% | Reverted |
97+
| Hoisted flat indexing plus one-pass hash insertion | Dense masks improved ~2.5%, but spheres repeatedly regressed 3–4% | Reverted |
98+
| Cached corner strengths and per-cell edge indices | Mostly 1–3% changes; dense Lorensen regressed ~2% | Reverted |
99+
| Two-slice face cache | Repeatable 17–22% sphere, 55% scalar, and 80–81% dense-mask reductions at 96³ | Retained |
100+
| Flat cube loads after the face cache | Helped sparse Lewiner, but repeat regressed scalar Lewiner 4% and left dense masks unchanged | Reverted |
101+
102+
The retained cache uses two `int32` arrays with four slots per `(y, x)`
103+
position. The upper edge layer becomes the next z-slice's lower layer; the new
104+
upper layer is cleared. MC33 center vertices remain cell-local. Deduplication
105+
memory is therefore `O(nx * ny)` instead of growing with the total surface.
106+
107+
## Profiling
108+
109+
Profiling used the repository's `BIOIMAGE_PROFILE` build on a 160³ Lewiner
110+
call. The baseline's global map destruction happened after the original core
111+
report, so it appears as the gap between core traversal/finalization and the
112+
binding's `core_call`.
113+
114+
| dense-mask phase | hash-map baseline | two-slice cache |
115+
|---|---:|---:|
116+
| cell traversal | 2.208 s | 0.449 s |
117+
| normal finalization | 0.009 s | 0.009 s |
118+
| cache/map cleanup | ~0.537 s | <0.001 s |
119+
| output orientation | 0.009 s | 0.016 s |
120+
| measured public core/orientation total | 2.762 s | 0.475 s |
121+
122+
After the cache change, traversal is still 98% of the measured core work;
123+
normalization, output orientation, cleanup, and NumPy handoff are individually
124+
too small to justify further single-threaded complexity in this pass.
125+
126+
Peak RSS for a single 160³ dense-mask Lewiner call fell from 323,048 KiB to
127+
235,760 KiB, a 27% reduction.
128+
129+
## Final results
130+
131+
`baseline / final` reports the speedup from this optimization pass.
132+
`skimage / final` above one means bioimage-cpp is faster.
133+
134+
| shape | workload | method | baseline ms | final ms | baseline / final | skimage / final |
135+
|---|---|---|---:|---:|---:|---:|
136+
| 48³ | sphere | lewiner | 2.56 | 1.85 | 1.38× | 1.34× |
137+
| 48³ | sphere | lorensen | 2.53 | 2.00 | 1.27× | 1.24× |
138+
| 48³ | dense mask | lewiner | 26.76 | 10.45 | 2.56× | 2.17× |
139+
| 48³ | dense mask | lorensen | 23.34 | 9.06 | 2.58× | 2.08× |
140+
| 48³ | scalar field | lewiner | 10.51 | 4.45 | 2.36× | 2.13× |
141+
| 48³ | scalar field | lorensen | 9.91 | 4.34 | 2.28× | 2.20× |
142+
| 96³ | sphere | lewiner | 15.15 | 12.15 | 1.25× | 1.35× |
143+
| 96³ | sphere | lorensen | 15.49 | 12.36 | 1.25× | 1.35× |
144+
| 96³ | dense mask | lewiner | 427.71 | 85.60 | 5.00× | 2.16× |
145+
| 96³ | dense mask | lorensen | 389.43 | 74.49 | 5.23× | 2.12× |
146+
| 96³ | scalar field | lewiner | 52.73 | 23.71 | 2.22× | 1.85× |
147+
| 96³ | scalar field | lorensen | 52.51 | 23.41 | 2.24× | 1.90× |
148+
| 160³ | sphere | lewiner | 63.65 | 56.15 | 1.13× | 1.14× |
149+
| 160³ | sphere | lorensen | 66.41 | 55.20 | 1.20× | 1.21× |
150+
| 160³ | dense mask | lewiner | 2,824.26 | 456.11 | 6.19× | 2.10× |
151+
| 160³ | dense mask | lorensen | 2,507.33 | 378.46 | 6.63× | 2.17× |
152+
| 160³ | scalar field | lewiner | 247.95 | 87.80 | 2.82× | 1.64× |
153+
| 160³ | scalar field | lorensen | 227.09 | 86.28 | 2.63× | 1.63× |
154+
155+
The size-dependent hash-map regression is gone. The implementation is faster
156+
than scikit-image in every measured case without threading, SIMD, or API
157+
changes.

0 commit comments

Comments
 (0)