From b1b5728cca854928bdfc553a91e07d2bb55fa794 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Fri, 10 Jul 2026 20:48:30 +0200 Subject: [PATCH 1/2] Implement marching cubes --- CMakeLists.txt | 1 + LICENSE | 30 + MIGRATION_GUIDE.md | 48 + README.md | 3 + THIRD_PARTY_NOTICES.md | 34 + development/mesh/PERFORMANCE_NOTES.md | 97 ++ development/mesh/_marching_cubes_reference.py | 129 +++ development/mesh/benchmark_marching_cubes.py | 235 +++++ development/mesh/check_marching_cubes.py | 125 +++ .../mesh/generate_marching_cubes_tables.py | 90 ++ .../bioimage_cpp/mesh/detail/mc33_luts.hxx | 393 ++++++++ include/bioimage_cpp/mesh/marching_cubes.hxx | 836 ++++++++++++++++++ src/bindings/mesh.cxx | 182 ++++ src/bindings/mesh.hxx | 9 + src/bindings/module.cxx | 2 + src/bioimage_cpp/__init__.py | 3 + src/bioimage_cpp/mesh/__init__.py | 169 ++++ tests/mesh/test_marching_cubes.py | 281 ++++++ 18 files changed, 2667 insertions(+) create mode 100644 THIRD_PARTY_NOTICES.md create mode 100644 development/mesh/PERFORMANCE_NOTES.md create mode 100644 development/mesh/_marching_cubes_reference.py create mode 100644 development/mesh/benchmark_marching_cubes.py create mode 100644 development/mesh/check_marching_cubes.py create mode 100644 development/mesh/generate_marching_cubes_tables.py create mode 100644 include/bioimage_cpp/mesh/detail/mc33_luts.hxx create mode 100644 include/bioimage_cpp/mesh/marching_cubes.hxx create mode 100644 src/bindings/mesh.cxx create mode 100644 src/bindings/mesh.hxx create mode 100644 src/bioimage_cpp/mesh/__init__.py create mode 100644 tests/mesh/test_marching_cubes.py diff --git a/CMakeLists.txt b/CMakeLists.txt index acb8ecb..7cc5323 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,6 +20,7 @@ nanobind_add_module(_core src/bindings/graph.cxx src/bindings/ground_truth.cxx src/bindings/label_multiset.cxx + src/bindings/mesh.cxx src/bindings/segmentation.cxx src/bindings/transformation.cxx src/bindings/util.cxx diff --git a/LICENSE b/LICENSE index 68ac1a2..b9e04cc 100644 --- a/LICENSE +++ b/LICENSE @@ -19,3 +19,33 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +Third-party notice — scikit-image Marching Cubes tables +-------------------------------------------------------- + +The source and binary distributions also include lookup-table data derived +from scikit-image 0.26.0. Copyright (c) 2009-2022 the scikit-image team. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index ace7291..7b32ad7 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -1872,6 +1872,54 @@ Notes: ## Skimage +### Marching Cubes + +`bioimage-cpp` provides dependency-free isosurface extraction under +`bic.mesh`, including both the topology-resolving Lewiner/MC33 method used by +default in scikit-image and the classic Lorensen lookup-table variant. + +```python +import bioimage_cpp as bic + +# Extract one object from a label image. `pad=True` closes objects that touch +# the volume boundary by adding a temporary zero-valued halo. +vertices, faces, normals, values = bic.mesh.marching_cubes( + labels == label_id, + level=0.5, + spacing=(z_spacing, y_spacing, x_spacing), + method="lewiner", + pad=True, +) +``` + +The signature follows `skimage.measure.marching_cubes`: `level`, `spacing`, +`gradient_direction`, `step_size`, `allow_degenerate`, `method`, and an +optional boolean `mask` have the same purpose. Coordinates use NumPy +`(z, y, x)` order. Vertices are `float32` at unit spacing and `float64` after +non-unit spacing; faces are consistently `int32`, and normals/values are +`float32`. + +Important details: + +- Inputs are converted to contiguous `float32` before extraction. Any real + numeric or boolean input dtype is accepted; complex inputs are rejected. +- `method="lewiner"` resolves ambiguous cases and is the default; + `method="lorensen"` selects the original 256-case algorithm. +- Normals and local-range values follow scikit-image semantics. As in + scikit-image, `gradient_direction` reverses face winding without changing + normals, and anisotropic spacing scales vertices without transforming + normals. +- `pad=False` matches scikit-image's open-boundary behavior. The additional + `pad=True` option uses a zero-valued halo and is intended for + foreground-positive segmentation masks. The iso-level is determined from + the original unpadded volume. +- Spacing entries must be positive and finite, and faces remain `int32` when + degenerate faces are removed. These validations/dtype choices are + intentional differences from scikit-image edge cases. + +See `development/mesh/check_marching_cubes.py` for reference comparisons and +`development/mesh/benchmark_marching_cubes.py` for reproducible timings. + ### Anti-Aliased Resampling `affine_transform` itself never pre-smooths the input; downsampling without diff --git a/README.md b/README.md index d61c4a9..73b6be8 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,9 @@ Image processing and segmentation functionality in C++ with light-weight python bindings through nanobind and minimal dependencies to enable distribution via pip. +The package includes dependency-free triangle-mesh extraction from 3D volumes +and segmentation masks under `bioimage_cpp.mesh`. + The `bioimage_cpp` python library can be installed via pip: ```bash pip install bioimage-cpp diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md new file mode 100644 index 0000000..ae10112 --- /dev/null +++ b/THIRD_PARTY_NOTICES.md @@ -0,0 +1,34 @@ +# Third-party notices + +## scikit-image Marching Cubes tables + +`include/bioimage_cpp/mesh/detail/mc33_luts.hxx` contains lookup-table data +derived from scikit-image 0.26.0's Marching Cubes implementation. The +corresponding MC33 control flow in `marching_cubes.hxx` follows that reference +port of the algorithm by Lewiner et al. The relevant scikit-image material is +licensed under the BSD 3-Clause License: + +Copyright (c) 2009-2022 the scikit-image team. All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. +3. Neither the name of the copyright holder nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/development/mesh/PERFORMANCE_NOTES.md b/development/mesh/PERFORMANCE_NOTES.md new file mode 100644 index 0000000..0537ff3 --- /dev/null +++ b/development/mesh/PERFORMANCE_NOTES.md @@ -0,0 +1,97 @@ +# Marching cubes performance + +`bic.mesh.marching_cubes` remains deterministic and single-threaded. The +optimization pass replaced the volume-growing edge hash map with Lewiner's +two-slice face cache. Several smaller candidates were measured first and +rejected when their gains did not survive repeated benchmarks. + +## Measurement setup + +- CPU: Intel Core i7-1185G7, 4 cores / 8 threads; the kernel uses one thread. +- OS: Linux 5.15 x86_64. +- Compiler/build: conda-forge GCC 14.3, editable release build (`-O3`). +- Python 3.13.13, NumPy 2.4.6, scikit-image 0.26.0. +- Workloads: a binary sphere, a reproducible 10%-foreground random binary + mask, and a deterministic smooth scalar field. +- Medium timings are medians of three seven-call batch medians. Small timings + use three batches of three calls; large timings use one batch of three calls. +- Every timed case passes the geometry/topology reference comparison first. + +Reproduce: + +```bash +python development/mesh/benchmark_marching_cubes.py --size small --repeats 3 --batches 3 +python development/mesh/benchmark_marching_cubes.py --size medium --repeats 7 --batches 3 +python development/mesh/benchmark_marching_cubes.py --size large --repeats 3 --batches 1 +``` + +## Candidate evaluation + +Each candidate was rebuilt, checked against the mesh tests and scikit-image +oracle, then timed twice at 96³. A candidate needed a repeatable improvement +of at least 3% without a regression above 3%. + +| candidate | observed result | decision | +|---|---|---| +| Slice-sized output reserves and scalar normal appends | First run helped spheres ~3%, repeat regressed spheres 3–4% and dense Lewiner 11% | Reverted | +| Hoisted flat indexing plus one-pass hash insertion | Dense masks improved ~2.5%, but spheres repeatedly regressed 3–4% | Reverted | +| Cached corner strengths and per-cell edge indices | Mostly 1–3% changes; dense Lorensen regressed ~2% | Reverted | +| Two-slice face cache | Repeatable 17–22% sphere, 55% scalar, and 80–81% dense-mask reductions at 96³ | Retained | +| Flat cube loads after the face cache | Helped sparse Lewiner, but repeat regressed scalar Lewiner 4% and left dense masks unchanged | Reverted | + +The retained cache uses two `int32` arrays with four slots per `(y, x)` +position. The upper edge layer becomes the next z-slice's lower layer; the new +upper layer is cleared. MC33 center vertices remain cell-local. Deduplication +memory is therefore `O(nx * ny)` instead of growing with the total surface. + +## Profiling + +Profiling used the repository's `BIOIMAGE_PROFILE` build on a 160³ Lewiner +call. The baseline's global map destruction happened after the original core +report, so it appears as the gap between core traversal/finalization and the +binding's `core_call`. + +| dense-mask phase | hash-map baseline | two-slice cache | +|---|---:|---:| +| cell traversal | 2.208 s | 0.449 s | +| normal finalization | 0.009 s | 0.009 s | +| cache/map cleanup | ~0.537 s | <0.001 s | +| output orientation | 0.009 s | 0.016 s | +| measured public core/orientation total | 2.762 s | 0.475 s | + +After the cache change, traversal is still 98% of the measured core work; +normalization, output orientation, cleanup, and NumPy handoff are individually +too small to justify further single-threaded complexity in this pass. + +Peak RSS for a single 160³ dense-mask Lewiner call fell from 323,048 KiB to +235,760 KiB, a 27% reduction. + +## Final results + +`baseline / final` reports the speedup from this optimization pass. +`skimage / final` above one means bioimage-cpp is faster. + +| shape | workload | method | baseline ms | final ms | baseline / final | skimage / final | +|---|---|---|---:|---:|---:|---:| +| 48³ | sphere | lewiner | 2.56 | 1.85 | 1.38× | 1.34× | +| 48³ | sphere | lorensen | 2.53 | 2.00 | 1.27× | 1.24× | +| 48³ | dense mask | lewiner | 26.76 | 10.45 | 2.56× | 2.17× | +| 48³ | dense mask | lorensen | 23.34 | 9.06 | 2.58× | 2.08× | +| 48³ | scalar field | lewiner | 10.51 | 4.45 | 2.36× | 2.13× | +| 48³ | scalar field | lorensen | 9.91 | 4.34 | 2.28× | 2.20× | +| 96³ | sphere | lewiner | 15.15 | 12.15 | 1.25× | 1.35× | +| 96³ | sphere | lorensen | 15.49 | 12.36 | 1.25× | 1.35× | +| 96³ | dense mask | lewiner | 427.71 | 85.60 | 5.00× | 2.16× | +| 96³ | dense mask | lorensen | 389.43 | 74.49 | 5.23× | 2.12× | +| 96³ | scalar field | lewiner | 52.73 | 23.71 | 2.22× | 1.85× | +| 96³ | scalar field | lorensen | 52.51 | 23.41 | 2.24× | 1.90× | +| 160³ | sphere | lewiner | 63.65 | 56.15 | 1.13× | 1.14× | +| 160³ | sphere | lorensen | 66.41 | 55.20 | 1.20× | 1.21× | +| 160³ | dense mask | lewiner | 2,824.26 | 456.11 | 6.19× | 2.10× | +| 160³ | dense mask | lorensen | 2,507.33 | 378.46 | 6.63× | 2.17× | +| 160³ | scalar field | lewiner | 247.95 | 87.80 | 2.82× | 1.64× | +| 160³ | scalar field | lorensen | 227.09 | 86.28 | 2.63× | 1.63× | + +The size-dependent hash-map regression is gone. The implementation is faster +than scikit-image in every measured case without threading, SIMD, or API +changes. diff --git a/development/mesh/_marching_cubes_reference.py b/development/mesh/_marching_cubes_reference.py new file mode 100644 index 0000000..8527c8b --- /dev/null +++ b/development/mesh/_marching_cubes_reference.py @@ -0,0 +1,129 @@ +"""scikit-image reference helpers for marching-cubes development checks. + +This module is intentionally kept outside the package and test dependency set. +It imports scikit-image lazily, mirrors the public ``pad`` extension, and is +used by the parity and benchmark scripts in this directory. +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np + + +def _skimage_marching_cubes(): + try: + from skimage.measure import marching_cubes + except ImportError as error: # pragma: no cover - development only + raise ImportError( + "scikit-image is required for the marching-cubes reference " + "(`pip install scikit-image`)" + ) from error + return marching_cubes + + +def reference_marching_cubes( + volume: np.ndarray, + level: float | None = None, + *, + spacing: Sequence[float] = (1.0, 1.0, 1.0), + gradient_direction: str = "descent", + step_size: int = 1, + allow_degenerate: bool = True, + method: str = "lewiner", + mask: np.ndarray | None = None, + pad: bool = False, +): + """Call scikit-image with the same public contract as ``bic.mesh``.""" + image = np.ascontiguousarray(volume, dtype=np.float32) + if level is None: + level = 0.5 * (float(image.min()) + float(image.max())) + mask_array = None if mask is None else np.ascontiguousarray(np.asarray(mask) != 0) + spacing_array = np.asarray(spacing, dtype=np.float64) + if pad: + image = np.pad(image, 1, mode="constant", constant_values=0) + if mask_array is not None: + mask_array = np.pad(mask_array, 1, mode="constant", constant_values=True) + result = _skimage_marching_cubes()( + image, + level, + spacing=spacing_array, + gradient_direction=gradient_direction, + step_size=step_size, + allow_degenerate=allow_degenerate, + method=method, + mask=mask_array, + ) + if pad: + vertices = result[0] - spacing_array.astype(result[0].dtype, copy=False) + result = (vertices, *result[1:]) + return result + + +def _sorted_rows(array: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + if len(array) == 0: + return array, np.empty(0, dtype=np.int64) + order = np.lexsort(tuple(array[:, axis] for axis in range(array.shape[1] - 1, -1, -1))) + return array[order], order + + +def _canonical_faces(vertices: np.ndarray, faces: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + unique_vertices, inverse = np.unique(vertices, axis=0, return_inverse=True) + canonical = np.sort(inverse[faces], axis=1) + canonical, _ = _sorted_rows(canonical) + return unique_vertices, canonical + + +def _surface_area(vertices: np.ndarray, faces: np.ndarray) -> float: + total = 0.0 + for begin in range(0, len(faces), 100_000): + triangles = vertices[faces[begin : begin + 100_000]].astype(np.float64, copy=False) + cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) + total += float(np.linalg.norm(cross, axis=1).sum()) * 0.5 + return total + + +def assert_mesh_matches(actual, reference, *, normal_atol: float = 1e-5) -> None: + """Assert geometry/topology parity independent of output ordering. + + Normals and local-range values are compared after sorting complete vertex + records by coordinates. Face winding is tested separately by the package + tests because canonical triangle comparison intentionally ignores it. + """ + actual_vertices, actual_faces, actual_normals, actual_values = actual + reference_vertices, reference_faces, reference_normals, reference_values = reference + + assert actual_vertices.shape == reference_vertices.shape + assert actual_faces.shape == reference_faces.shape + assert actual_normals.shape == reference_normals.shape + assert actual_values.shape == reference_values.shape + + actual_unique, actual_triangles = _canonical_faces(actual_vertices, actual_faces) + reference_unique, reference_triangles = _canonical_faces(reference_vertices, reference_faces) + np.testing.assert_allclose(actual_unique, reference_unique, rtol=0.0, atol=1e-6) + np.testing.assert_array_equal(actual_triangles, reference_triangles) + np.testing.assert_allclose( + _surface_area(actual_vertices, actual_faces), + _surface_area(reference_vertices, reference_faces), + rtol=1e-6, + atol=1e-8, + ) + + actual_records = np.column_stack( + (actual_vertices, actual_normals, actual_values) + ).astype(np.float64, copy=False) + reference_records = np.column_stack( + (reference_vertices, reference_normals, reference_values) + ).astype(np.float64, copy=False) + actual_records, _ = _sorted_rows(actual_records) + reference_records, _ = _sorted_rows(reference_records) + np.testing.assert_allclose( + actual_records[:, :3], reference_records[:, :3], rtol=0.0, atol=1e-6 + ) + np.testing.assert_allclose( + actual_records[:, 3:6], reference_records[:, 3:6], rtol=0.0, atol=normal_atol + ) + np.testing.assert_allclose( + actual_records[:, 6], reference_records[:, 6], rtol=0.0, atol=1e-6 + ) diff --git a/development/mesh/benchmark_marching_cubes.py b/development/mesh/benchmark_marching_cubes.py new file mode 100644 index 0000000..dbd5ecc --- /dev/null +++ b/development/mesh/benchmark_marching_cubes.py @@ -0,0 +1,235 @@ +"""Benchmark bioimage-cpp marching cubes against scikit-image. + +The benchmark runs a direct parity preflight before timing every workload. + +Examples +-------- +python development/mesh/benchmark_marching_cubes.py --size medium --repeats 7 +python development/mesh/benchmark_marching_cubes.py --size large --method lorensen +""" + +from __future__ import annotations + +import argparse +import json +import platform +from statistics import median +import sys +from time import perf_counter + +import numpy as np +import skimage + +import bioimage_cpp as bic + +from _marching_cubes_reference import assert_mesh_matches, reference_marching_cubes + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--size", choices=("small", "medium", "large"), default="medium") + parser.add_argument("--method", choices=("lewiner", "lorensen", "all"), default="all") + parser.add_argument( + "--workload", + choices=("binary_sphere", "dense_binary_mask", "scalar_field", "all"), + default="all", + ) + parser.add_argument( + "--backend", choices=("both", "bic", "skimage"), default="both" + ) + parser.add_argument("--repeats", type=int, default=7) + parser.add_argument("--warmup", type=int, default=1) + parser.add_argument("--batches", type=int, default=3) + parser.add_argument("--seed", type=int, default=20260709) + parser.add_argument("--json", default="", help="optional JSON result path") + parser.add_argument("--baseline", default="", help="optional prior benchmark JSON for bic relative changes") + return parser.parse_args() + + +def shape_for(size: str) -> tuple[int, int, int]: + return {"small": (48, 48, 48), "medium": (96, 96, 96), "large": (160, 160, 160)}[size] + + +def workloads(shape: tuple[int, int, int], seed: int): + z, y, x = np.ogrid[: shape[0], : shape[1], : shape[2]] + center = np.asarray(shape, dtype=np.float32) / 2.0 + radius = min(shape) * 0.28 + sphere = ((z - center[0]) ** 2 + (y - center[1]) ** 2 + (x - center[2]) ** 2 < radius**2).astype(np.uint8) + rng = np.random.default_rng(seed) + dense_mask = (rng.random(shape) < 0.10).astype(np.uint8) + zf = z.astype(np.float32) / max(shape[0] - 1, 1) + yf = y.astype(np.float32) / max(shape[1] - 1, 1) + xf = x.astype(np.float32) / max(shape[2] - 1, 1) + scalar = ( + np.sin(4.0 * np.pi * zf) + + np.cos(6.0 * np.pi * yf) + + np.sin(5.0 * np.pi * xf) + ).astype(np.float32) + return [ + ("binary_sphere", sphere, 0.5), + ("dense_binary_mask", dense_mask, 0.5), + ("scalar_field", scalar, 0.0), + ] + + +def time_call(function, repeats: int, warmup: int, batches: int) -> dict[str, object]: + raw_timings = [] + batch_medians = [] + for _ in range(batches): + for _ in range(warmup): + function() + timings = [] + for _ in range(repeats): + start = perf_counter() + function() + timings.append(perf_counter() - start) + raw_timings.extend(timings) + batch_medians.append(median(timings)) + return { + "raw_s": raw_timings, + "batch_medians_s": batch_medians, + "median_s": median(batch_medians), + "p10_s": float(np.percentile(raw_timings, 10)), + "p90_s": float(np.percentile(raw_timings, 90)), + "min_s": min(raw_timings), + } + + +def load_baseline(path: str) -> dict[tuple[str, str, tuple[int, int, int]], float]: + if not path: + return {} + with open(path) as file: + rows = json.load(file)["results"] + return { + (row["workload"], row["method"], tuple(row["shape"])): row["bic_median_s"] + for row in rows + } + + +def main() -> int: + args = parse_args() + if args.repeats < 1 or args.warmup < 0 or args.batches < 1: + raise SystemExit("repeats and batches must be >= 1 and warmup must be >= 0") + methods = ("lewiner", "lorensen") if args.method == "all" else (args.method,) + shape = shape_for(args.size) + n_voxels = int(np.prod(shape)) + baseline = load_baseline(args.baseline) + rows = [] + + print(f"shape={shape} voxels={n_voxels} repeats={args.repeats} batches={args.batches}") + print(f"{'workload/method':<28} {'V':>8} {'F':>8} {'bic ms':>11} {'p10-p90 ms':>15} {'skimage ms':>12} {'speed':>9} {'Mvox/s':>10} {'delta':>9}") + print("-" * 125) + selected_workloads = workloads(shape, args.seed) + if args.workload != "all": + selected_workloads = [row for row in selected_workloads if row[0] == args.workload] + for workload, volume, level in selected_workloads: + for method in methods: + kwargs = {"method": method} + start = perf_counter() + actual = bic.mesh.marching_cubes(volume, level, **kwargs) + bic_first_call = perf_counter() - start + start = perf_counter() + reference = reference_marching_cubes(volume, level, **kwargs) + reference_first_call = perf_counter() - start + assert_mesh_matches(actual, reference) + bic_times = None + if args.backend in ("both", "bic"): + bic_times = time_call( + lambda: bic.mesh.marching_cubes(volume, level, **kwargs), + args.repeats, + args.warmup, + args.batches, + ) + reference_times = None + if args.backend in ("both", "skimage"): + reference_times = time_call( + lambda: reference_marching_cubes(volume, level, **kwargs), + args.repeats, + args.warmup, + args.batches, + ) + bic_median = None if bic_times is None else bic_times["median_s"] + reference_median = ( + None if reference_times is None else reference_times["median_s"] + ) + baseline_time = baseline.get((workload, method, shape)) + relative_change = ( + None + if baseline_time is None or bic_median is None + else bic_median / baseline_time - 1.0 + ) + speedup = ( + None + if bic_median is None or reference_median is None + else reference_median / bic_median + ) + row = { + "workload": workload, + "method": method, + "shape": shape, + "vertices": len(actual[0]), + "faces": len(actual[1]), + "bic_first_call_s": bic_first_call, + "reference_first_call_s": reference_first_call, + "bic_raw_s": None if bic_times is None else bic_times["raw_s"], + "bic_batch_medians_s": None if bic_times is None else bic_times["batch_medians_s"], + "bic_median_s": bic_median, + "bic_min_s": None if bic_times is None else bic_times["min_s"], + "bic_p10_s": None if bic_times is None else bic_times["p10_s"], + "bic_p90_s": None if bic_times is None else bic_times["p90_s"], + "reference_raw_s": None if reference_times is None else reference_times["raw_s"], + "reference_batch_medians_s": None if reference_times is None else reference_times["batch_medians_s"], + "reference_median_s": reference_median, + "reference_min_s": None if reference_times is None else reference_times["min_s"], + "speedup": speedup, + "mvox_per_s": None if bic_median is None else n_voxels / bic_median / 1e6, + "baseline_bic_median_s": baseline_time, + "relative_change": relative_change, + } + rows.append(row) + delta = "-" if row["relative_change"] is None else f"{100.0 * row['relative_change']:+.1f}%" + bic_text = "-" if bic_median is None else f"{bic_median * 1e3:>11.2f}" + spread_text = ( + " - " + if bic_times is None + else f"{bic_times['p10_s'] * 1e3:>6.2f}-{bic_times['p90_s'] * 1e3:>6.2f}" + ) + reference_text = ( + "-" if reference_median is None else f"{reference_median * 1e3:>12.2f}" + ) + speed_text = "-" if speedup is None else f"{speedup:>8.2f}x" + throughput_text = ( + "-" if row["mvox_per_s"] is None else f"{row['mvox_per_s']:>10.2f}" + ) + print( + f"{workload + '/' + method:<28} {row['vertices']:>8} {row['faces']:>8} " + f"{bic_text:>11} {spread_text:>15} {reference_text:>12} " + f"{speed_text:>9} {throughput_text:>10} {delta:>9}" + ) + if args.json: + with open(args.json, "w") as file: + json.dump( + { + "shape": shape, + "repeats": args.repeats, + "warmup": args.warmup, + "batches": args.batches, + "seed": args.seed, + "workload": args.workload, + "backend": args.backend, + "environment": { + "python": sys.version, + "platform": platform.platform(), + "numpy": np.__version__, + "scikit_image": skimage.__version__, + }, + "results": rows, + }, + file, + indent=2, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/development/mesh/check_marching_cubes.py b/development/mesh/check_marching_cubes.py new file mode 100644 index 0000000..b8e3a24 --- /dev/null +++ b/development/mesh/check_marching_cubes.py @@ -0,0 +1,125 @@ +"""Compare ``bic.mesh.marching_cubes`` to ``skimage.measure.marching_cubes``. + +Run from the repository root: + + python development/mesh/check_marching_cubes.py +""" + +from __future__ import annotations + +import argparse +import sys + +import numpy as np + +import bioimage_cpp as bic + +from _marching_cubes_reference import assert_mesh_matches, reference_marching_cubes + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--random-cases", + type=int, + default=1024, + help="number of deterministic random scalar one-cube configurations to compare (default: 1024)", + ) + return parser.parse_args() + + +def cases(): + box = np.zeros((7, 8, 9), dtype=np.uint8) + box[1:6, 2:7, 1:8] = 1 + + boundary = np.zeros((6, 7, 8), dtype=np.uint8) + boundary[:3, 2:6, 2:7] = 1 + + rng = np.random.default_rng(20260709) + scalar = rng.normal(size=(8, 9, 10)).astype(np.float32) + binary = rng.integers(0, 2, size=(8, 9, 10), dtype=np.uint8) + roi = np.ones_like(binary, dtype=bool) + roi[:, :2] = False + degenerate = rng.integers(0, 3, size=(8, 9, 10), dtype=np.uint8) + + positive_background = np.full((6, 7, 8), 2.0, dtype=np.float32) + positive_background[1:5, 2:6, 2:7] = 4.0 + + ambiguous = np.array([(6 >> bit) & 1 for bit in range(8)], dtype=np.uint8).reshape(2, 2, 2) + return [ + ("box/lewiner", box, 0.5, {"method": "lewiner"}), + ("box/lorensen", box, 0.5, {"method": "lorensen"}), + ("scalar/lewiner", scalar, 0.0, {"method": "lewiner"}), + ("binary/lorensen-step", binary, 0.5, {"method": "lorensen", "step_size": 2}), + ("masked/lewiner", binary, 0.5, {"mask": roi}), + ("masked/lewiner-step3", binary, 0.5, {"mask": roi, "step_size": 3}), + ("degenerate/removed", degenerate, 1.0, {"allow_degenerate": False}), + ("anisotropic/ascent", box, 0.5, {"spacing": (2.0, 0.5, 3.0), "gradient_direction": "ascent"}), + ("ambiguous/lewiner", ambiguous, 0.5, {"method": "lewiner"}), + ("ambiguous/lorensen", ambiguous, 0.5, {"method": "lorensen"}), + ("boundary/padded", boundary, 0.5, {"pad": True}), + ("padded/default-level", positive_background, None, {"pad": True}), + ] + + +def assert_case(volume: np.ndarray, level: float | None, kwargs: dict[str, object]) -> None: + actual = bic.mesh.marching_cubes(volume, level, **kwargs) + reference = reference_marching_cubes(volume, level, **kwargs) + assert_mesh_matches(actual, reference) + + +def check_all_binary_cube_configurations() -> None: + for configuration in range(1, 255): + volume = np.array( + [(configuration >> bit) & 1 for bit in range(8)], dtype=np.uint8 + ).reshape(2, 2, 2) + for method in ("lewiner", "lorensen"): + assert_case(volume, 0.5, {"method": method}) + + +def check_random_scalar_cubes(n_cases: int) -> None: + rng = np.random.default_rng(20260710) + checked = 0 + while checked < n_cases: + volume = rng.normal(size=(2, 2, 2)).astype(np.float32) + if np.all(volume <= 0.0) or np.all(volume > 0.0): + continue + for method in ("lewiner", "lorensen"): + assert_case(volume, 0.0, {"method": method}) + checked += 1 + + +def main() -> int: + args = parse_args() + if args.random_cases < 0: + raise SystemExit("random-cases must be >= 0") + failed = False + print(f"{'case':<26} {'vertices':>10} {'faces':>10} {'status':>8}") + print("-" * 60) + for name, volume, level, kwargs in cases(): + try: + actual = bic.mesh.marching_cubes(volume, level, **kwargs) + reference = reference_marching_cubes(volume, level, **kwargs) + assert_mesh_matches(actual, reference) + except Exception as error: + print(f"{name:<26} {'-':>10} {'-':>10} {'FAIL':>8} {error}") + failed = True + continue + print(f"{name:<26} {len(actual[0]):>10} {len(actual[1]):>10} {'OK':>8}") + if not failed: + try: + check_all_binary_cube_configurations() + print(f"{'all binary one-cube cases':<26} {'254':>10} {'508':>10} {'OK':>8}") + check_random_scalar_cubes(args.random_cases) + print(f"{'random scalar one-cube cases':<26} {args.random_cases:>10} {2 * args.random_cases:>10} {'OK':>8}") + except Exception as error: + print(f"{'randomized parity':<26} {'-':>10} {'-':>10} {'FAIL':>8} {error}") + failed = True + if failed: + print("marching-cubes parity check failed", file=sys.stderr) + return 1 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/development/mesh/generate_marching_cubes_tables.py b/development/mesh/generate_marching_cubes_tables.py new file mode 100644 index 0000000..fffdfa5 --- /dev/null +++ b/development/mesh/generate_marching_cubes_tables.py @@ -0,0 +1,90 @@ +"""Regenerate the committed Marching Cubes lookup-table header. + +The table payloads come from scikit-image's BSD-3-Clause licensed +``_marching_cubes_lewiner_luts`` module. The generated header is committed, so +scikit-image remains a development-only dependency. + +Run from the repository root:: + + python development/mesh/generate_marching_cubes_tables.py +""" + +from __future__ import annotations + +from pathlib import Path +from textwrap import wrap + +import skimage +import skimage.measure._marching_cubes_lewiner_luts as luts + + +TABLE_NAMES = ( + "CASESCLASSIC", "CASES", "TILING1", "TILING2", "TILING3_1", + "TILING3_2", "TILING4_1", "TILING4_2", "TILING5", "TILING6_1_1", + "TILING6_1_2", "TILING6_2", "TILING7_1", "TILING7_2", "TILING7_3", + "TILING7_4_1", "TILING7_4_2", "TILING8", "TILING9", + "TILING10_1_1", "TILING10_1_1_", "TILING10_1_2", "TILING10_2", + "TILING10_2_", "TILING11", "TILING12_1_1", "TILING12_1_1_", + "TILING12_1_2", "TILING12_2", "TILING12_2_", "TILING13_1", + "TILING13_1_", "TILING13_2", "TILING13_2_", "TILING13_3", + "TILING13_3_", "TILING13_4", "TILING13_5_1", "TILING13_5_2", + "TILING14", "TEST3", "TEST4", "TEST6", "TEST7", "TEST10", + "TEST12", "TEST13", "SUBCONFIG13", +) + + +def generate_header() -> str: + lines = [ + "#pragma once", + "", + "// AUTO-GENERATED by development/mesh/generate_marching_cubes_tables.py.", + "// Lookup tables derived from scikit-image 0.26.0, BSD-3-Clause.", + "// The original implementation is based on Lewiner et al. (2003).", + "// See THIRD_PARTY_NOTICES.md for the complete license notice.", + "", + "#include ", + "#include ", + "", + "namespace bioimage_cpp::mesh::detail {", + "", + "struct EncodedLut {", + " std::array shape;", + " int ndim;", + " std::string_view data;", + "};", + "", + ] + for name in TABLE_NAMES: + shape, payload = getattr(luts, name) + padded_shape = tuple(shape) + (0,) * (3 - len(shape)) + shape_text = ", ".join(str(value) for value in padded_shape) + compact_payload = "".join(payload.split()) + lines.append( + f"inline constexpr EncodedLut k{name}{{{{{shape_text}}}, {len(shape)}," + ) + chunks = wrap(compact_payload, width=112) + for chunk in chunks: + lines.append(f' "{chunk}"') + lines.append("};") + lines.append("") + lines.append("} // namespace bioimage_cpp::mesh::detail") + lines.append("") + return "\n".join(lines) + + +def main() -> None: + if skimage.__version__ != "0.26.0": + raise RuntimeError( + "table regeneration is pinned to scikit-image 0.26.0, got " + f"{skimage.__version__}" + ) + output = ( + Path(__file__).resolve().parents[2] + / "include/bioimage_cpp/mesh/detail/mc33_luts.hxx" + ) + output.write_text(generate_header()) + print(f"wrote {output}") + + +if __name__ == "__main__": + main() diff --git a/include/bioimage_cpp/mesh/detail/mc33_luts.hxx b/include/bioimage_cpp/mesh/detail/mc33_luts.hxx new file mode 100644 index 0000000..5eb1f33 --- /dev/null +++ b/include/bioimage_cpp/mesh/detail/mc33_luts.hxx @@ -0,0 +1,393 @@ +#pragma once + +// AUTO-GENERATED by development/mesh/generate_marching_cubes_tables.py. +// Lookup tables derived from scikit-image 0.26.0, BSD-3-Clause. +// The original implementation is based on Lewiner et al. (2003). +// See THIRD_PARTY_NOTICES.md for the complete license notice. + +#include +#include + +namespace bioimage_cpp::mesh::detail { + +struct EncodedLut { + std::array shape; + int ndim; + std::string_view data; +}; + +inline constexpr EncodedLut kCASESCLASSIC{{256, 16, 0}, 2, + "/////////////////////wAIA/////////////////8AAQn/////////////////AQgDCQgB/////////////wECCv////////////////8ACAMB" + "Agr/////////////CQIKAAIJ/////////////wIIAwIKCAoJCP////////8DCwL/////////////////AAsCCAsA/////////////wEJAAIDC///" + "//////////8BCwIBCQsJCAv/////////AwoBCwoD/////////////wAKAQAICggLCv////////8DCQADCwkLCgn/////////CQgKCggL////////" + "/////wQHCP////////////////8EAwAHAwT/////////////AAEJCAQH/////////////wQBCQQHAQcDAf////////8BAgoIBAf/////////////" + "AwQHAwAEAQIK/////////wkCCgkAAggEB/////////8CCgkCCQcCBwMHCQT/////CAQHAwsC/////////////wsEBwsCBAIABP////////8JAAEI" + "BAcCAwv/////////BAcLCQQLCQsCCQIB/////wMKAQMLCgcIBP////////8BCwoBBAsBAAQHCwT/////BAcICQALCQsKCwAD/////wQHCwQLCQkL" + "Cv////////8JBQT/////////////////CQUEAAgD/////////////wAFBAEFAP////////////8IBQQIAwUDAQX/////////AQIKCQUE////////" + "/////wMACAECCgQJBf////////8FAgoFBAIEAAL/////////AgoFAwIFAwUEAwQI/////wkFBAIDC/////////////8ACwIACAsECQX/////////" + "AAUEAAEFAgML/////////wIBBQIFCAIICwQIBf////8KAwsKAQMJBQT/////////BAkFAAgBCAoBCAsK/////wUEAAUACwULCgsAA/////8FBAgF" + "CAoKCAv/////////CQcIBQcJ/////////////wkDAAkFAwUHA/////////8ABwgAAQcBBQf/////////AQUDAwUH/////////////wkHCAkFBwoB" + "Av////////8KAQIJBQAFAwAFBwP/////CAACCAIFCAUHCgUC/////wIKBQIFAwMFB/////////8HCQUHCAkDCwL/////////CQUHCQcCCQIAAgcL" + "/////wIDCwABCAEHCAEFB/////8LAgELAQcHAQX/////////CQUICAUHCgEDCgML/////wUHAAUACQcLAAEACgsKAP8LCgALAAMKBQAIAAcFBwD/" + "CwoFBwsF/////////////woGBf////////////////8ACAMFCgb/////////////CQABBQoG/////////////wEIAwEJCAUKBv////////8BBgUC" + "BgH/////////////AQYFAQIGAwAI/////////wkGBQkABgACBv////////8FCQgFCAIFAgYDAgj/////AgMLCgYF/////////////wsACAsCAAoG" + "Bf////////8AAQkCAwsFCgb/////////BQoGAQkCCQsCCQgL/////wYDCwYFAwUBA/////////8ACAsACwUABQEFCwb/////AwsGAAMGAAYFAAUJ" + "/////wYFCQYJCwsJCP////////8FCgYEBwj/////////////BAMABAcDBgUK/////////wEJAAUKBggEB/////////8KBgUBCQcBBwMHCQT/////" + "BgECBgUBBAcI/////////wECBQUCBgMABAMEB/////8IBAcJAAUABgUAAgb/////BwMJBwkEAwIJBQkGAgYJ/wMLAgcIBAoGBf////////8FCgYE" + "BwIEAgACBwv/////AAEJBAcIAgMLBQoG/////wkCAQkLAgkECwcLBAUKBv8IBAcDCwUDBQEFCwb/////BQELBQsGAQALBwsEAAQL/wAFCQAGBQAD" + "BgsGAwgEB/8GBQkGCQsEBwkHCwn/////CgQJBgQK/////////////wQKBgQJCgAIA/////////8KAAEKBgAGBAD/////////CAMBCAEGCAYEBgEK" + "/////wEECQECBAIGBP////////8DAAgBAgkCBAkCBgT/////AAIEBAIG/////////////wgDAggCBAQCBv////////8KBAkKBgQLAgP/////////" + "AAgCAggLBAkKBAoG/////wMLAgABBgAGBAYBCv////8GBAEGAQoECAECAQsICwH/CQYECQMGCQEDCwYD/////wgLAQgBAAsGAQkBBAYEAf8DCwYD" + "BgAABgT/////////BgQICwYI/////////////wcKBgcICggJCv////////8ABwMACgcACQoGBwr/////CgYHAQoHAQcIAQgA/////woGBwoHAQEH" + "A/////////8BAgYBBggBCAkIBgf/////AgYJAgkBBgcJAAkDBwMJ/wcIAAcABgYAAv////////8HAwIGBwL/////////////AgMLCgYICggJCAYH" + "/////wIABwIHCwAJBwYHCgkKB/8BCAABBwgBCgcGBwoCAwv/CwIBCwEHCgYBBgcB/////wgJBggGBwkBBgsGAwEDBv8ACQELBgf/////////////" + "BwgABwAGAwsACwYA/////wcLBv////////////////8HBgv/////////////////AwAICwcG/////////////wABCQsHBv////////////8IAQkI" + "AwELBwb/////////CgECBgsH/////////////wECCgMACAYLB/////////8CCQACCgkGCwf/////////BgsHAgoDCggDCgkI/////wcCAwYCB///" + "//////////8HAAgHBgAGAgD/////////AgcGAgMHAAEJ/////////wEGAgEIBgEJCAgHBv////8KBwYKAQcBAwf/////////CgcGAQcKAQgHAQAI" + "/////wADBwAHCgAKCQYKB/////8HBgoHCggICgn/////////BggECwgG/////////////wMGCwMABgAEBv////////8IBgsIBAYJAAH/////////" + "CQQGCQYDCQMBCwMG/////wYIBAYLCAIKAf////////8BAgoDAAsABgsABAb/////BAsIBAYLAAIJAgoJ/////woJAwoDAgkEAwsDBgQGA/8IAgMI" + "BAIEBgL/////////AAQCBAYC/////////////wEJAAIDBAIEBgQDCP////8BCQQBBAICBAb/////////CAEDCAYBCAQGBgoB/////woBAAoABgYA" + "BP////////8EBgMEAwgGCgMAAwkKCQP/CgkEBgoE/////////////wQJBQcGC/////////////8ACAMECQULBwb/////////BQABBQQABwYL////" + "/////wsHBggDBAMFBAMBBf////8JBQQKAQIHBgv/////////BgsHAQIKAAgDBAkF/////wcGCwUECgQCCgQAAv////8DBAgDBQQDAgUKBQILBwb/" + "BwIDBwYCBQQJ/////////wkFBAAIBgAGAgYIB/////8DBgIDBwYBBQAFBAD/////BgIIBggHAgEIBAgFAQUI/wkFBAoBBgEHBgEDB/////8BBgoB" + "BwYBAAcIBwAJBQT/BAAKBAoFAAMKBgoHAwcK/wcGCgcKCAUECgQICv////8GCQUGCwkLCAn/////////AwYLAAYDAAUGAAkF/////wALCAAFCwAB" + "BQUGC/////8GCwMGAwUFAwH/////////AQIKCQULCQsICwUG/////wALAwAGCwAJBgUGCQECCv8LCAULBQYIAAUKBQIAAgX/BgsDBgMFAgoDCgUD" + "/////wUICQUCCAUGAgMIAv////8JBQYJBgAABgL/////////AQUIAQgABQYIAwgCBgII/wEFBgIBBv////////////8BAwYBBgoDCAYFBgkICQb/" + "CgEACgAGCQUABQYA/////wADCAUGCv////////////8KBQb/////////////////CwUKBwUL/////////////wsFCgsHBQgDAP////////8FCwcF" + "CgsBCQD/////////CgcFCgsHCQgBCAMB/////wsBAgsHAQcFAf////////8ACAMBAgcBBwUHAgv/////CQcFCQIHCQACAgsH/////wcFAgcCCwUJ" + "AgMCCAkIAv8CBQoCAwUDBwX/////////CAIACAUCCAcFCgIF/////wkAAQUKAwUDBwMKAv////8JCAIJAgEIBwIKAgUHBQL/AQMFAwcF////////" + "/////wAIBwAHAQEHBf////////8JAAMJAwUFAwf/////////CQgHBQkH/////////////wUIBAUKCAoLCP////////8FAAQFCwAFCgsLAwD/////" + "AAEJCAQKCAoLCgQF/////woLBAoEBQsDBAkEAQMBBP8CBQECCAUCCwgEBQj/////AAQLAAsDBAULAgsBBQEL/wACBQAFCQILBQQFCAsIBf8JBAUC" + "CwP/////////////AgUKAwUCAwQFAwgE/////wUKAgUCBAQCAP////////8DCgIDBQoDCAUEBQgAAQn/BQoCBQIEAQkCCQQC/////wgEBQgFAwMF" + "Af////////8ABAUBAAX/////////////CAQFCAUDCQAFAAMF/////wkEBf////////////////8ECwcECQsJCgv/////////AAgDBAkHCQsHCQoL" + "/////wEKCwELBAEEAAcEC/////8DAQQDBAgBCgQHBAsKCwT/BAsHCQsECQILCQEC/////wkHBAkLBwkBCwILAQAIA/8LBwQLBAICBAD/////////" + "CwcECwQCCAMEAwIE/////wIJCgIHCQIDBwcECf////8JCgcJBwQKAgcIBwACAAf/AwcKAwoCBwQKAQoABAAK/wEKAggHBP////////////8ECQEE" + "AQcHAQP/////////BAkBBAEHAAgBCAcB/////wQAAwcEA/////////////8ECAf/////////////////CQoICgsI/////////////wMACQMJCwsJ" + "Cv////////8AAQoACggICgv/////////AwEKCwMK/////////////wECCwELCQkLCP////////8DAAkDCQsBAgkCCwn/////AAILCAAL////////" + "/////wMCC/////////////////8CAwgCCAoKCAn/////////CQoCAAkC/////////////wIDCAIICgABCAEKCP////8BCgL/////////////////" + "AQMICQEI/////////////wAJAf////////////////8AAwj//////////////////////////////////////w==" +}; + +inline constexpr EncodedLut kCASES{{256, 2, 0}, 2, + "AP8BAAEBAgABAgMAAgMFAAEDAgEDAwUBAgUFBAUJCAABBAICAwQFAgQCBgIGCQsAAwgFBQcDCQEGEA4DDAwFGAEFAwECBAUDAwYHAAUKCQAEAwYE" + "BgsOAQYRDAQLBgUZAggFBwUMCAEGEgwFDgcFHAYVCwQMDwUeCgUGIAYnAgwBBgQAAwUGAAIGBgMFCw4AAwkGBQcEDAEFDgsDCQQFGgMKBgYHBQwC" + "BhMKAQwNBhgHBwwJDQEHCQwUBiEHDQMMAgoGBwUNCwIFEAwHCAMFHQYWCgIMEQYbDgkGIgUnAg4FFA4FCQUFIAsKBiMFKQIQDBcGJQcOAxAGLgQG" + "AxUBCAEHAwIEAQYBAwcHAQYKDAACBwUGBgwLAQUPCQIOBgUbAgkFCAYNDgIGFAwGCgMGGQUSCAIMEAUfCwkFIgYoAg0DCwcCBg4MAwcGDQAMDgcI" + "BhcMCgoEBhwMFQcKBikDDQUVCQMLCAUhDBYHCwYqAw4OCwUkBiwCEQYvAxIEBwEJAgsGCAYPCgAFEQwICwcGGgUTDgQMEgYdCAQFIwUoAg8FFgsF" + "DBMGHg4KBiQGKwQECQcFJQcPAxEFLAITAxYBCgUXDAsOCAYfCQYHDAUqAw8LCwYmBi0EBQUtAxMCFQELCAUFJgUrAhIFLgMUAhYBDAUvAhQDFwEN" + "AhcBDgEPAP8=" +}; + +inline constexpr EncodedLut kTILING1{{16, 3, 0}, 2, + "AAgDAAEJAQIKAwsCBAcICQUECgYFBwYLBwsGCgUGCQQFBAgHAwILAQoCAAkBAAMI" +}; + +inline constexpr EncodedLut kTILING2{{24, 6, 0}, 2, + "AQgDCQgBAAsCCAsABAMABwMECQIKAAIJAAUEAQUAAwoBCwoDAQYFAgYBBwIDBgIHCQcIBQcJBggECwgGCgQJBgQKCwUKBwULCwoFBwsFCgkEBgoE" + "BgQICwYICQgHBQkHBwMCBgcCAQUGAgEGAwEKCwMKAAQFAQAFCQoCAAkCBAADBwQDAAILCAALAQMICQEI" +}; + +inline constexpr EncodedLut kTILING3_1{{24, 6, 0}, 2, + "AAgDAQIKCQUEAAgDAwAICwcGAQkAAgMLAAEJCAQHCQABBQoGAQIKCQUECgECBgsHCAQHAwsCAgMLCgYFBQoGBAcIBAkFBwYLBQkECwYHBgoFCAcE" + "CwMCBQYKBwQIAgsDAgEKBwsGCgIBBAUJAQAJBgoFCQEABwQIAAkBCwMCCAADBgcLBAUJAwgAAwgACgIB" +}; + +inline constexpr EncodedLut kTILING3_2{{24, 12, 0}, 2, + "CgMCCggDCgEACAoAAwQIAwUEAwAJBQMJBggHBgAIBgsDAAYDCwADCwkACwIBCQsBBwkEBwEJBwgAAQcABgEKBgABCQAGCQYFBAoFBAIKBAkBAgQB" + "BwILBwECBwYKAQcKAgcLAgQHAgMIBAIIBQsGBQMLBQoCAwUCCAYHCAoGCAQFCggFCwUGCwkFCwcECQsEBgULBQkLBAcLBAsJBwYIBgoIBQQIBQgK" + "BgsFCwMFAgoFAgUDCwcCBwQCCAMCCAIECwIHAgEHCgYHCgcBBQoECgIEAQkEAQQCCgEGAQAGBgAJBQYJBAkHCQEHAAgHAAcBAwALAAkLAQILAQsJ" + "BwgGCAAGAwsGAwYACAQDBAUDCQADCQMFAgMKAwgKAAEKAAoI" +}; + +inline constexpr EncodedLut kTILING4_1{{8, 6, 0}, 2, + "AAgDBQoGAAEJCwcGAQIKCAQHCQUEAgMLBAUJCwMCCgIBBwQICQEABgcLAwgABgoF" +}; + +inline constexpr EncodedLut kTILING4_2{{8, 18, 0}, 2, + "CAUABQgGAwYIBgMKAAoDCgAFCQYBBgkHAAcJBwALAQsACwEGCgcCBwoEAQQKBAEIAggBCAIHCwQDBAsFAgULBQIJAwkCCQMEAwQLBQsECwUCCQIF" + "AgkDBAMJAgcKBAoHCgQBCAEEAQgCBwIIAQYJBwkGCQcACwAHAAsBBgELAAUIBggFCAYDCgMGAwoABQAK" +}; + +inline constexpr EncodedLut kTILING5{{48, 9, 0}, 2, + "AggDAgoICgkIAQsCAQkLCQgLBAEJBAcBBwMBCAUECAMFAwEFAAoBAAgKCAsKCwQHCwIEAgAEBwAIBwYABgIACQMACQUDBQcDAwYLAwAGAAQGAwkA" + "AwsJCwoJBQIKBQQCBAACCQYFCQAGAAIGAAcIAAEHAQUHCgABCgYABgQABgMLBgUDBQEDCgcGCgEHAQMHAQQJAQIEAgYECwECCwcBBwUBCAIDCAQC" + "BAYCAgUKAgMFAwcFBwoGBwgKCAkKBgkFBgsJCwgJBQgEBQoICgsIBAsHBAkLCQoLBAcLBAsJCQsKBQQIBQgKCggLBgUJBgkLCwkIBwYKBwoICAoJ" + "AgoFAgUDAwUHCAMCCAIEBAIGCwIBCwEHBwEFAQkEAQQCAgQGCgYHCgcBAQcDBgsDBgMFBQMBCgEACgAGBgAEAAgHAAcBAQcFCQUGCQYAAAYCBQoC" + "BQIEBAIAAwAJAwkLCwkKAwsGAwYAAAYECQADCQMFBQMHBwgABwAGBgACCwcECwQCAgQAAAEKAAoICAoLCAQFCAUDAwUBBAkBBAEHBwEDAQILAQsJ" + "CQsIAgMIAggKCggJ" +}; + +inline constexpr EncodedLut kTILING6_1_1{{48, 9, 0}, 2, + "BgUKAwEICQgBCwcGCQMBAwkIAQIKBwAEAAcDAwAIBQIGAgUBBQQJAgALCAsACgYFCAIAAggLCgYFAAQDBwMEAwAIBgQKCQoECAMACgcFBwoLCAQH" + "CgACAAoJBwYLAAIJCgkCAgMLBAEFAQQAAAEJBgMHAwYCCQABCwQGBAsICwcGAQUABAAFAAEJBwULCgsFBAcIAQMKCwoDCQUECwEDAQsKCgECCAUH" + "BQgJCAQHAgYBBQEGAQIKBAYICwgGAgMLBQcJCAkHCwIDCQYEBgkKCQUEAwcCBgIHBAUJAgcDBwIGAwILBAYJCgkGCwMCCQcFBwkICgIBCAYEBggL" + "BwQIAQYCBgEFAgEKBwUICQgFBAUJAwELCgsBCAcECgMBAwoLCQEACwUHBQsKBgcLAAUBBQAEAQAJBgQLCAsECQEABwMGAgYDCwMCBQEEAAQBCwYH" + "CQIAAgkKBwQIAgAKCQoAAAMIBQcKCwoHCAADCgQGBAoJBQYKAwQABAMHBQYKAAIICwgCCQQFCwACAAsICAADBgIFAQUCCgIBBAAHAwcABgcLAQMJ" + "CAkDCgUGCAEDAQgJ" +}; + +inline constexpr EncodedLut kTILING6_1_2{{48, 27, 0}, 2, + "AQwDDAoDBgMKAwYIBQgGCAUMDAkIAQkMDAUKAQwDAQsMCwEGCQYBBgkHDAcJCQgMDAgDCwcMBAwABAEMAQQKBwoECgcCDAIHBwMMDAMAAQIMBgwC" + "BgMMAwYIBQgGCAUADAAFBQEMDAECAwAMAAwCDAkCBQIJAgULBAsFCwQMDAgLAAgMDAQJAAwCAAoMCgAFCAUABQgGDAYICAsMDAsCCgYMBAwADAUA" + "CgAFAAoDBgMKAwYMDAcDBAcMDAYFBAwGDAgGAwYIBgMKAAoDCgAMDAkKBAkMDAAIBQwHBQgMCAUACgAFAAoDDAMKCgsMDAsHCAMMAgwAAggMCAIH" + "CgcCBwoEDAQKCgkMDAkACAQMAgwADAsABwALAAcJBgkHCQYMDAoJAgoMDAYLBQwBBQIMAgULBAsFCwQDDAMEBAAMDAABAgMMBwwDBwAMAAcJBgkH" + "CQYBDAEGBgIMDAIDAAEMBgwEBgkMCQYBCwEGAQsADAALCwgMDAgECQAMBQwBDAYBCwEGAQsABwALAAcMDAQABQQMDAcGBQwHDAkHAAcJBwALAQsA" + "CwEMDAoLBQoMDAEJAwwBDAgBBAEIAQQKBwoECgcMDAsKAwsMDAcIAwwBAwkMCQMECwQDBAsFDAULCwoMDAoBCQUMBwwFBwoMCgcCCAIHAggBDAEI" + "CAkMDAkFCgEMBgwCDAcCCAIHAggBBAEIAQQMDAUBBgUMDAQHBgwEDAoEAQQKBAEIAggBCAIMDAsIBgsMDAIKBwwFDAsFAgULBQIJAwkCCQMMDAgJ" + "BwgMDAMLBAwGBAsMCwQDCQMEAwkCDAIJCQoMDAoGCwIMBwwDDAQDCQMEAwkCBQIJAgUMDAYCBwYMDAUEAwwHAwQMBAMJAgkDCQIFDAUCAgYMDAYH" + "BAUMBgwEDAsEAwQLBAMJAgkDCQIMDAoJBgoMDAILBQwHBQsMCwUCCQIFAgkDDAMJCQgMDAgHCwMMBAwGBAoMCgQBCAEEAQgCDAIICAsMDAsGCgIM" + "AgwGAgcMBwIIAQgCCAEEDAQBAQUMDAUGBwQMBQwHDAoHAgcKBwIIAQgCCAEMDAkIBQkMDAEKAQwDDAkDBAMJAwQLBQsECwUMDAoLAQoMDAUJAQwD" + "AQgMCAEECgQBBAoHDAcKCgsMDAsDCAcMBwwFBwkMCQcACwAHAAsBDAELCwoMDAoFCQEMAQwFAQYMBgELAAsBCwAHDAcAAAQMDAQFBgcMBAwGDAkG" + "AQYJBgELAAsBCwAMDAgLBAgMDAAJAwwHDAAHCQcABwkGAQYJBgEMDAIGAwIMDAEAAQwFDAIFCwUCBQsEAwQLBAMMDAAEAQAMDAMCAAwCAAsMCwAH" + "CQcABwkGDAYJCQoMDAoCCwYMAAwCDAgCBwIIAgcKBAoHCgQMDAkKAAkMDAQIBwwFDAgFAAUIBQAKAwoACgMMDAsKBwsMDAMIBgwEBggMCAYDCgMG" + "AwoADAAKCgkMDAkECAAMAAwEAAUMBQAKAwoACgMGDAYDAwcMDAcEBQYMAgwADAoABQAKAAUIBggFCAYMDAsIAgsMDAYKAgwAAgkMCQIFCwUCBQsE" + "DAQLCwgMDAgACQQMAgwGDAMGCAYDBggFAAUIBQAMDAEFAgEMDAADAAwEDAEECgQBBAoHAgcKBwIMDAMHAAMMDAIBAwwBDAsBBgELAQYJBwkGCQcM" + "DAgJAwgMDAcLAwwBAwoMCgMGCAYDBggFDAUICAkMDAkBCgUM" +}; + +inline constexpr EncodedLut kTILING6_2{{48, 15, 0}, 2, + "AQoDBgMKAwYIBQgGCAUJAQsDCwEGCQYBBgkHCAcJBAEAAQQKBwoECgcCAwIHBgMCAwYIBQgGCAUAAQAFAAkCBQIJAgULBAsFCwQIAAoCCgAFCAUA" + "BQgGCwYIBAUACgAFAAoDBgMKAwYHBAgGAwYIBgMKAAoDCgAJBQgHCAUACgAFAAoDCwMKAggACAIHCgcCBwoECQQKAgsABwALAAcJBgkHCQYKBQIB" + "AgULBAsFCwQDAAMEBwADAAcJBgkHCQYBAgEGBgkECQYBCwEGAQsACAALBQYBCwEGAQsABwALAAcEBQkHAAcJBwALAQsACwEKAwgBBAEIAQQKBwoE" + "CgcLAwkBCQMECwQDBAsFCgULBwoFCgcCCAIHAggBCQEIBgcCCAIHAggBBAEIAQQFBgoEAQQKBAEIAggBCAILBwsFAgULBQIJAwkCCQMIBAsGCwQD" + "CQMEAwkCCgIJBwQDCQMEAwkCBQIJAgUGAwQHBAMJAgkDCQIFBgUCBgsEAwQLBAMJAgkDCQIKBQsHCwUCCQIFAgkDCAMJBAoGCgQBCAEEAQgCCwII" + "AgcGBwIIAQgCCAEEBQQBBQoHAgcKBwIIAQgCCAEJAQkDBAMJAwQLBQsECwUKAQgDCAEECgQBBAoHCwcKBwkFCQcACwAHAAsBCgELAQYFBgELAAsB" + "CwAHBAcABAkGAQYJBgELAAsBCwAIAwAHCQcABwkGAQYJBgECAQIFCwUCBQsEAwQLBAMAAAsCCwAHCQcABwkGCgYJAAgCBwIIAgcKBAoHCgQJBwgF" + "AAUIBQAKAwoACgMLBggECAYDCgMGAwoACQAKAAUEBQAKAwoACgMGBwYDAgoABQAKAAUIBggFCAYLAgkACQIFCwUCBQsECAQLAgMGCAYDBggFAAUI" + "BQABAAEECgQBBAoHAgcKBwIDAwsBBgELAQYJBwkGCQcIAwoBCgMGCAYDBggFCQUI" +}; + +inline constexpr EncodedLut kTILING7_1{{16, 9, 0}, 2, + "CQUECgECCAMACwcGCAMACgECAwAIBQQJBwYLCAQHCQABCwIDCgYFCwIDCQABAAEJBgUKBAcIAQIKBwYLBQQJAgMLBAcIBgUKCwMCCAcECgUGCgIB" + "CwYHCQQFCQEACgUGCAcEBQYKAwILAQAJBwQIAQAJAwILCAADCQQFCwYHBgcLAAMIAgEKBAUJAgEKAAMI" +}; + +inline constexpr EncodedLut kTILING7_2{{16, 3, 15}, 3, + "AQIKAwQIBAMFAAUDBQAJAwAICQEEAgQBBAIFCgUCCQUEAAoBCgAICggCAwIIAwAIAQYKBgEHAgcBBwILAQIKCwMGAAYDBgAHCAcACwcGAggDCAIK" + "CAoAAQAKCQUECwMGAAYDBgAHCAcACwcGAwQIBAMFAAUDBQAJAwAIBAkHCwcJBQsJCwUGAAEJAgcLBwIEAwQCBAMIAgMLCAAHAQcABwEECQQBCAQH" + "AwkACQMLCQsBAgELAgMLAAUJBQAGAQYABgEKAAEJCgIFAwUCBQMGCwYDBgUKAQsCCwEJCwkDAAMJBgUKCAAHAQcABwEECQQBCAQHAAUJBQAGAQYA" + "BgEKAAEJBQoECAQKBggKCAYHCwcGCQEEAgQBBAIFCgUCCQUEAQYKBgEHAgcBBwILAQIKBgsFCQULBwkLCQcECAQHCgIFAwUCBQMGCwYDBgUKAgcL" + "BwIEAwQCBAMIAgMLBwgGCgYIBAoICgQFBwQIBQIKAgUDBgMFAwYLCgUGCwcCBAIHAgQDCAMECwMCBggHCAYKCAoEBQQKBgcLBAEJAQQCBQIEAgUK" + "BAUJCgYBBwEGAQcCCwIHCgIBBQsGCwUJCwkHBAcJCgUGBwAIAAcBBAEHAQQJBwQICQUABgAFAAYBCgEGCQEABAoFCgQICggGBwYICwMCCQUABgAF" + "AAYBCgEGCQEABQIKAgUDBgMFAwYLCgUGAgsBCQELAwkLCQMACQEACwcCBAIHAgQDCAMECwMCBwAIAAcBBAEHAQQJBwQIAAkDCwMJAQsJCwECBAUJ" + "BgMLAwYABwAGAAcIBgcLCAQDBQMEAwUACQAFCAADBwkECQcLCQsFBgULCAADCgYBBwEGAQcCCwIHCgIBBgMLAwYABwAGAAcIBgcLAwgCCgIIAAoI" + "CgABCgIBCAQDBQMEAwUACQAFCAADBAEJAQQCBQIEAgUKBAUJAQoACAAKAggKCAID" +}; + +inline constexpr EncodedLut kTILING7_3{{16, 3, 27}, 3, + "DAIKDAoFDAUEDAQIDAgDDAMADAAJDAkBDAECDAUEDAQIDAgDDAMCDAIKDAoBDAEADAAJDAkFBQQMCgUMAgoMAwIMCAMMAAgMAQAMCQEMBAkMDAAI" + "DAgHDAcGDAYKDAoBDAECDAILDAsDDAMADAcGDAYKDAoBDAEADAAIDAgDDAMCDAILDAsHBwYMCAcMAAgMAQAMCgEMAgoMAwIMCwMMBgsMCQUMAAkM" + "AwAMCwMMBgsMBwYMCAcMBAgMBQQMAwAMCwMMBgsMBQYMCQUMBAkMBwQMCAcMAAgMDAMADAAJDAkFDAUGDAYLDAsHDAcEDAQIDAgDDAEJDAkEDAQH" + "DAcLDAsCDAIDDAMIDAgADAABDAQHDAcLDAsCDAIBDAEJDAkADAADDAMIDAgEBAcMCQQMAQkMAgEMCwIMAwsMAAMMCAAMBwgMDAMLDAsGDAYFDAUJ" + "DAkADAABDAEKDAoCDAIDDAYFDAUJDAkADAADDAMLDAsCDAIBDAEKDAoGBgUMCwYMAwsMAAMMCQAMAQkMAgEMCgIMBQoMCgYMAQoMAAEMCAAMBwgM" + "BAcMCQQMBQkMBgUMAAEMCAAMBwgMBgcMCgYMBQoMBAUMCQQMAQkMDAABDAEKDAoGDAYHDAcIDAgEDAQFDAUJDAkACwcMAgsMAQIMCQEMBAkMBQQM" + "CgUMBgoMBwYMAQIMCQEMBAkMBwQMCwcMBgsMBQYMCgUMAgoMDAECDAILDAsHDAcEDAQJDAkFDAUGDAYKDAoBCAQMAwgMAgMMCgIMBQoMBgUMCwYM" + "BwsMBAcMAgMMCgIMBQoMBAUMCAQMBwgMBgcMCwYMAwsMDAIDDAMIDAgEDAQFDAUKDAoGDAYHDAcLDAsCDAQIDAgDDAMCDAIKDAoFDAUGDAYLDAsH" + "DAcEDAMCDAIKDAoFDAUEDAQIDAgHDAcGDAYLDAsDAwIMCAMMBAgMBQQMCgUMBgoMBwYMCwcMAgsMDAcLDAsCDAIBDAEJDAkEDAQFDAUKDAoGDAYH" + "DAIBDAEJDAkEDAQHDAcLDAsGDAYFDAUKDAoCAgEMCwIMBwsMBAcMCQQMBQkMBgUMCgYMAQoMDAYKDAoBDAEADAAIDAgHDAcEDAQJDAkFDAUGDAEA" + "DAAIDAgHDAcGDAYKDAoFDAUEDAQJDAkBAQAMCgEMBgoMBwYMCAcMBAgMBQQMCQUMAAkMCwMMBgsMBQYMCQUMAAkMAQAMCgEMAgoMAwIMBQYMCQUM" + "AAkMAwAMCwMMAgsMAQIMCgEMBgoMDAUGDAYLDAsDDAMADAAJDAkBDAECDAIKDAoFCQEMBAkMBwQMCwcMAgsMAwIMCAMMAAgMAQAMBwQMCwcMAgsM" + "AQIMCQEMAAkMAwAMCAMMBAgMDAcEDAQJDAkBDAECDAILDAsDDAMADAAIDAgHDAUJDAkADAADDAMLDAsGDAYHDAcIDAgEDAQFDAADDAMLDAsGDAYF" + "DAUJDAkEDAQHDAcIDAgAAAMMCQAMBQkMBgUMCwYMBwsMBAcMCAQMAwgMCAAMBwgMBgcMCgYMAQoMAgEMCwIMAwsMAAMMBgcMCgYMAQoMAAEMCAAM" + "AwgMAgMMCwIMBwsMDAYHDAcIDAgADAABDAEKDAoCDAIDDAMLDAsGCgIMBQoMBAUMCAQMAwgMAAMMCQAMAQkMAgEMBAUMCAQMAwgMAgMMCgIMAQoM" + "AAEMCQAMBQkMDAQFDAUKDAoCDAIDDAMIDAgADAABDAEJDAkE" +}; + +inline constexpr EncodedLut kTILING7_4_1{{16, 15, 0}, 2, + "AwQIBAMKAgoDBAoFCQEAAQYKBgEIAAgBBggHCwMCCwMGCQYDBgkFAAkDBwQIAgcLBwIJAQkCBwkECAADAAUJBQALAwsABQsGCgIBCAAHCgcABwoG" + "AQoABAUJCQEECwQBBAsHAgsBBQYKCgIFCAUCBQgEAwgCBgcLBQIKAgUIBAgFAggDCwcGBAEJAQQLBwsEAQsCCgYFBwAIAAcKBgoHAAoBCQUECQUA" + "CwAFAAsDBgsFAQIKCwcCCQIHAgkBBAkHAwAIBgMLAwYJBQkGAwkACAQHCgYBCAEGAQgABwgGAgMLCAQDCgMEAwoCBQoEAAEJ" +}; + +inline constexpr EncodedLut kTILING7_4_2{{16, 27, 0}, 2, + "CQQIBAkFCgUJAQoJCgECAAIBAgADCAMACQgACwYKBgsHCAcLAwgLCAMAAgADAAIBCgECCwoCCwMIAAgDCAAJCAkEBQQJBAUHBgcFBwYLBwsICAcL" + "BwgECQQIAAkICQABAwEAAQMCCwIDCAsDCgUJBQoGCwYKAgsKCwIDAQMCAwEACQABCgkBCAAJAQkACQEKCQoFBgUKBQYEBwQGBAcIBAgJCQEKAgoB" + "CgILCgsGBwYLBgcFBAUHBQQJBQkKCgILAwsCCwMICwgHBAcIBwQGBQYEBgUKBgoLCwIKAgsDCAMLBwgLCAcEBgQHBAYFCgUGCwoGCgEJAQoCCwIK" + "BgsKCwYHBQcGBwUECQQFCgkFCQAIAAkBCgEJBQoJCgUGBAYFBgQHCAcECQgECQUKBgoFCgYLCgsCAwILAgMBAAEDAQAJAQkKCwcIBAgHCAQJCAkA" + "AQAJAAEDAgMBAwILAwsICAMLAwgACQAIBAkICQQFBwUEBQcGCwYHCAsHCgYLBwsGCwcICwgDAAMIAwACAQIAAgEKAgoLCAQJBQkECQUKCQoBAgEK" + "AQIAAwACAAMIAAgJ" +}; + +inline constexpr EncodedLut kTILING8{{6, 6, 0}, 2, + "CQgKCggLAQUDAwUHAAQCBAYCAAIEBAIGAQMFAwcFCQoICgsI" +}; + +inline constexpr EncodedLut kTILING9{{8, 12, 0}, 2, + "AgoFAwIFAwUEAwQIBAcLCQQLCQsCCQIBCgcGAQcKAQgHAQAIAwYLAAYDAAUGAAkFAwsGAAMGAAYFAAUJCgYHAQoHAQcIAQgABAsHCQsECQILCQEC" + "AgUKAwUCAwQFAwgE" +}; + +inline constexpr EncodedLut kTILING10_1_1{{6, 12, 0}, 2, + "BQoHCwcKCAEJAQgDAQIFBgUCBAMAAwQHCwAIAAsCBAkGCgYJCQAKAgoABggECAYLBwIDAgcGAAEEBQQBBwkFCQcICgELAwsB" +}; + +inline constexpr EncodedLut kTILING10_1_1_{{6, 12, 0}, 2, + "BQkHCAcJCwEKAQsDAwIHBgcCBAEAAQQFCgAJAAoCBAgGCwYICAALAgsABgkECQYKBQIBAgUGAAMEBwQDBwoFCgcLCQEIAwgB" +}; + +inline constexpr EncodedLut kTILING10_1_2{{6, 24, 0}, 2, + "AwsHAwcICQgHBQkHCQUKCQoBAwEKCwMKBwYFBwUEAAQFAQAFAAECAAIDBwMCBgcCCwIKBgsKCwYECwQIAAgECQAEAAkKAAoCCwIKCwoGBAYKCQQK" + "BAkABAAICwgAAgsABwYFBAcFBwQABwADAgMAAQIAAgEFAgUGBwgDCwcDBwsKBwoFCQUKAQkKCQEDCQMI" +}; + +inline constexpr EncodedLut kTILING10_2{{6, 24, 0}, 2, + "DAUJDAkIDAgDDAMBDAEKDAoLDAsHDAcFDAEADAAEDAQHDAcDDAMCDAIGDAYFDAUBBAgMBgQMCgYMCQoMAAkMAgAMCwIMCAsMDAkEDAQGDAYLDAsI" + "DAgADAACDAIKDAoJAAMMBAAMBQQMAQUMAgEMBgIMBwYMAwcMCgUMCwoMAwsMAQMMCQEMCAkMBwgMBQcM" +}; + +inline constexpr EncodedLut kTILING10_2_{{6, 24, 0}, 2, + "CAcMCQgMAQkMAwEMCwMMCgsMBQoMBwUMBAUMAAQMAwAMBwMMBgcMAgYMAQIMBQEMDAsGDAYEDAQJDAkKDAoCDAIADAAIDAgLBgoMBAYMCAQMCwgM" + "AgsMAAIMCQAMCgkMDAcEDAQADAABDAEFDAUGDAYCDAIDDAMHDAcLDAsKDAoBDAEDDAMIDAgJDAkFDAUH" +}; + +inline constexpr EncodedLut kTILING11{{12, 12, 0}, 2, + "AgoJAgkHAgcDBwkEAQYCAQgGAQkICAcGCAMBCAEGCAYEBgEKAAgLAAsFAAUBBQsGCQUHCQcCCQIAAgcLBQAEBQsABQoLCwMABQQABQALBQsKCwAD" + "CQcFCQIHCQACAgsHAAsIAAULAAEFBQYLCAEDCAYBCAQGBgoBAQIGAQYIAQgJCAYHAgkKAgcJAgMHBwQJ" +}; + +inline constexpr EncodedLut kTILING12_1_1{{24, 12, 0}, 2, + "BwYLCgMCAwoICQgKBgUKCQIBAgkLCAsJCgYFBwkECQcBAwEHBwYLBAgFAwUIBQMBBQQJCAEAAQgKCwoIAQIKAAkDBQMJAwUHCgECAAsDCwAGBAYA" + "CAMAAgkBCQIEBgQCAwAIAgsBBwELAQcFBgUKBwsEAgQLBAIACQUEBggHCAYAAgAGCAMABwQLCQsECwkKBAcICwADAAsJCgkLBAcIBQkGAAYJBgAC" + "CwcGBAoFCgQCAAIECwIDAQgACAEHBQcBAAEJAwgCBAIIAgQGAgMLAQoABgAKAAYECQABAwoCCgMFBwUDCQABBAUICggFCAoLCAQHBQsGCwUDAQMF" + "BQQJBgoHAQcKBwEDCgECBQYJCwkGCQsICwIDBgcKCAoHCggJ" +}; + +inline constexpr EncodedLut kTILING12_1_1_{{24, 12, 0}, 2, + "AwILCgcGBwoICQgKAgEKCQYFBgkLCAsJCQQFBwoGCgcBAwEHBwQIBgsFAwULBQMBAQAJCAUEBQgKCwoIAQAJAgoDBQMKAwUHCwMCAAoBCgAGBAYA" + "CQEAAggDCAIEBgQCAwILAAgBBwEIAQcFBgcLBQoEAgQKBAIACAcEBgkFCQYAAgAGCAcEAwALCQsACwkKAAMICwQHBAsJCgkLBAUJBwgGAAYIBgAC" + "CgUGBAsHCwQCAAIECAADAQsCCwEHBQcBAAMIAQkCBAIJAgQGAgEKAwsABgALAAYECgIBAwkACQMFBwUDCQQFAAEICggBCAoLCwYHBQgECAUDAQMF" + "BQYKBAkHAQcJBwEDCgUGAQIJCwkCCQsICwYHAgMKCAoDCggJ" +}; + +inline constexpr EncodedLut kTILING12_1_2{{24, 24, 0}, 2, + "BwMLAwcICQgHBgkHCQYKAgoGCwIGAgsDBgIKAgYLCAsGBQgGCAUJAQkFCgEFAQoCCgkFCQoBAwEKBgMKAwYHBAcGBQQGBAUJBwgLAwsICwMBCwEG" + "BQYBBgUEBgQHCAcEBQEJAQUKCwoFBAsFCwQIAAgECQAEAAkBAQkKBQoJCgUHCgcCAwIHAgMAAgABCQEACgsCCwoGBAYKAQQKBAEAAwABAgMBAwIL" + "CAkACQgEBgQIAwYIBgMCAQIDAAEDAQAJAwsIBwgLCAcFCAUAAQAFAAECAAIDCwMCBgsKAgoLCgIACgAFBAUABQQHBQcGCwYHCQgECAkAAgAJBQIJ" + "AgUGBwYFBAcFBwQICAQACQAEAAkKAAoDCwMKAwsHAwcIBAgHBAAIAAQJCgkEBwoECgcLAwsHCAMHAwgABAkIAAgJCAACCAIHBgcCBwYFBwUECQQF" + "CwoGCgsCAAILBwALAAcEBQQHBgUHBQYKCwgDCAsHBQcLAgULBQIBAAECAwACAAMIAAgJBAkICQQGCQYBAgEGAQIDAQMACAADAgoLBgsKCwYECwQD" + "AAMEAwABAwECCgIBCQoBCgkFBwUJAAcJBwADAgMAAQIAAgEKCQUBCgEFAQoLAQsACAALAAgEAAQJBQkECAsHCwgDAQMIBAEIAQQFBgUEBwYEBgcL" + "BQoJAQkKCQEDCQMEBwQDBAcGBAYFCgUGCgYCCwIGAgsIAggBCQEIAQkFAQUKBgoFCwcDCAMHAwgJAwkCCgIJAgoGAgYLBwsG" +}; + +inline constexpr EncodedLut kTILING12_2{{24, 24, 0}, 2, + "CQgMCgkMAgoMAwIMCwMMBgsMBwYMCAcMCAsMCQgMAQkMAgEMCgIMBQoMBgUMCwYMAwEMBwMMBAcMCQQMBQkMBgUMCgYMAQoMDAMBDAEFDAUGDAYL" + "DAsHDAcEDAQIDAgDCwoMCAsMAAgMAQAMCQEMBAkMBQQMCgUMDAUHDAcDDAMCDAIKDAoBDAEADAAJDAkFBAYMAAQMAQAMCgEMAgoMAwIMCwMMBgsM" + "BgQMAgYMAwIMCAMMAAgMAQAMCQEMBAkMDAcFDAUBDAEADAAIDAgDDAMCDAILDAsHDAIADAAEDAQFDAUKDAoGDAYHDAcLDAsCAgAMBgIMBwYMCAcM" + "BAgMBQQMCQUMAAkMDAkKDAoLDAsHDAcEDAQIDAgDDAMADAAJCgkMCwoMBwsMBAcMCAQMAwgMAAMMCQAMDAACDAIGDAYHDAcIDAgEDAQFDAUJDAkA" + "AAIMBAAMBQQMCgUMBgoMBwYMCwcMAgsMBQcMAQUMAAEMCAAMAwgMAgMMCwIMBwsMDAQGDAYCDAIDDAMIDAgADAABDAEJDAkEDAYEDAQADAABDAEK" + "DAoCDAIDDAMLDAsGBwUMAwcMAgMMCgIMAQoMAAEMCQAMBQkMDAoLDAsIDAgADAABDAEJDAkEDAQFDAUKAQMMBQEMBgUMCwYMBwsMBAcMCAQMAwgM" + "DAEDDAMHDAcEDAQJDAkFDAUGDAYKDAoBDAsIDAgJDAkBDAECDAIKDAoFDAUGDAYLDAgJDAkKDAoCDAIDDAMLDAsGDAYHDAcI" +}; + +inline constexpr EncodedLut kTILING12_2_{{24, 24, 0}, 2, + "DAILDAsHDAcGDAYKDAoJDAkIDAgDDAMCDAEKDAoGDAYFDAUJDAkIDAgLDAsCDAIBDAQFDAUKDAoGDAYHDAcDDAMBDAEJDAkEBwYMCAcMBAgMBQQM" + "AQUMAwEMCwMMBgsMDAAJDAkFDAUEDAQIDAgLDAsKDAoBDAEAAQIMCQEMAAkMAwAMBwMMBQcMCgUMAgoMDAECDAILDAsDDAMADAAEDAQGDAYKDAoB" + "DAMADAAJDAkBDAECDAIGDAYEDAQIDAgDAwAMCwMMAgsMAQIMBQEMBwUMCAcMAAgMBgUMCwYMBwsMBAcMAAQMAgAMCgIMBQoMDAcEDAQJDAkFDAUG" + "DAYCDAIADAAIDAgHCAcMAAgMAwAMCwMMCgsMCQoMBAkMBwQMDAcIDAgADAADDAMLDAsKDAoJDAkEDAQHBAcMCQQMBQkMBgUMAgYMAAIMCAAMBwgM" + "DAUGDAYLDAsHDAcEDAQADAACDAIKDAoFDAADDAMLDAsCDAIBDAEFDAUHDAcIDAgAAAMMCQAMAQkMAgEMBgIMBAYMCAQMAwgMAgEMCwIMAwsMAAMM" + "BAAMBgQMCgYMAQoMDAIBDAEJDAkADAADDAMHDAcFDAUKDAoCCQAMBQkMBAUMCAQMCwgMCgsMAQoMAAEMDAYHDAcIDAgEDAQFDAUBDAEDDAMLDAsG" + "BQQMCgUMBgoMBwYMAwcMAQMMCQEMBAkMCgEMBgoMBQYMCQUMCAkMCwgMAgsMAQIMCwIMBwsMBgcMCgYMCQoMCAkMAwgMAgMM" +}; + +inline constexpr EncodedLut kTILING13_1{{2, 12, 0}, 2, + "CwcGAQIKCAMACQUECAQHAgMLCQABCgYF" +}; + +inline constexpr EncodedLut kTILING13_1_{{2, 12, 0}, 2, + "BwQICwMCAQAJBQYKBgcLCgIBAAMIBAUJ" +}; + +inline constexpr EncodedLut kTILING13_2{{2, 6, 18}, 3, + "AQIKCwcGAwQIBAMFAAUDBQAJCAMACwcGCQEEAgQBBAIFCgUCCQUECAMAAQYKBgEHAgcBBwILCQUEAQIKCwMGAAYDBgAHCAcACQUECwcGAAoBCgAI" + "CggCAwIIAQIKAwAIBAkHCwcJBQsJCwUGAgMLCAQHAAUJBQAGAQYABgEKCQABCAQHCgIFAwUCBQMGCwYDBgUKCQABAgcLBwIEAwQCBAMIBgUKAgML" + "CAAHAQcABwEECQQBBgUKCAQHAQsCCwEJCwkDAAMJAgMLAAEJBQoECAQKBggKCAYH" +}; + +inline constexpr EncodedLut kTILING13_2_{{2, 6, 18}, 3, + "CgUGCwMCBwAIAAcBBAEHAQQJCwMCBwQICQUABgAFAAYBCgEGAQAJBwQIBQIKAgUDBgMFAwYLCgUGAQAJCwcCBAIHAgQDCAMECgUGBwQIAgsBCQEL" + "AwkLCQMACwMCCQEABAoFCgQICggGBwYIBgcLCAADBAEJAQQCBQIEAgUKCAADBAUJCgYBBwEGAQcCCwIHAgEKBAUJBgMLAwYABwAGAAcIBgcLAgEK" + "CAQDBQMEAwUACQAFBgcLBAUJAwgCCgIIAAoICgABCAADCgIBBQsGCwUJCwkHBAcJ" +}; + +inline constexpr EncodedLut kTILING13_3{{2, 12, 30}, 3, + "CwcGDAIKDAoFDAUEDAQIDAgDDAMADAAJDAkBDAECAQIKCQUMAAkMAwAMCwMMBgsMBwYMCAcMBAgMBQQMCwcGDAUEDAQIDAgDDAMCDAIKDAoBDAEA" + "DAAJDAkFAQIKDAMADAAJDAkFDAUGDAYLDAsHDAcEDAQIDAgDCAMACwcMAgsMAQIMCQEMBAkMBQQMCgUMBgoMBwYMCwcGBQQMCgUMAgoMAwIMCAMM" + "AAgMAQAMCQEMBAkMCAMAAQIMCQEMBAkMBwQMCwcMBgsMBQYMCgUMAgoMCQUEDAAIDAgHDAcGDAYKDAoBDAECDAILDAsDDAMACQUEDAcGDAYKDAoB" + "DAEADAAIDAgDDAMCDAILDAsHCAMADAECDAILDAsHDAcEDAQJDAkFDAUGDAYKDAoBCQUEBwYMCAcMAAgMAQAMCgEMAgoMAwIMCwMMBgsMAQIKAwAM" + "CwMMBgsMBQYMCQUMBAkMBwQMCAcMAAgMCAQHDAMLDAsGDAYFDAUJDAkADAABDAEKDAoCDAIDAgMLCgYMAQoMAAEMCAAMBwgMBAcMCQQMBQkMBgUM" + "CAQHDAYFDAUJDAkADAADDAMLDAsCDAIBDAEKDAoGAgMLDAABDAEKDAoGDAYHDAcIDAgEDAQFDAUJDAkAAAEJCAQMAwgMAgMMCgIMBQoMBgUMCwYM" + "BwsMBAcMCAQHBgUMCwYMAwsMAAMMCQAMAQkMAgEMCgIMBQoMCQABAgMMCgIMBQoMBAUMCAQMBwgMBgcMCwYMAwsMBgUKDAEJDAkEDAQHDAcLDAsC" + "DAIDDAMIDAgADAABBgUKDAQHDAcLDAsCDAIBDAEJDAkADAADDAMIDAgECQABDAIDDAMIDAgEDAQFDAUKDAoGDAYHDAcLDAsCBgUKBAcMCQQMAQkM" + "AgEMCwIMAwsMAAMMCAAMBwgMAgMLAAEMCAAMBwgMBgcMCgYMBQoMBAUMCQQMAQkM" +}; + +inline constexpr EncodedLut kTILING13_3_{{2, 12, 30}, 3, + "AwILCAcMAAgMAQAMCgEMBgoMBQYMCQUMBAkMBwQMBQYKDAILDAsHDAcEDAQJDAkBDAEADAAIDAgDDAMCCgUGDAcEDAQJDAkBDAECDAILDAsDDAMA" + "DAAIDAgHCwMCDAEADAAIDAgHDAcGDAYKDAoFDAUEDAQJDAkBBwQICwMMBgsMBQYMCQUMAAkMAQAMCgEMAgoMAwIMBwQIBQYMCQUMAAkMAwAMCwMM" + "AgsMAQIMCgEMBgoMCwMCAQAMCgEMBgoMBwYMCAcMBAgMBQQMCQUMAAkMAQAJDAQIDAgDDAMCDAIKDAoFDAUGDAYLDAsHDAcEBwQIDAUGDAYLDAsD" + "DAMADAAJDAkBDAECDAIKDAoFAQAJDAMCDAIKDAoFDAUEDAQIDAgHDAcGDAYLDAsDCgUGBwQMCwcMAgsMAQIMCQEMAAkMAwAMCAMMBAgMCQEAAwIM" + "CAMMBAgMBQQMCgUMBgoMBwYMCwcMAgsMAAMICQQMAQkMAgEMCwIMBwsMBgcMCgYMBQoMBAUMCwYHDAMIDAgEDAQFDAUKDAoCDAIBDAEJDAkADAAD" + "BgcLDAQFDAUKDAoCDAIDDAMIDAgADAABDAEJDAkECAADDAIBDAEJDAkEDAQHDAcLDAsGDAYFDAUKDAoCBAUJCAAMBwgMBgcMCgYMAQoMAgEMCwIM" + "AwsMAAMMBAUJBgcMCgYMAQoMAAEMCAAMAwgMAgMMCwIMBwsMCAADAgEMCwIMBwsMBAcMCQQMBQkMBgUMCgYMAQoMAgEKDAUJDAkADAADDAMLDAsG" + "DAYHDAcIDAgEDAQFBAUJDAYHDAcIDAgADAABDAEKDAoCDAIDDAMLDAsGAgEKDAADDAMLDAsGDAYFDAUJDAkEDAQHDAcIDAgABgcLBAUMCAQMAwgM" + "AgMMCgIMAQoMAAEMCQAMBQkMCgIBAAMMCQAMBQkMBgUMCwYMBwsMBAcMCAQMAwgM" +}; + +inline constexpr EncodedLut kTILING13_4{{2, 4, 36}, 3, + "DAIKDAoFDAUGDAYLDAsHDAcEDAQIDAgDDAMADAAJDAkBDAECCwMMBgsMBwYMCAcMBAgMBQQMCQUMAAkMAQAMCgEMAgoMAwIMCQEMBAkMBQQMCgUM" + "BgoMBwYMCwcMAgsMAwIMCAMMAAgMAQAMDAAIDAgHDAcEDAQJDAkFDAUGDAYKDAoBDAECDAILDAsDDAMADAMLDAsGDAYHDAcIDAgEDAQFDAUJDAkA" + "DAABDAEKDAoCDAIDCAAMBwgMBAcMCQQMBQkMBgUMCgYMAQoMAgEMCwIMAwsMAAMMCgIMBQoMBgUMCwYMBwsMBAcMCAQMAwgMAAMMCQAMAQkMAgEM" + "DAEJDAkEDAQFDAUKDAoGDAYHDAcLDAsCDAIDDAMIDAgADAAB" +}; + +inline constexpr EncodedLut kTILING13_5_1{{2, 4, 18}, 3, + "BwYLAQAJCgMCAwoFAwUIBAgFAQIKBwQIAwALBgsACQYABgkFAwAIBQYKAQIJBAkCCwQCBAsHBQQJAwILCAEAAQgHAQcKBgoHBAcIAgEKCwADAAsG" + "AAYJBQkGAgMLBAUJAAEIBwgBCgcBBwoGAAEJBgcLAgMKBQoDCAUDBQgEBgUKAAMICQIBAgkEAgQLBwsE" +}; + +inline constexpr EncodedLut kTILING13_5_2{{2, 4, 30}, 3, + "AQAJBwQIBwgDBwMLAgsDCwIKCwoGBQYKBgUHBAcFBwQICwMCBgsCCgYCBgoFCQUKAQkKCQEAAgABAAIDBQYKCQEABAkACAQABAgHCwcIAwsICwMC" + "AAIDAgABAwILBQYKBQoBBQEJAAkBCQAICQgEBAgHBAcFBgUHAgEKBAUJBAkABAAIAwgACAMLCAsHBgcLBwYEBQQGBAUJCAADBwgDCwcDBwsGCgYL" + "AgoLCgIBAwECAQMABgcLCgIBBQoBCQUBBQkECAQJAAgJCAADAQMAAwECAAMIBgcLBgsCBgIKAQoCCgEJCgkFBQkEBQQGBwYE" +}; + +inline constexpr EncodedLut kTILING14{{12, 12, 0}, 2, + "BQkIBQgCBQIGAwIIAgEFAgUIAggLBAgFCQQGCQYDCQMBCwMGAQsKAQQLAQAEBwsECAIACAUCCAcFCgIFAAcDAAoHAAkKBgcKAAMHAAcKAAoJBgoH" + "CAACCAIFCAUHCgUCAQoLAQsEAQQABwQLCQYECQMGCQEDCwYDAgUBAggFAgsIBAUIBQgJBQIIBQYCAwgC" +}; + +inline constexpr EncodedLut kTEST3{{24, 0, 0}, 1, + "BQEEBQECAgMEAwYG+vr9/P3+/v/7/P/7" +}; + +inline constexpr EncodedLut kTEST4{{8, 0, 0}, 1, + "BwcHB/n5+fk=" +}; + +inline constexpr EncodedLut kTEST6{{48, 3, 0}, 2, + "AgcKBAcLBQcBBQcDAQcJAwcKBgcFAQcIBAcIAQcIAwcLBQcCBQcAAQcJBgcGAgcJBAcIAgcJAgcKBgcHAwcKBAcLAwcLBgcE+vkE/fkL/PkL/fkK" + "+vkH/vkK/vkJ/PkI/vkJ+vkG//kJ+/kA+/kC/fkL//kI/PkI//kI+vkF/fkK//kJ+/kD+/kB/PkL/vkK" +}; + +inline constexpr EncodedLut kTEST7{{16, 5, 0}, 2, + "AQIFBwEDBAUHAwQBBgcEBAEFBwACAwUHAgECBgcFAgMGBwYDBAYHB/38+vkH/v36+Qb//vr5Bf79+/kC/P/7+QD8//r5BP38+/kD//77+QE=" +}; + +inline constexpr EncodedLut kTEST10{{6, 3, 0}, 2, + "AgQHBQYHAQMHAQMHBQYHAgQH" +}; + +inline constexpr EncodedLut kTEST12{{24, 4, 0}, 2, + "BAMHCwMCBwoCBgcFBgQHBwIBBwkFAgcBBQMHAgUBBwAFBAcDBgMHBgEGBwQBBAcIBAEHCAYBBwQDBgcGBAUHAwEFBwADBQcCAgUHAQECBwkEBgcH" + "BgIHBQIDBwoDBAcL" +}; + +inline constexpr EncodedLut kTEST13{{2, 7, 0}, 2, + "AQIDBAUGBwIDBAEFBgc=" +}; + +inline constexpr EncodedLut kSUBCONFIG13{{64, 0, 0}, 1, + "AAECBwP/C/8ECP//Dv///wUJDBcP/xUmERT/JBohHiwGCg0TEP8ZJRIY/yMWIB0r////Iv//HCr/H/8pGygnLQ==" +}; + +} // namespace bioimage_cpp::mesh::detail diff --git a/include/bioimage_cpp/mesh/marching_cubes.hxx b/include/bioimage_cpp/mesh/marching_cubes.hxx new file mode 100644 index 0000000..701082d --- /dev/null +++ b/include/bioimage_cpp/mesh/marching_cubes.hxx @@ -0,0 +1,836 @@ +#pragma once + +// Marching Cubes 33 / Lewiner implementation. +// +// The MC33 lookup tables and the case-selection structure below are derived +// from scikit-image 0.26.0 (BSD-3-Clause), whose implementation credits the +// original algorithm by Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and +// Geovan Tavares, "Efficient implementation of Marching Cubes' cases with +// topological guarantees", Journal of Graphics Tools 8(2), 2003. The +// scikit-image BSD-3-Clause notice must accompany redistributions of this +// derived table data. + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/profile.hxx" +#include "bioimage_cpp/mesh/detail/mc33_luts.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::mesh { + +enum class MarchingCubesMethod { + Lewiner, + Lorensen, +}; + +struct MarchingCubesResult { + // All vector-valued arrays are flat C-order arrays with a trailing size-3 + // component axis. Vertices/normals use the reference kernel's x/y/z order; + // the binding converts them to NumPy z/y/x order before returning them. + std::vector vertices; + std::vector faces; + std::vector normals; + std::vector values; +}; + +namespace detail::marching_cubes { + +constexpr double kEpsilon = std::numeric_limits::epsilon(); + +class Lut { +public: + explicit Lut(const detail::EncodedLut &encoded) + : shape_(encoded.shape), values_(decode(encoded.data)) {} + + [[nodiscard]] int get1(const int i0) const { + return values_[static_cast(i0)]; + } + + [[nodiscard]] int get2(const int i0, const int i1) const { + return values_[static_cast(i0 * shape_[1] + i1)]; + } + + [[nodiscard]] int get3(const int i0, const int i1, const int i2) const { + return values_[static_cast((i0 * shape_[1] + i1) * shape_[2] + i2)]; + } + +private: + static int base64_value(const char value) { + if (value >= 'A' && value <= 'Z') return value - 'A'; + if (value >= 'a' && value <= 'z') return value - 'a' + 26; + if (value >= '0' && value <= '9') return value - '0' + 52; + if (value == '+') return 62; + if (value == '/') return 63; + return -1; + } + + static std::vector decode(const std::string_view encoded) { + std::vector decoded; + decoded.reserve(encoded.size() * 3 / 4); + int accumulator = 0; + int bits = -8; + for (const char value : encoded) { + if (value == '=') break; + const int digit = base64_value(value); + if (digit < 0) continue; + accumulator = (accumulator << 6) | digit; + bits += 6; + if (bits >= 0) { + decoded.push_back(static_cast((accumulator >> bits) & 0xff)); + bits -= 8; + } + } + return decoded; + } + + std::array shape_{}; + std::vector values_; +}; + +struct Luts { + Lut cases_classic{detail::kCASESCLASSIC}; + Lut cases{detail::kCASES}; + Lut tiling1{detail::kTILING1}; + Lut tiling2{detail::kTILING2}; + Lut tiling3_1{detail::kTILING3_1}; + Lut tiling3_2{detail::kTILING3_2}; + Lut tiling4_1{detail::kTILING4_1}; + Lut tiling4_2{detail::kTILING4_2}; + Lut tiling5{detail::kTILING5}; + Lut tiling6_1_1{detail::kTILING6_1_1}; + Lut tiling6_1_2{detail::kTILING6_1_2}; + Lut tiling6_2{detail::kTILING6_2}; + Lut tiling7_1{detail::kTILING7_1}; + Lut tiling7_2{detail::kTILING7_2}; + Lut tiling7_3{detail::kTILING7_3}; + Lut tiling7_4_1{detail::kTILING7_4_1}; + Lut tiling7_4_2{detail::kTILING7_4_2}; + Lut tiling8{detail::kTILING8}; + Lut tiling9{detail::kTILING9}; + Lut tiling10_1_1{detail::kTILING10_1_1}; + Lut tiling10_1_1_alt{detail::kTILING10_1_1_}; + Lut tiling10_1_2{detail::kTILING10_1_2}; + Lut tiling10_2{detail::kTILING10_2}; + Lut tiling10_2_alt{detail::kTILING10_2_}; + Lut tiling11{detail::kTILING11}; + Lut tiling12_1_1{detail::kTILING12_1_1}; + Lut tiling12_1_1_alt{detail::kTILING12_1_1_}; + Lut tiling12_1_2{detail::kTILING12_1_2}; + Lut tiling12_2{detail::kTILING12_2}; + Lut tiling12_2_alt{detail::kTILING12_2_}; + Lut tiling13_1{detail::kTILING13_1}; + Lut tiling13_1_alt{detail::kTILING13_1_}; + Lut tiling13_2{detail::kTILING13_2}; + Lut tiling13_2_alt{detail::kTILING13_2_}; + Lut tiling13_3{detail::kTILING13_3}; + Lut tiling13_3_alt{detail::kTILING13_3_}; + Lut tiling13_4{detail::kTILING13_4}; + Lut tiling13_5_1{detail::kTILING13_5_1}; + Lut tiling13_5_2{detail::kTILING13_5_2}; + Lut tiling14{detail::kTILING14}; + Lut test3{detail::kTEST3}; + Lut test4{detail::kTEST4}; + Lut test6{detail::kTEST6}; + Lut test7{detail::kTEST7}; + Lut test10{detail::kTEST10}; + Lut test12{detail::kTEST12}; + Lut test13{detail::kTEST13}; + Lut subconfig13{detail::kSUBCONFIG13}; +}; + +inline const Luts &luts() { + static const Luts instance; + return instance; +} + +constexpr std::array, 12> kEdgeRelativeX{{ + {{0, 1}}, {{1, 1}}, {{1, 0}}, {{0, 0}}, {{0, 1}}, {{1, 1}}, + {{1, 0}}, {{0, 0}}, {{0, 0}}, {{1, 1}}, {{1, 1}}, {{0, 0}}, +}}; +constexpr std::array, 12> kEdgeRelativeY{{ + {{0, 0}}, {{0, 1}}, {{1, 1}}, {{1, 0}}, {{0, 0}}, {{0, 1}}, + {{1, 1}}, {{1, 0}}, {{0, 0}}, {{0, 0}}, {{1, 1}}, {{1, 1}}, +}}; +constexpr std::array, 12> kEdgeRelativeZ{{ + {{0, 0}}, {{0, 0}}, {{0, 0}}, {{0, 0}}, {{1, 1}}, {{1, 1}}, + {{1, 1}}, {{1, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, {{0, 1}}, +}}; + +class Cell { +public: + Cell(const int nx, const int ny) + : nx_(nx), ny_(ny), + face_layer1_(static_cast(nx) * static_cast(ny) * 4, -1), + face_layer2_(static_cast(nx) * static_cast(ny) * 4, -1) {} + + void new_z_value() { + face_layer1_.swap(face_layer2_); + std::fill(face_layer2_.begin(), face_layer2_.end(), -1); + } + + void set_cube( + const double isovalue, const int x, const int y, const int z, const int step, + const float v0, const float v1, const float v2, const float v3, + const float v4, const float v5, const float v6, const float v7 + ) { + x_ = x; + y_ = y; + z_ = z; + step_ = step; + v0_ = static_cast(v0) - isovalue; + v1_ = static_cast(v1) - isovalue; + v2_ = static_cast(v2) - isovalue; + v3_ = static_cast(v3) - isovalue; + v4_ = static_cast(v4) - isovalue; + v5_ = static_cast(v5) - isovalue; + v6_ = static_cast(v6) - isovalue; + v7_ = static_cast(v7) - isovalue; + index_ = (v0_ > 0.0 ? 1 : 0) | (v1_ > 0.0 ? 2 : 0) | (v2_ > 0.0 ? 4 : 0) + | (v3_ > 0.0 ? 8 : 0) | (v4_ > 0.0 ? 16 : 0) | (v5_ > 0.0 ? 32 : 0) + | (v6_ > 0.0 ? 64 : 0) | (v7_ > 0.0 ? 128 : 0); + center_calculated_ = false; + center_vertex_ = -1; + } + + [[nodiscard]] int index() const { return index_; } + [[nodiscard]] double v0() const { return v0_; } + [[nodiscard]] double v1() const { return v1_; } + [[nodiscard]] double v2() const { return v2_; } + [[nodiscard]] double v3() const { return v3_; } + [[nodiscard]] double v4() const { return v4_; } + [[nodiscard]] double v5() const { return v5_; } + [[nodiscard]] double v6() const { return v6_; } + [[nodiscard]] double v7() const { return v7_; } + + void add_triangles(const Lut &lut, const int lut_index, const int n_triangles) { + prepare_for_triangles(); + for (int triangle = 0; triangle < n_triangles; ++triangle) { + for (int corner = 0; corner < 3; ++corner) { + add_face_from_edge(lut.get2(lut_index, triangle * 3 + corner)); + } + } + } + + void add_triangles2(const Lut &lut, const int lut_index, const int lut_index2, const int n_triangles) { + prepare_for_triangles(); + for (int triangle = 0; triangle < n_triangles; ++triangle) { + for (int corner = 0; corner < 3; ++corner) { + add_face_from_edge(lut.get3(lut_index, lut_index2, triangle * 3 + corner)); + } + } + } + + [[nodiscard]] MarchingCubesResult take_result() { + for (std::size_t vertex = 0; vertex < values_.size(); ++vertex) { + const std::size_t base = vertex * 3; + const double length = std::sqrt( + static_cast(normals_[base]) * normals_[base] + + static_cast(normals_[base + 1]) * normals_[base + 1] + + static_cast(normals_[base + 2]) * normals_[base + 2] + ); + const double scale = length == 0.0 ? 0.0 : 1.0 / length; + normals_[base] = static_cast(normals_[base] * scale); + normals_[base + 1] = static_cast(normals_[base + 1] * scale); + normals_[base + 2] = static_cast(normals_[base + 2] * scale); + } + return { + std::move(vertices_), + std::move(faces_), + std::move(normals_), + std::move(values_), + }; + } + +private: + int add_vertex(const double x, const double y, const double z) { + if (values_.size() >= static_cast(std::numeric_limits::max())) { + throw std::runtime_error("marching cubes produced too many vertices for int32 faces"); + } + const int index = static_cast(values_.size()); + vertices_.push_back(static_cast(x)); + vertices_.push_back(static_cast(y)); + vertices_.push_back(static_cast(z)); + normals_.insert(normals_.end(), {0.0F, 0.0F, 0.0F}); + values_.push_back(0.0F); + return index; + } + + void add_face(const int index) { + faces_.push_back(static_cast(index)); + values_[static_cast(index)] = std::max( + values_[static_cast(index)], static_cast(vmax_) + ); + } + + void add_gradient(const int vertex, const double x, const double y, const double z) { + const std::size_t base = static_cast(vertex) * 3; + normals_[base] += static_cast(x); + normals_[base + 1] += static_cast(y); + normals_[base + 2] += static_cast(z); + } + + void add_gradient_from_corner(const int vertex, const int corner, const double strength) { + const std::size_t base = static_cast(corner) * 3; + add_gradient(vertex, gradient_[base] * strength, gradient_[base + 1] * strength, + gradient_[base + 2] * strength); + } + + std::pair *, std::size_t> face_layer_index(int edge) { + std::size_t cube = static_cast(nx_) * static_cast(y_) + + static_cast(x_); + std::size_t slot = 0; + std::vector *layer = nullptr; + if (edge < 8) { + if (edge < 4) { + layer = &face_layer1_; + } else { + edge -= 4; + layer = &face_layer2_; + } + if (edge == 1) { + cube += static_cast(step_); + slot = 1; + } else if (edge == 2) { + cube += static_cast(nx_) * static_cast(step_); + } else if (edge == 3) { + slot = 1; + } + } else { + layer = &face_layer1_; + slot = 2; + if (edge == 9) cube += static_cast(step_); + else if (edge == 10) { + cube += (static_cast(nx_) + 1) + * static_cast(step_); + } else if (edge == 11) { + cube += static_cast(nx_) * static_cast(step_); + } + } + return {layer, 4 * cube + slot}; + } + + void add_face_from_edge(const int edge) { + if (edge == 12) { + if (!center_calculated_) calculate_center_vertex(); + if (center_vertex_ < 0) { + center_vertex_ = add_vertex(center_x_, center_y_, center_z_); + } + add_face(center_vertex_); + add_gradient(center_vertex_, center_gx_, center_gy_, center_gz_); + return; + } + + const int dx1 = kEdgeRelativeX[static_cast(edge)][0]; + const int dx2 = kEdgeRelativeX[static_cast(edge)][1]; + const int dy1 = kEdgeRelativeY[static_cast(edge)][0]; + const int dy2 = kEdgeRelativeY[static_cast(edge)][1]; + const int dz1 = kEdgeRelativeZ[static_cast(edge)][0]; + const int dz2 = kEdgeRelativeZ[static_cast(edge)][1]; + const int corner1 = dz1 * 4 + dy1 * 2 + dx1; + const int corner2 = dz2 * 4 + dy2 * 2 + dx2; + const double strength1 = 1.0 / (kEpsilon + std::abs(values_by_corner_[corner1])); + const double strength2 = 1.0 / (kEpsilon + std::abs(values_by_corner_[corner2])); + const auto [layer, layer_index] = face_layer_index(edge); + int vertex = (*layer)[layer_index]; + if (vertex < 0) { + const double weight = strength1 + strength2; + vertex = add_vertex( + x_ + step_ * (dx1 * strength1 + dx2 * strength2) / weight, + y_ + step_ * (dy1 * strength1 + dy2 * strength2) / weight, + z_ + step_ * (dz1 * strength1 + dz2 * strength2) / weight + ); + (*layer)[layer_index] = vertex; + } + add_face(vertex); + add_gradient_from_corner(vertex, corner1, strength1); + add_gradient_from_corner(vertex, corner2, strength2); + } + + void prepare_for_triangles() { + values_by_corner_ = {v0_, v1_, v3_, v2_, v4_, v5_, v7_, v6_}; + double minimum = 0.0; + double maximum = 0.0; + for (const double value : values_by_corner_) { + minimum = std::min(minimum, value); + maximum = std::max(maximum, value); + } + vmax_ = maximum - minimum; + auto set_gradient = [this](const int corner, const double x, const double y, const double z) { + const std::size_t base = static_cast(corner) * 3; + gradient_[base] = x; + gradient_[base + 1] = y; + gradient_[base + 2] = z; + }; + set_gradient(0, v0_ - v1_, v0_ - v3_, v0_ - v4_); + set_gradient(1, v0_ - v1_, v1_ - v2_, v1_ - v5_); + set_gradient(2, v3_ - v2_, v1_ - v2_, v2_ - v6_); + set_gradient(3, v3_ - v2_, v0_ - v3_, v3_ - v7_); + set_gradient(4, v4_ - v5_, v4_ - v7_, v0_ - v4_); + set_gradient(5, v4_ - v5_, v5_ - v6_, v1_ - v5_); + set_gradient(6, v7_ - v6_, v5_ - v6_, v2_ - v6_); + set_gradient(7, v7_ - v6_, v4_ - v7_, v3_ - v7_); + } + + void calculate_center_vertex() { + const std::array strengths{{ + 1.0 / (kEpsilon + std::abs(v0_)), 1.0 / (kEpsilon + std::abs(v1_)), + 1.0 / (kEpsilon + std::abs(v2_)), 1.0 / (kEpsilon + std::abs(v3_)), + 1.0 / (kEpsilon + std::abs(v4_)), 1.0 / (kEpsilon + std::abs(v5_)), + 1.0 / (kEpsilon + std::abs(v6_)), 1.0 / (kEpsilon + std::abs(v7_)), + }}; + const std::array xs{{0, 1, 1, 0, 0, 1, 1, 0}}; + const std::array ys{{0, 0, 1, 1, 0, 0, 1, 1}}; + const std::array zs{{0, 0, 0, 0, 1, 1, 1, 1}}; + double x = 0.0; + double y = 0.0; + double z = 0.0; + double sum = 0.0; + center_gx_ = 0.0; + center_gy_ = 0.0; + center_gz_ = 0.0; + for (int corner = 0; corner < 8; ++corner) { + const double strength = strengths[static_cast(corner)]; + x += xs[static_cast(corner)] * strength; + y += ys[static_cast(corner)] * strength; + z += zs[static_cast(corner)] * strength; + sum += strength; + const std::size_t base = static_cast(corner) * 3; + center_gx_ += strength * gradient_[base]; + center_gy_ += strength * gradient_[base + 1]; + center_gz_ += strength * gradient_[base + 2]; + } + center_x_ = x_ + step_ * x / sum; + center_y_ = y_ + step_ * y / sum; + center_z_ = z_ + step_ * z / sum; + center_calculated_ = true; + } + + int nx_ = 0; + int ny_ = 0; + int x_ = 0; + int y_ = 0; + int z_ = 0; + int step_ = 1; + int index_ = 0; + double v0_ = 0.0; + double v1_ = 0.0; + double v2_ = 0.0; + double v3_ = 0.0; + double v4_ = 0.0; + double v5_ = 0.0; + double v6_ = 0.0; + double v7_ = 0.0; + double vmax_ = 0.0; + std::array values_by_corner_{}; + std::array gradient_{}; + bool center_calculated_ = false; + int center_vertex_ = -1; + double center_x_ = 0.0; + double center_y_ = 0.0; + double center_z_ = 0.0; + double center_gx_ = 0.0; + double center_gy_ = 0.0; + double center_gz_ = 0.0; + std::vector face_layer1_; + std::vector face_layer2_; + std::vector vertices_; + std::vector faces_; + std::vector normals_; + std::vector values_; +}; + +inline bool test_face(const Cell &cell, const int face) { + const int absolute_face = std::abs(face); + double a = 0.0; + double b = 0.0; + double c = 0.0; + double d = 0.0; + if (absolute_face == 1) { a = cell.v0(); b = cell.v4(); c = cell.v5(); d = cell.v1(); } + else if (absolute_face == 2) { a = cell.v1(); b = cell.v5(); c = cell.v6(); d = cell.v2(); } + else if (absolute_face == 3) { a = cell.v2(); b = cell.v6(); c = cell.v7(); d = cell.v3(); } + else if (absolute_face == 4) { a = cell.v3(); b = cell.v7(); c = cell.v4(); d = cell.v0(); } + else if (absolute_face == 5) { a = cell.v0(); b = cell.v3(); c = cell.v2(); d = cell.v1(); } + else if (absolute_face == 6) { a = cell.v4(); b = cell.v7(); c = cell.v6(); d = cell.v5(); } + const double determinant = a * c - b * d; + if (determinant > -kEpsilon && determinant < kEpsilon) return face >= 0; + return face * a * determinant >= 0.0; +} + +inline bool test_internal(const Cell &cell, const Luts &luts, const int case_, const int config, + const int subconfig, const int sign); +inline void select_mc33_tiling(const Luts &luts, Cell &cell, int case_, int config); + +inline void remove_degenerate_faces(MarchingCubesResult &result) { + const std::size_t n_vertices = result.values.size(); + std::vector map(n_vertices); + for (std::size_t i = 0; i < n_vertices; ++i) map[i] = static_cast(i); + std::vector keep_face(result.faces.size() / 3, true); + const auto equal_vertex = [&result](const std::int32_t a, const std::int32_t b) { + const std::size_t first = static_cast(a) * 3; + const std::size_t second = static_cast(b) * 3; + return result.vertices[first] == result.vertices[second] + && result.vertices[first + 1] == result.vertices[second + 1] + && result.vertices[first + 2] == result.vertices[second + 2]; + }; + for (std::size_t face = 0; face < keep_face.size(); ++face) { + const std::int32_t a = result.faces[face * 3]; + const std::int32_t b = result.faces[face * 3 + 1]; + const std::int32_t c = result.faces[face * 3 + 2]; + if (equal_vertex(a, b)) { map[static_cast(a)] = map[static_cast(b)] = std::min(map[static_cast(a)], map[static_cast(b)]); keep_face[face] = false; } + if (equal_vertex(a, c)) { map[static_cast(a)] = map[static_cast(c)] = std::min(map[static_cast(a)], map[static_cast(c)]); keep_face[face] = false; } + if (equal_vertex(b, c)) { map[static_cast(b)] = map[static_cast(c)] = std::min(map[static_cast(b)], map[static_cast(c)]); keep_face[face] = false; } + } + std::vector compact(n_vertices, -1); + MarchingCubesResult output; + for (std::size_t vertex = 0; vertex < n_vertices; ++vertex) { + if (map[vertex] == static_cast(vertex)) { + compact[vertex] = static_cast(output.values.size()); + output.vertices.insert(output.vertices.end(), result.vertices.begin() + static_cast(vertex * 3), result.vertices.begin() + static_cast(vertex * 3 + 3)); + output.normals.insert(output.normals.end(), result.normals.begin() + static_cast(vertex * 3), result.normals.begin() + static_cast(vertex * 3 + 3)); + output.values.push_back(result.values[vertex]); + } + } + for (std::size_t face = 0; face < keep_face.size(); ++face) { + if (!keep_face[face]) continue; + for (int corner = 0; corner < 3; ++corner) { + const auto old = static_cast(result.faces[face * 3 + static_cast(corner)]); + output.faces.push_back(compact[static_cast(map[old])]); + } + } + result = std::move(output); +} + +inline bool test_internal( + const Cell &cell, + const Luts &luts, + const int case_, + const int config, + const int subconfig, + const int sign +) { + double at = 0.0; + double bt = 0.0; + double ct = 0.0; + double dt = 0.0; + if (case_ == 4 || case_ == 10) { + const double a = (cell.v4() - cell.v0()) * (cell.v6() - cell.v2()) + - (cell.v7() - cell.v3()) * (cell.v5() - cell.v1()); + const double b = cell.v2() * (cell.v4() - cell.v0()) + + cell.v0() * (cell.v6() - cell.v2()) - cell.v1() * (cell.v7() - cell.v3()) + - cell.v3() * (cell.v5() - cell.v1()); + const double t = -b / (2.0 * a + kEpsilon); + if (t < 0.0 || t > 1.0) return sign > 0; + at = cell.v0() + (cell.v4() - cell.v0()) * t; + bt = cell.v3() + (cell.v7() - cell.v3()) * t; + ct = cell.v2() + (cell.v6() - cell.v2()) * t; + dt = cell.v1() + (cell.v5() - cell.v1()) * t; + } else { + int edge = -1; + if (case_ == 6) edge = luts.test6.get2(config, 2); + else if (case_ == 7) edge = luts.test7.get2(config, 4); + else if (case_ == 12) edge = luts.test12.get2(config, 3); + else if (case_ == 13) edge = luts.tiling13_5_1.get3(config, subconfig, 0); + else return false; + double t = 0.0; + switch (edge) { + case 0: + t = cell.v0() / (cell.v0() - cell.v1() + kEpsilon); + bt = cell.v3() + (cell.v2() - cell.v3()) * t; + ct = cell.v7() + (cell.v6() - cell.v7()) * t; + dt = cell.v4() + (cell.v5() - cell.v4()) * t; + break; + case 1: + t = cell.v1() / (cell.v1() - cell.v2() + kEpsilon); + bt = cell.v0() + (cell.v3() - cell.v0()) * t; + ct = cell.v4() + (cell.v7() - cell.v4()) * t; + dt = cell.v5() + (cell.v6() - cell.v5()) * t; + break; + case 2: + t = cell.v2() / (cell.v2() - cell.v3() + kEpsilon); + bt = cell.v1() + (cell.v0() - cell.v1()) * t; + ct = cell.v5() + (cell.v4() - cell.v5()) * t; + dt = cell.v6() + (cell.v7() - cell.v6()) * t; + break; + case 3: + t = cell.v3() / (cell.v3() - cell.v0() + kEpsilon); + bt = cell.v2() + (cell.v1() - cell.v2()) * t; + ct = cell.v6() + (cell.v5() - cell.v6()) * t; + dt = cell.v7() + (cell.v4() - cell.v7()) * t; + break; + case 4: + t = cell.v4() / (cell.v4() - cell.v5() + kEpsilon); + bt = cell.v7() + (cell.v6() - cell.v7()) * t; + ct = cell.v3() + (cell.v2() - cell.v3()) * t; + dt = cell.v0() + (cell.v1() - cell.v0()) * t; + break; + case 5: + t = cell.v5() / (cell.v5() - cell.v6() + kEpsilon); + bt = cell.v4() + (cell.v7() - cell.v4()) * t; + ct = cell.v0() + (cell.v3() - cell.v0()) * t; + dt = cell.v1() + (cell.v2() - cell.v1()) * t; + break; + case 6: + t = cell.v6() / (cell.v6() - cell.v7() + kEpsilon); + bt = cell.v5() + (cell.v4() - cell.v5()) * t; + ct = cell.v1() + (cell.v0() - cell.v1()) * t; + dt = cell.v2() + (cell.v3() - cell.v2()) * t; + break; + case 7: + t = cell.v7() / (cell.v7() - cell.v4() + kEpsilon); + bt = cell.v6() + (cell.v5() - cell.v6()) * t; + ct = cell.v2() + (cell.v1() - cell.v2()) * t; + dt = cell.v3() + (cell.v0() - cell.v3()) * t; + break; + case 8: + t = cell.v0() / (cell.v0() - cell.v4() + kEpsilon); + bt = cell.v3() + (cell.v7() - cell.v3()) * t; + ct = cell.v2() + (cell.v6() - cell.v2()) * t; + dt = cell.v1() + (cell.v5() - cell.v1()) * t; + break; + case 9: + t = cell.v1() / (cell.v1() - cell.v5() + kEpsilon); + bt = cell.v0() + (cell.v4() - cell.v0()) * t; + ct = cell.v3() + (cell.v7() - cell.v3()) * t; + dt = cell.v2() + (cell.v6() - cell.v2()) * t; + break; + case 10: + t = cell.v2() / (cell.v2() - cell.v6() + kEpsilon); + bt = cell.v1() + (cell.v5() - cell.v1()) * t; + ct = cell.v0() + (cell.v4() - cell.v0()) * t; + dt = cell.v3() + (cell.v7() - cell.v3()) * t; + break; + case 11: + t = cell.v3() / (cell.v3() - cell.v7() + kEpsilon); + bt = cell.v2() + (cell.v6() - cell.v2()) * t; + ct = cell.v1() + (cell.v5() - cell.v1()) * t; + dt = cell.v0() + (cell.v4() - cell.v0()) * t; + break; + default: + return false; + } + } + int test = 0; + if (at >= 0.0) test += 1; + if (bt >= 0.0) test += 2; + if (ct >= 0.0) test += 4; + if (dt >= 0.0) test += 8; + switch (test) { + case 0: case 1: case 2: case 3: case 4: case 6: case 8: case 9: + return sign > 0; + case 5: + return at * ct - bt * dt < kEpsilon ? sign > 0 : false; + case 7: case 11: case 13: case 14: case 15: + return sign < 0; + case 12: + return sign > 0; + case 10: + return at * ct - bt * dt >= kEpsilon ? sign > 0 : false; + default: + return sign < 0; + } +} + +inline void select_mc33_tiling(const Luts &luts, Cell &cell, const int case_, const int config) { + int subconfig = 0; + switch (case_) { + case 1: cell.add_triangles(luts.tiling1, config, 1); break; + case 2: cell.add_triangles(luts.tiling2, config, 2); break; + case 3: { + const bool split = test_face(cell, luts.test3.get1(config)); + cell.add_triangles(split ? luts.tiling3_2 : luts.tiling3_1, config, split ? 4 : 2); + break; + } + case 4: { + const bool connected = test_internal(cell, luts, case_, config, subconfig, luts.test4.get1(config)); + cell.add_triangles(connected ? luts.tiling4_1 : luts.tiling4_2, config, connected ? 2 : 6); + break; + } + case 5: cell.add_triangles(luts.tiling5, config, 3); break; + case 6: + if (test_face(cell, luts.test6.get2(config, 0))) cell.add_triangles(luts.tiling6_2, config, 5); + else if (test_internal(cell, luts, case_, config, subconfig, luts.test6.get2(config, 1))) cell.add_triangles(luts.tiling6_1_1, config, 3); + else cell.add_triangles(luts.tiling6_1_2, config, 9); + break; + case 7: + if (test_face(cell, luts.test7.get2(config, 0))) subconfig += 1; + if (test_face(cell, luts.test7.get2(config, 1))) subconfig += 2; + if (test_face(cell, luts.test7.get2(config, 2))) subconfig += 4; + if (subconfig == 0) cell.add_triangles(luts.tiling7_1, config, 3); + else if (subconfig == 1 || subconfig == 2 || subconfig == 4) cell.add_triangles2(luts.tiling7_2, config, subconfig == 1 ? 0 : subconfig == 2 ? 1 : 2, 5); + else if (subconfig == 3 || subconfig == 5 || subconfig == 6) cell.add_triangles2(luts.tiling7_3, config, subconfig == 3 ? 0 : subconfig == 5 ? 1 : 2, 9); + else if (test_internal(cell, luts, case_, config, subconfig, luts.test7.get2(config, 3))) cell.add_triangles(luts.tiling7_4_2, config, 9); + else cell.add_triangles(luts.tiling7_4_1, config, 5); + break; + case 8: cell.add_triangles(luts.tiling8, config, 2); break; + case 9: cell.add_triangles(luts.tiling9, config, 4); break; + case 10: + if (test_face(cell, luts.test10.get2(config, 0))) { + if (test_face(cell, luts.test10.get2(config, 1))) cell.add_triangles(luts.tiling10_1_1_alt, config, 4); + else cell.add_triangles(luts.tiling10_2, config, 8); + } else if (test_face(cell, luts.test10.get2(config, 1))) cell.add_triangles(luts.tiling10_2_alt, config, 8); + else if (test_internal(cell, luts, case_, config, subconfig, luts.test10.get2(config, 2))) cell.add_triangles(luts.tiling10_1_1, config, 4); + else cell.add_triangles(luts.tiling10_1_2, config, 8); + break; + case 11: cell.add_triangles(luts.tiling11, config, 4); break; + case 12: + if (test_face(cell, luts.test12.get2(config, 0))) { + if (test_face(cell, luts.test12.get2(config, 1))) cell.add_triangles(luts.tiling12_1_1_alt, config, 4); + else cell.add_triangles(luts.tiling12_2, config, 8); + } else if (test_face(cell, luts.test12.get2(config, 1))) cell.add_triangles(luts.tiling12_2_alt, config, 8); + else if (test_internal(cell, luts, case_, config, subconfig, luts.test12.get2(config, 2))) cell.add_triangles(luts.tiling12_1_1, config, 4); + else cell.add_triangles(luts.tiling12_1_2, config, 8); + break; + case 13: + for (int face = 0; face < 6; ++face) { + if (test_face(cell, luts.test13.get2(config, face))) subconfig |= 1 << face; + } + subconfig = luts.subconfig13.get1(subconfig); + if (subconfig == 0) cell.add_triangles(luts.tiling13_1, config, 4); + else if (subconfig >= 1 && subconfig <= 6) cell.add_triangles2(luts.tiling13_2, config, subconfig - 1, 6); + else if (subconfig >= 7 && subconfig <= 18) cell.add_triangles2(luts.tiling13_3, config, subconfig - 7, 10); + else if (subconfig >= 19 && subconfig <= 22) cell.add_triangles2(luts.tiling13_4, config, subconfig - 19, 12); + else if (subconfig >= 23 && subconfig <= 26) { + const int local = subconfig - 23; + if (test_internal(cell, luts, case_, config, local, luts.test13.get2(config, 6))) cell.add_triangles2(luts.tiling13_5_1, config, local, 6); + else cell.add_triangles2(luts.tiling13_5_2, config, local, 10); + } else if (subconfig >= 27 && subconfig <= 38) cell.add_triangles2(luts.tiling13_3_alt, config, subconfig - 27, 10); + else if (subconfig >= 39 && subconfig <= 44) cell.add_triangles2(luts.tiling13_2_alt, config, subconfig - 39, 6); + else if (subconfig == 45) cell.add_triangles(luts.tiling13_1_alt, config, 4); + break; + case 14: cell.add_triangles(luts.tiling14, config, 4); break; + default: break; + } +} + +template +inline void traverse_cells( + const ConstArrayView &volume, + const double level, + const int step_size, + const MarchingCubesMethod method, + const ConstArrayView *mask, + const Luts &tables, + Cell &cell +) { + const int nz = static_cast(volume.shape[0]); + const int ny = static_cast(volume.shape[1]); + const int nx = static_cast(volume.shape[2]); + const auto at = [nx, ny, &volume](const int z, const int y, const int x) { + return volume.data[(static_cast(z) * ny + y) * nx + x]; + }; + const int max_x = nx - 2 * step_size; + const int max_y = ny - 2 * step_size; + const int max_z = nz - 2 * step_size; + for (int z = -step_size; z < max_z;) { + z += step_size; + cell.new_z_value(); + const int z_next = z + step_size; + for (int y = -step_size; y < max_y;) { + y += step_size; + const int y_next = y + step_size; + for (int x = -step_size; x < max_x;) { + x += step_size; + const int x_next = x + step_size; + if constexpr (UseMask) { + const auto mask_index = + (static_cast(z_next) * ny + y_next) * nx + x_next; + if (mask->data[mask_index] == 0) continue; + } + cell.set_cube( + level, x, y, z, step_size, + at(z, y, x), at(z, y, x_next), + at(z, y_next, x_next), at(z, y_next, x), + at(z_next, y, x), at(z_next, y, x_next), + at(z_next, y_next, x_next), at(z_next, y_next, x) + ); + if (method == MarchingCubesMethod::Lorensen) { + int triangles = 0; + while (triangles < 5 && tables.cases_classic.get2(cell.index(), 3 * triangles) != -1) ++triangles; + if (triangles != 0) cell.add_triangles(tables.cases_classic, cell.index(), triangles); + } else { + const int case_ = tables.cases.get2(cell.index(), 0); + if (case_ > 0) select_mc33_tiling(tables, cell, case_, tables.cases.get2(cell.index(), 1)); + } + } + } + } +} + +} // namespace detail::marching_cubes + +// Extract an isosurface from a 3-D float32 image. The core deliberately does +// not know about Python layout conventions beyond the C-order ArrayView; the +// binding handles argument validation, output axis order, and spacing. +inline MarchingCubesResult marching_cubes( + const ConstArrayView &volume, + const double level, + const int step_size, + const MarchingCubesMethod method, + const ConstArrayView *mask, + const bool allow_degenerate +) { + BIOIMAGE_PROFILE_INIT(profiler); + if (volume.shape.size() != 3) { + throw std::invalid_argument("volume must have ndim=3"); + } + for (const auto axis : volume.shape) { + if (axis < 2) throw std::invalid_argument("volume dimensions must all be at least 2"); + if (axis > std::numeric_limits::max()) { + throw std::invalid_argument("volume dimensions exceed marching cubes index range"); + } + } + if (step_size < 1) throw std::invalid_argument("step_size must be at least 1"); + if (mask != nullptr && mask->shape != volume.shape) { + throw std::invalid_argument("mask must have the same shape as volume"); + } + const int ny = static_cast(volume.shape[1]); + const int nx = static_cast(volume.shape[2]); + const detail::marching_cubes::Luts *tables_ptr = nullptr; + { + BIOIMAGE_PROFILE_SCOPE(profiler, "lookup_tables"); + tables_ptr = &detail::marching_cubes::luts(); + } + const auto &tables = *tables_ptr; + std::optional cell; + cell.emplace(nx, ny); + { + BIOIMAGE_PROFILE_SCOPE(profiler, "cell_traversal"); + if (mask == nullptr) { + detail::marching_cubes::traverse_cells( + volume, level, step_size, method, nullptr, tables, *cell + ); + } else { + detail::marching_cubes::traverse_cells( + volume, level, step_size, method, mask, tables, *cell + ); + } + } + MarchingCubesResult result; + { + BIOIMAGE_PROFILE_SCOPE(profiler, "finalize_normals"); + result = cell->take_result(); + } + { + BIOIMAGE_PROFILE_SCOPE(profiler, "cell_cleanup"); + cell.reset(); + } + if (!allow_degenerate) { + BIOIMAGE_PROFILE_SCOPE(profiler, "remove_degenerate_faces"); + detail::marching_cubes::remove_degenerate_faces(result); + } + BIOIMAGE_PROFILE_REPORT(profiler); + return result; +} + +} // namespace bioimage_cpp::mesh diff --git a/src/bindings/mesh.cxx b/src/bindings/mesh.cxx new file mode 100644 index 0000000..2e9000c --- /dev/null +++ b/src/bindings/mesh.cxx @@ -0,0 +1,182 @@ +#include "mesh.hxx" + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/profile.hxx" +#include "bioimage_cpp/mesh/marching_cubes.hxx" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace nb = nanobind; + +namespace bioimage_cpp::bindings { +namespace { + +using FloatInput = nb::ndarray; +using UInt8Input = nb::ndarray; +using FloatOutput = nb::ndarray; +using Int32Output = nb::ndarray; + +template +Output output_array(std::vector &&values, const std::vector &shape) { + if (values.empty()) { + // A zero-size allocation keeps the data pointer valid for nanobind's + // ndarray constructor without paying an elementwise copy. + auto allocation = std::make_unique(0); + auto *data = allocation.get(); + nb::capsule owner(data, [](void *p) noexcept { delete[] static_cast(p); }); + allocation.release(); + return Output(data, shape.size(), shape.data(), owner); + } + + // NumPy owns the heap-allocated vector through the capsule. The vector is + // never resized after this point, so data() remains valid for the full + // lifetime of the returned ndarray. + auto allocation = std::make_unique>(std::move(values)); + auto *data = allocation->data(); + nb::capsule owner(allocation.get(), [](void *p) noexcept { + delete static_cast *>(p); + }); + allocation.release(); + return Output(data, shape.size(), shape.data(), owner); +} + +std::vector shape_of(const FloatInput &array) { + std::vector shape(array.ndim()); + for (std::size_t axis = 0; axis < array.ndim(); ++axis) { + shape[axis] = static_cast(array.shape(axis)); + } + return shape; +} + +nb::tuple marching_cubes_float32( + FloatInput volume, + const double level, + const int step_size, + const bool classic, + const bool descent, + std::optional mask, + const bool allow_degenerate +) { + BIOIMAGE_PROFILE_INIT(profiler); + if (volume.ndim() != 3) { + throw std::invalid_argument( + "volume must have ndim=3, got ndim=" + std::to_string(volume.ndim()) + ); + } + for (std::size_t axis = 0; axis < 3; ++axis) { + if (volume.shape(axis) < 2) { + throw std::invalid_argument("volume dimensions must all be at least 2"); + } + } + if (step_size < 1) { + throw std::invalid_argument("step_size must be at least 1"); + } + if (mask.has_value()) { + if (mask->ndim() != 3) { + throw std::invalid_argument( + "mask must have ndim=3, got ndim=" + std::to_string(mask->ndim()) + ); + } + for (std::size_t axis = 0; axis < 3; ++axis) { + if (mask->shape(axis) != volume.shape(axis)) { + throw std::invalid_argument("mask must have the same shape as volume"); + } + } + } + + const auto shape = shape_of(volume); + ConstArrayView volume_view{volume.data(), shape, {}}; + ConstArrayView mask_view; + const ConstArrayView *mask_ptr = nullptr; + if (mask.has_value()) { + std::vector mask_shape(shape.begin(), shape.end()); + mask_view = ConstArrayView{mask->data(), std::move(mask_shape), {}}; + mask_ptr = &mask_view; + } + + mesh::MarchingCubesResult result; + { + BIOIMAGE_PROFILE_SCOPE(profiler, "core_call"); + nb::gil_scoped_release release; + result = mesh::marching_cubes( + volume_view, + level, + step_size, + classic ? mesh::MarchingCubesMethod::Lorensen : mesh::MarchingCubesMethod::Lewiner, + mask_ptr, + allow_degenerate + ); + } + if (result.values.empty()) { + throw std::runtime_error("No surface found at the given iso value."); + } + + { + BIOIMAGE_PROFILE_SCOPE(profiler, "orient_output"); + // The core mirrors the reference kernel's x/y/z convention. Public arrays + // use NumPy's z/y/x axis order, matching skimage.measure.marching_cubes. + for (std::size_t vertex = 0; vertex < result.values.size(); ++vertex) { + const std::size_t base = vertex * 3; + std::swap(result.vertices[base], result.vertices[base + 2]); + std::swap(result.normals[base], result.normals[base + 2]); + } + if (descent) { + for (std::size_t face = 0; face < result.faces.size(); face += 3) { + std::swap(result.faces[face], result.faces[face + 2]); + } + } + } + + const std::size_t n_vertices = result.values.size(); + const std::size_t n_faces = result.faces.size() / 3; + FloatOutput vertices; + Int32Output faces; + FloatOutput normals; + FloatOutput values; + { + BIOIMAGE_PROFILE_SCOPE(profiler, "numpy_handoff"); + vertices = output_array( + std::move(result.vertices), {n_vertices, 3} + ); + faces = output_array( + std::move(result.faces), {n_faces, 3} + ); + normals = output_array( + std::move(result.normals), {n_vertices, 3} + ); + values = output_array( + std::move(result.values), {n_vertices} + ); + } + BIOIMAGE_PROFILE_REPORT(profiler); + return nb::make_tuple(std::move(vertices), std::move(faces), std::move(normals), std::move(values)); +} + +} // namespace + +void bind_mesh(nb::module_ &m) { + m.def( + "_marching_cubes_float32", + &marching_cubes_float32, + nb::arg("volume"), + nb::arg("level"), + nb::arg("step_size"), + nb::arg("classic"), + nb::arg("descent"), + nb::arg("mask") = nb::none(), + nb::arg("allow_degenerate") = true, + "Marching Cubes 33/Lorensen extraction from a contiguous float32 3-D volume." + ); +} + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/mesh.hxx b/src/bindings/mesh.hxx new file mode 100644 index 0000000..2c2c71f --- /dev/null +++ b/src/bindings/mesh.hxx @@ -0,0 +1,9 @@ +#pragma once + +#include + +namespace bioimage_cpp::bindings { + +void bind_mesh(nanobind::module_ &m); + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/module.cxx b/src/bindings/module.cxx index 909524f..72b6e2d 100644 --- a/src/bindings/module.cxx +++ b/src/bindings/module.cxx @@ -6,6 +6,7 @@ #include "graph.hxx" #include "ground_truth.hxx" #include "label_multiset.hxx" +#include "mesh.hxx" #include "segmentation.hxx" #include "transformation.hxx" #include "util.hxx" @@ -25,6 +26,7 @@ NB_MODULE(_core, m) { bioimage_cpp::bindings::bind_graph(m); bioimage_cpp::bindings::bind_ground_truth(m); bioimage_cpp::bindings::bind_label_multiset(m); + bioimage_cpp::bindings::bind_mesh(m); bioimage_cpp::bindings::bind_segmentation(m); bioimage_cpp::bindings::bind_transformation(m); bioimage_cpp::bindings::bind_util(m); diff --git a/src/bioimage_cpp/__init__.py b/src/bioimage_cpp/__init__.py index 631d732..e5ce86f 100644 --- a/src/bioimage_cpp/__init__.py +++ b/src/bioimage_cpp/__init__.py @@ -56,6 +56,7 @@ - `distance`: distance transform functionality. - `filters`: efficient implementation of convolutional image filters. - `graph`: graph creation and graph (partitioning) algorithms. +- `mesh`: triangle-mesh extraction from 3D volumes and segmentation masks. - `segmentation`: image segmentation functionality. - `transformation`: affine transformations. - `utils`: misc utility functionality. @@ -97,6 +98,7 @@ from . import flow from . import graph from . import label_multiset +from . import mesh from . import segmentation from . import transformation from . import utils @@ -109,6 +111,7 @@ "flow", "graph", "label_multiset", + "mesh", "segmentation", "transformation", "utils", diff --git a/src/bioimage_cpp/mesh/__init__.py b/src/bioimage_cpp/mesh/__init__.py new file mode 100644 index 0000000..d0c0ffd --- /dev/null +++ b/src/bioimage_cpp/mesh/__init__.py @@ -0,0 +1,169 @@ +"""Triangle-mesh extraction from 3-D scalar volumes.""" + +from __future__ import annotations + +from collections.abc import Sequence +import operator + +import numpy as np + +from .. import _core + + +def _as_spacing(spacing: Sequence[float]) -> np.ndarray: + try: + values = np.asarray(spacing, dtype=np.float64) + except (TypeError, ValueError) as error: + raise ValueError("spacing must consist of three floats") from error + if values.shape != (3,): + raise ValueError("spacing must consist of three floats") + if not np.all(np.isfinite(values)) or np.any(values <= 0.0): + raise ValueError("spacing entries must be positive and finite") + return values + + +def marching_cubes( + volume: np.ndarray, + level: float | None = None, + *, + spacing: Sequence[float] = (1.0, 1.0, 1.0), + gradient_direction: str = "descent", + step_size: int = 1, + allow_degenerate: bool = True, + method: str = "lewiner", + mask: np.ndarray | None = None, + pad: bool = False, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Extract an isosurface from a three-dimensional scalar volume. + + This follows :func:`skimage.measure.marching_cubes` for its volume, + level, spacing, winding, stride, degenerate-face, method, and mask + semantics. ``method="lewiner"`` is the topology-resolving Marching Cubes + 33 implementation and is the default; ``method="lorensen"`` uses the + original classic 256-case lookup table. + + Parameters + ---------- + volume: + Three-dimensional numeric array. It is converted to C-contiguous + ``float32`` before entering the C++ kernel. + level: + Iso-value to extract. ``None`` uses the midpoint of the data range. + spacing: + Three physical spacings in NumPy axis order ``(z, y, x)``. + gradient_direction: + ``"descent"`` (the default) treats objects as values larger than the + exterior; ``"ascent"`` reverses triangle winding. + step_size: + Sampling stride in voxels. Values above one produce a coarser mesh. + allow_degenerate: + If false, remove faces with repeated vertex coordinates. + method: + ``"lewiner"`` or ``"lorensen"``. + mask: + Optional boolean region of the same shape. Cubes are emitted only in + its true region; this can intentionally produce open surfaces. + pad: + If true, add a one-voxel zero-valued halo before extraction and shift + returned coordinates back to the original volume origin. This closes + foreground objects that touch the input boundary. The padded halo is + enabled in a supplied mask so that its boundary cells can be emitted. + + Returns + ------- + vertices, faces, normals, values: + Vertices and normals have shape ``(V, 3)`` in NumPy ``(z, y, x)`` + coordinate order; faces has shape ``(F, 3)`` and dtype ``int32``; + values has shape ``(V,)``. At unit spacing vertices are ``float32``; + non-unit spacing follows skimage and produces ``float64`` vertices. + Normals are normalized gradients accumulated from incident cells and + ``values`` stores the largest local data range seen at each vertex, as in + scikit-image. ``gradient_direction`` changes face winding only; spacing + scales vertices but does not transform normals. + + Raises + ------ + ValueError + If shapes or options are invalid, or the requested level is outside + the input data range. + RuntimeError + If no surface intersects the chosen level. + """ + volume_array = np.asarray(volume) + if volume_array.ndim != 3: + raise ValueError(f"Input volume should be a 3D numpy array, got ndim={volume_array.ndim}.") + if any(size < 2 for size in volume_array.shape): + raise ValueError("Input array must be at least 2x2x2.") + if not np.issubdtype(volume_array.dtype, np.number) and volume_array.dtype != np.bool_: + raise TypeError(f"volume must have a numeric dtype, got dtype={volume_array.dtype}") + if np.issubdtype(volume_array.dtype, np.complexfloating): + raise TypeError(f"volume must have a real numeric dtype, got dtype={volume_array.dtype}") + volume_float = np.ascontiguousarray(volume_array, dtype=np.float32) + + volume_min = float(volume_float.min()) + volume_max = float(volume_float.max()) + if level is None: + level_float = 0.5 * (volume_min + volume_max) + else: + level_float = float(level) + if not np.isfinite(level_float): + raise ValueError("level must be finite") + if level_float < volume_min or level_float > volume_max: + raise ValueError("Surface level must be within volume data range.") + + spacing_array = _as_spacing(spacing) + try: + step = operator.index(step_size) + except TypeError as error: + raise TypeError("step_size must be an integer") from error + if step < 1: + raise ValueError("step_size must be at least one.") + if method == "lewiner": + classic = False + elif method == "lorensen": + classic = True + else: + raise ValueError("method should be either 'lewiner' or 'lorensen'") + if gradient_direction == "descent": + descent = True + elif gradient_direction == "ascent": + descent = False + else: + raise ValueError( + "Incorrect input gradient_direction, see marching_cubes documentation." + ) + + mask_array: np.ndarray | None = None + if mask is not None: + mask_input = np.asarray(mask) + if mask_input.shape != volume_array.shape: + raise ValueError("volume and mask must have the same shape.") + mask_array = np.ascontiguousarray(mask_input != 0, dtype=np.uint8) + + if bool(pad): + volume_float = np.pad(volume_float, 1, mode="constant", constant_values=0) + if mask_array is not None: + mask_array = np.pad(mask_array, 1, mode="constant", constant_values=True) + + vertices, faces, normals, values = _core._marching_cubes_float32( + volume_float, + level_float, + step, + classic, + descent, + mask_array, + bool(allow_degenerate), + ) + + if np.array_equal(spacing_array, (1.0, 1.0, 1.0)): + if bool(pad): + vertices = vertices - np.ones(3, dtype=np.float32) + return vertices, faces, normals, values + + vertices = vertices.astype(np.float64, copy=False) * spacing_array + if bool(pad): + vertices -= spacing_array + return vertices, faces, normals, values + + +__all__ = ["marching_cubes"] diff --git a/tests/mesh/test_marching_cubes.py b/tests/mesh/test_marching_cubes.py new file mode 100644 index 0000000..a143e12 --- /dev/null +++ b/tests/mesh/test_marching_cubes.py @@ -0,0 +1,281 @@ +import gc + +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def _interior_box() -> np.ndarray: + volume = np.zeros((5, 6, 7), dtype=np.uint8) + volume[1:4, 1:5, 1:6] = 1 + return volume + + +def _interior_ball(n: int = 32, radius: float = 11.0) -> np.ndarray: + center = (n - 1) / 2.0 + z, y, x = np.ogrid[:n, :n, :n] + return ( + (z - center) ** 2 + (y - center) ** 2 + (x - center) ** 2 + <= radius**2 + ).astype(np.uint8) + + +def _edge_counts(faces: np.ndarray) -> dict[tuple[int, int], int]: + counts: dict[tuple[int, int], int] = {} + for face in faces: + for first, second in ( + (face[0], face[1]), + (face[1], face[2]), + (face[2], face[0]), + ): + edge = tuple(sorted((int(first), int(second)))) + counts[edge] = counts.get(edge, 0) + 1 + return counts + + +def _zero_area_faces(vertices: np.ndarray, faces: np.ndarray) -> int: + triangles = vertices[faces] + cross = np.cross( + triangles[:, 1] - triangles[:, 0], + triangles[:, 2] - triangles[:, 0], + ) + return int(np.count_nonzero(np.linalg.norm(cross, axis=1) <= 1e-12)) + + +def _collapsed_coordinate_faces(vertices: np.ndarray, faces: np.ndarray) -> int: + triangles = vertices[faces] + first_second = np.all(triangles[:, 0] == triangles[:, 1], axis=1) + first_third = np.all(triangles[:, 0] == triangles[:, 2], axis=1) + second_third = np.all(triangles[:, 1] == triangles[:, 2], axis=1) + return int(np.count_nonzero(first_second | first_third | second_third)) + + +def test_box_mesh_has_reference_dimensions_and_dtypes(): + vertices, faces, normals, values = bic.mesh.marching_cubes(_interior_box(), 0.5) + + assert vertices.shape == (94, 3) + assert faces.shape == (184, 3) + assert normals.shape == vertices.shape + assert values.shape == (94,) + assert vertices.dtype == np.float32 + assert faces.dtype == np.int32 + assert normals.dtype == np.float32 + assert values.dtype == np.float32 + np.testing.assert_allclose(np.linalg.norm(normals, axis=1), 1.0, rtol=2e-6) + assert np.all((faces >= 0) & (faces < len(vertices))) + np.testing.assert_array_equal(vertices.min(axis=0), [0.5, 0.5, 0.5]) + np.testing.assert_array_equal(vertices.max(axis=0), [3.5, 4.5, 5.5]) + np.testing.assert_array_equal(values, np.ones_like(values)) + + +def test_output_arrays_are_contiguous_writable_and_survive_collection(): + outputs = bic.mesh.marching_cubes(_interior_box(), 0.5) + expected = tuple(output.copy() for output in outputs) + + for output in outputs: + assert output.flags.c_contiguous + assert output.flags.writeable + + # The binding transfers ownership to NumPy. Exercise the arrays after the + # C++ result has gone out of scope and after allocation churn/collection. + for _ in range(8): + np.empty((1024, 1024), dtype=np.float32) + gc.collect() + for output, want in zip(outputs, expected, strict=True): + np.testing.assert_array_equal(output, want) + + outputs[3][0] += 1.0 + assert outputs[3][0] == expected[3][0] + 1.0 + + +def test_lorensen_and_lewiner_choose_different_ambiguous_topologies(): + # The bit pattern 0b00000110 is an ambiguous one-cube configuration. + volume = np.array([(6 >> bit) & 1 for bit in range(8)], dtype=np.uint8).reshape(2, 2, 2) + lewiner = bic.mesh.marching_cubes(volume, 0.5, method="lewiner") + lorensen = bic.mesh.marching_cubes(volume, 0.5, method="lorensen") + + assert lewiner[0].shape == lorensen[0].shape == (6, 3) + assert lewiner[1].shape == (4, 3) + assert lorensen[1].shape == (2, 3) + + +def test_spacing_preserves_mesh_and_uses_float64_vertices(): + volume = _interior_box() + unit_vertices, unit_faces, unit_normals, unit_values = bic.mesh.marching_cubes(volume, 0.5) + spacing = (2.0, 0.5, 3.0) + vertices, faces, normals, values = bic.mesh.marching_cubes(volume, 0.5, spacing=spacing) + + assert vertices.dtype == np.float64 + np.testing.assert_allclose(vertices, unit_vertices * np.asarray(spacing)) + np.testing.assert_array_equal(faces, unit_faces) + np.testing.assert_array_equal(normals, unit_normals) + np.testing.assert_array_equal(values, unit_values) + + +@pytest.mark.parametrize("dtype", [np.bool_, np.uint16, np.float64]) +def test_numeric_dtype_conversion_matches_uint8_input(dtype): + volume = _interior_box() + if dtype is np.bool_: + converted = volume.astype(bool) + else: + converted = volume.astype(dtype) + expected = bic.mesh.marching_cubes(volume, 0.5) + actual = bic.mesh.marching_cubes(converted, 0.5) + for got, want in zip(actual, expected, strict=True): + np.testing.assert_array_equal(got, want) + + +def test_non_contiguous_input_is_normalised_at_the_python_boundary(): + volume = _interior_box() + repeated = np.repeat(volume, 2, axis=2) + non_contiguous = repeated[:, :, ::2] + assert not non_contiguous.flags.c_contiguous + + actual = bic.mesh.marching_cubes(non_contiguous, 0.5) + expected = bic.mesh.marching_cubes(volume, 0.5) + for got, want in zip(actual, expected, strict=True): + np.testing.assert_array_equal(got, want) + + +def test_gradient_direction_only_reverses_face_winding(): + descent = bic.mesh.marching_cubes(_interior_box(), 0.5, gradient_direction="descent") + ascent = bic.mesh.marching_cubes(_interior_box(), 0.5, gradient_direction="ascent") + + np.testing.assert_array_equal(descent[0], ascent[0]) + np.testing.assert_array_equal(descent[2], ascent[2]) + np.testing.assert_array_equal(descent[3], ascent[3]) + np.testing.assert_array_equal(descent[1], ascent[1][:, ::-1]) + + +def test_padding_closes_a_boundary_object_and_restores_coordinates(): + volume = np.zeros((4, 4, 4), dtype=np.uint8) + volume[:2, 1:3, 1:3] = 1 + open_vertices, open_faces, _, _ = bic.mesh.marching_cubes(volume, 0.5) + padded_vertices, padded_faces, _, _ = bic.mesh.marching_cubes(volume, 0.5, pad=True) + + assert open_vertices.shape == (20, 3) + assert open_faces.shape == (30, 3) + assert padded_vertices.shape == (24, 3) + assert padded_faces.shape == (44, 3) + np.testing.assert_array_equal(open_vertices.min(axis=0), [0.0, 0.5, 0.5]) + np.testing.assert_array_equal(padded_vertices.min(axis=0), [-0.5, 0.5, 0.5]) + assert set(_edge_counts(padded_faces).values()) == {2} + + +def test_mask_and_step_size_are_honoured(): + rng = np.random.default_rng(9) + volume = rng.integers(0, 2, size=(8, 9, 10), dtype=np.uint8) + mask = np.ones_like(volume, dtype=bool) + mask[:, 0] = False + vertices, faces, normals, values = bic.mesh.marching_cubes( + volume, 0.5, method="lorensen", mask=mask, step_size=2 + ) + + assert vertices.shape == normals.shape + assert vertices.shape[1] == 3 + assert len(vertices) > 0 + assert faces.shape[1] == 3 + assert len(faces) > 0 + assert values.shape == (len(vertices),) + + +@pytest.mark.parametrize("method", ["lewiner", "lorensen"]) +@pytest.mark.parametrize("step_size", [1, 2, 3]) +def test_vertex_cache_with_step_size_and_mask_is_deterministic(method, step_size): + rng = np.random.default_rng(27) + volume = rng.random((11, 12, 13), dtype=np.float32) + mask = np.ones(volume.shape, dtype=bool) + mask[1::3, 2::4, :] = False + mask[:, 1::4, 2::3] = False + + first = bic.mesh.marching_cubes( + volume, + 0.5, + method=method, + step_size=step_size, + mask=mask, + ) + second = bic.mesh.marching_cubes( + volume, + 0.5, + method=method, + step_size=step_size, + mask=mask, + ) + + for actual, expected in zip(first, second, strict=True): + np.testing.assert_array_equal(actual, expected) + assert len(np.unique(first[0], axis=0)) == len(first[0]) + assert np.all((first[1] >= 0) & (first[1] < len(first[0]))) + + +def test_lewiner_closed_sphere_is_watertight_with_euler_characteristic_two(): + vertices, faces, _, _ = bic.mesh.marching_cubes( + _interior_ball(n=34, radius=12.0), 0.5, method="lewiner" + ) + edge_counts = _edge_counts(faces) + assert set(edge_counts.values()) == {2} + assert len(vertices) - len(edge_counts) + len(faces) == 2 + + +@pytest.mark.parametrize("method", ["lewiner", "lorensen"]) +def test_output_is_deterministic(method): + volume = _interior_ball() + first = bic.mesh.marching_cubes(volume, 0.5, method=method) + second = bic.mesh.marching_cubes(volume, 0.5, method=method) + for actual, expected in zip(first, second, strict=True): + np.testing.assert_array_equal(actual, expected) + + +def test_allow_degenerate_false_removes_collapsed_faces(): + rng = np.random.default_rng(3) + volume = rng.integers(0, 3, size=(16, 16, 16), dtype=np.uint8) + kept = bic.mesh.marching_cubes(volume, 1.0, allow_degenerate=True) + removed = bic.mesh.marching_cubes(volume, 1.0, allow_degenerate=False) + + assert _collapsed_coordinate_faces(kept[0], kept[1]) > 0 + assert _collapsed_coordinate_faces(removed[0], removed[1]) == 0 + assert _zero_area_faces(removed[0], removed[1]) < _zero_area_faces(kept[0], kept[1]) + assert len(removed[1]) < len(kept[1]) + assert removed[1].dtype == np.int32 + + +def test_default_level_is_midpoint_of_original_volume_when_padding(): + volume = np.full((5, 5, 5), 2.0, dtype=np.float32) + volume[1:4, 1:4, 1:4] = 4.0 + default = bic.mesh.marching_cubes(volume, pad=True) + explicit = bic.mesh.marching_cubes(volume, 3.0, pad=True) + for actual, expected in zip(default, explicit, strict=True): + np.testing.assert_array_equal(actual, expected) + + +def test_invalid_arguments_and_missing_surface(): + with pytest.raises(ValueError, match="3D"): + bic.mesh.marching_cubes(np.zeros((3, 3), dtype=np.uint8)) + with pytest.raises(ValueError, match="at least 2x2x2"): + bic.mesh.marching_cubes(np.zeros((1, 2, 2), dtype=np.uint8)) + with pytest.raises(ValueError, match="Surface level"): + bic.mesh.marching_cubes(_interior_box(), 2.0) + with pytest.raises(ValueError, match="method should"): + bic.mesh.marching_cubes(_interior_box(), 0.5, method="other") + with pytest.raises(ValueError, match="gradient_direction"): + bic.mesh.marching_cubes(_interior_box(), 0.5, gradient_direction="other") + with pytest.raises(ValueError, match="spacing"): + bic.mesh.marching_cubes(_interior_box(), 0.5, spacing=(1.0, 1.0)) + with pytest.raises(ValueError, match="positive and finite"): + bic.mesh.marching_cubes(_interior_box(), 0.5, spacing=(1.0, 0.0, 1.0)) + with pytest.raises(ValueError, match="positive and finite"): + bic.mesh.marching_cubes(_interior_box(), 0.5, spacing=(1.0, np.nan, 1.0)) + with pytest.raises(ValueError, match="step_size"): + bic.mesh.marching_cubes(_interior_box(), 0.5, step_size=0) + with pytest.raises(TypeError, match="step_size"): + bic.mesh.marching_cubes(_interior_box(), 0.5, step_size=1.5) + with pytest.raises(ValueError, match="finite"): + bic.mesh.marching_cubes(_interior_box(), np.nan) + with pytest.raises(ValueError, match="same shape"): + bic.mesh.marching_cubes(_interior_box(), 0.5, mask=np.ones((3, 3, 3), bool)) + with pytest.raises(RuntimeError, match="No surface"): + bic.mesh.marching_cubes(np.zeros((3, 3, 3), dtype=np.uint8), 0.0) + with pytest.raises(TypeError, match="real numeric"): + bic.mesh.marching_cubes(np.ones((3, 3, 3), dtype=np.complex64), 1.0) From c524c990bae80f8ce4220ca013325f20d928f749 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 11 Jul 2026 13:09:48 +0200 Subject: [PATCH 2/2] Improve implementation, tests, and benchmark --- MIGRATION_GUIDE.md | 10 +- development/mesh/PERFORMANCE_NOTES.md | 60 ++ development/mesh/benchmark_marching_cubes.py | 562 +++++++++++++------ development/mesh/check_marching_cubes.py | 73 ++- include/bioimage_cpp/mesh/marching_cubes.hxx | 98 +++- src/bindings/mesh.cxx | 17 +- src/bioimage_cpp/mesh/__init__.py | 28 +- tests/mesh/test_marching_cubes.py | 117 ++++ 8 files changed, 745 insertions(+), 220 deletions(-) diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 7b32ad7..2cbf5c3 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -1905,6 +1905,8 @@ Important details: numeric or boolean input dtype is accepted; complex inputs are rejected. - `method="lewiner"` resolves ambiguous cases and is the default; `method="lorensen"` selects the original 256-case algorithm. +- `spacing` accepts either one positive finite scalar for isotropic data or a + length-three `(z, y, x)` sequence. - Normals and local-range values follow scikit-image semantics. As in scikit-image, `gradient_direction` reverses face winding without changing normals, and anisotropic spacing scales vertices without transforming @@ -1914,8 +1916,12 @@ Important details: foreground-positive segmentation masks. The iso-level is determined from the original unpadded volume. - Spacing entries must be positive and finite, and faces remain `int32` when - degenerate faces are removed. These validations/dtype choices are - intentional differences from scikit-image edge cases. + degenerate faces are removed. Duplicate vertices in a collapsed face are + merged transitively with the first vertex as representative, and faces that + still collapse after remapping are discarded. This guarantees in-range face + indices and intentionally avoids a rare scikit-image negative-index + remapping quirk. These validation and cleanup choices are intentional + differences from scikit-image edge cases. See `development/mesh/check_marching_cubes.py` for reference comparisons and `development/mesh/benchmark_marching_cubes.py` for reproducible timings. diff --git a/development/mesh/PERFORMANCE_NOTES.md b/development/mesh/PERFORMANCE_NOTES.md index 0537ff3..9520ba6 100644 --- a/development/mesh/PERFORMANCE_NOTES.md +++ b/development/mesh/PERFORMANCE_NOTES.md @@ -1,5 +1,65 @@ # Marching cubes performance +## Codex-Sol / Claude consolidation (2026-07-11) + +The consolidation pass kept the Codex-Sol float32 MC33 kernel, exact +scikit-image normal/value semantics, and zero-copy NumPy handoff. It moved +axis conversion and winding into triangle emission, added robust transitive +degenerate-vertex merging with the shared `UnionFind`, and broadened the +correctness and benchmark suites. Scalar spacing is now accepted as an +isotropic convenience. + +The benchmark now has a backward-compatible single-size mode and a full +scaling suite matching the independent implementation comparison: + +```bash +python development/mesh/benchmark_marching_cubes.py --size medium +python development/mesh/benchmark_marching_cubes.py \ + --suite scaling --repeats 5 --batches 1 --warmup 1 --memory \ + --json /tmp/marching-cubes.json +``` + +The scaling suite alternates paired bioimage-cpp/scikit-image calls and covers +Lewiner spheres/scalar fields through 512³, dense 10%-foreground masks through +256³, and all three Lorensen workloads at 128³. Fresh subprocesses measure +peak RSS before the timing process allocates any large volumes. + +Same-harness before/after results on the machine described below: + +| case | before ms | after ms | change | scikit-image / after | +|---|---:|---:|---:|---:| +| sphere, Lewiner, 512³ | 1,898.68 | 1,786.24 | -5.9% | 1.06× | +| dense mask, Lewiner, 256³ | 2,119.26 | 2,063.38 | -2.6% | 1.98× | +| scalar field, Lewiner, 512³ | 2,268.58 | 2,295.50 | +1.2% | 1.30× | +| sphere, Lorensen, 128³ | 30.31 | 28.39 | -6.3% | 1.31× | +| dense mask, Lorensen, 128³ | 192.19 | 212.75 | +10.7% | 2.16× | +| scalar field, Lorensen, 128³ | 48.72 | 49.80 | +2.2% | 1.84× | + +Absolute timings moved with CPU state: across all 15 cases the geometric-mean +wall time improved 2.4%, while normalization by each paired scikit-image time +showed a 1.1% improvement. Targeted same-session A/B repeats against the exact +pre-change commit resolved the apparent small-scalar and dense-Lorensen +regressions: scalar 64³ retained a ~2.04–2.06× paired speedup, while dense +Lorensen improved from ~2.01–2.02× to ~2.08×. No reproducible regression above +3% remained. + +Peak RSS was unchanged within measurement noise: + +| case | before MiB | after MiB | scikit-image MiB | +|---|---:|---:|---:| +| sphere 512³ | 725.8 | 725.7 | 734.6 | +| dense mask 256³ | 643.5 | 643.5 | 1,108.6 | +| scalar field 512³ | 703.3 | 703.4 | 843.3 | + +Correctness validation finished with 29 mesh tests, 1,024 full-suite tests, +all 254 nontrivial binary cube configurations for both methods, 1,024 random +scalar cubes for both methods, and 1,000 randomized multi-cube mask/stride/ +degeneracy cases. Twelve `allow_degenerate=False` cases intentionally differed +from scikit-image's negative-index remapping quirk; all returned faces were +valid and contained no collapsed-coordinate triangles. + +## Two-slice face-cache optimization + `bic.mesh.marching_cubes` remains deterministic and single-threaded. The optimization pass replaced the volume-growing edge hash map with Lewiner's two-slice face cache. Several smaller candidates were measured first and diff --git a/development/mesh/benchmark_marching_cubes.py b/development/mesh/benchmark_marching_cubes.py index dbd5ecc..dbdb985 100644 --- a/development/mesh/benchmark_marching_cubes.py +++ b/development/mesh/benchmark_marching_cubes.py @@ -1,48 +1,72 @@ """Benchmark bioimage-cpp marching cubes against scikit-image. -The benchmark runs a direct parity preflight before timing every workload. +The default ``single`` suite preserves the historical one-size benchmark. The +``scaling`` suite reproduces the independent comparison matrix used to select +the implementation: full Lewiner scaling, representative Lorensen cases, and +optional fresh-process peak-memory measurements. Examples -------- -python development/mesh/benchmark_marching_cubes.py --size medium --repeats 7 -python development/mesh/benchmark_marching_cubes.py --size large --method lorensen +python development/mesh/benchmark_marching_cubes.py --size medium +python development/mesh/benchmark_marching_cubes.py --suite scaling --repeats 5 --batches 1 +python development/mesh/benchmark_marching_cubes.py --suite scaling --memory --json /tmp/mc.json """ from __future__ import annotations import argparse +import gc import json +import os import platform +import random from statistics import median +import subprocess import sys from time import perf_counter import numpy as np import skimage +try: + import resource +except ImportError: # pragma: no cover - unavailable on Windows + resource = None + import bioimage_cpp as bic from _marching_cubes_reference import assert_mesh_matches, reference_marching_cubes +WORKLOADS = ("binary_sphere", "dense_binary_mask", "scalar_field") +METHODS = ("lewiner", "lorensen") +MEMORY_CASES = ( + ("binary_sphere", 512), + ("dense_binary_mask", 256), + ("scalar_field", 512), +) + + def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--suite", choices=("single", "scaling"), default="single") parser.add_argument("--size", choices=("small", "medium", "large"), default="medium") - parser.add_argument("--method", choices=("lewiner", "lorensen", "all"), default="all") - parser.add_argument( - "--workload", - choices=("binary_sphere", "dense_binary_mask", "scalar_field", "all"), - default="all", - ) - parser.add_argument( - "--backend", choices=("both", "bic", "skimage"), default="both" - ) + parser.add_argument("--method", choices=(*METHODS, "all"), default="all") + parser.add_argument("--workload", choices=(*WORKLOADS, "all"), default="all") + parser.add_argument("--backend", choices=("both", "bic", "skimage"), default="both") parser.add_argument("--repeats", type=int, default=7) parser.add_argument("--warmup", type=int, default=1) parser.add_argument("--batches", type=int, default=3) - parser.add_argument("--seed", type=int, default=20260709) + parser.add_argument("--seed", type=int, default=20260710) + parser.add_argument("--memory", action="store_true", help="run fresh-process largest-case RSS probes") + parser.add_argument("--memory-only", action="store_true", help="run only the fresh-process RSS probes") parser.add_argument("--json", default="", help="optional JSON result path") - parser.add_argument("--baseline", default="", help="optional prior benchmark JSON for bic relative changes") + parser.add_argument("--baseline", default="", help="optional prior JSON for bic relative changes") + # Private subprocess interface used by --memory. + parser.add_argument("--memory-worker", action="store_true", help=argparse.SUPPRESS) + parser.add_argument("--memory-backend", choices=("bic", "skimage"), help=argparse.SUPPRESS) + parser.add_argument("--memory-workload", choices=WORKLOADS, help=argparse.SUPPRESS) + parser.add_argument("--memory-size", type=int, help=argparse.SUPPRESS) return parser.parse_args() @@ -50,51 +74,125 @@ def shape_for(size: str) -> tuple[int, int, int]: return {"small": (48, 48, 48), "medium": (96, 96, 96), "large": (160, 160, 160)}[size] -def workloads(shape: tuple[int, int, int], seed: int): - z, y, x = np.ogrid[: shape[0], : shape[1], : shape[2]] - center = np.asarray(shape, dtype=np.float32) / 2.0 - radius = min(shape) * 0.28 - sphere = ((z - center[0]) ** 2 + (y - center[1]) ** 2 + (x - center[2]) ** 2 < radius**2).astype(np.uint8) - rng = np.random.default_rng(seed) - dense_mask = (rng.random(shape) < 0.10).astype(np.uint8) - zf = z.astype(np.float32) / max(shape[0] - 1, 1) - yf = y.astype(np.float32) / max(shape[1] - 1, 1) - xf = x.astype(np.float32) / max(shape[2] - 1, 1) - scalar = ( - np.sin(4.0 * np.pi * zf) - + np.cos(6.0 * np.pi * yf) - + np.sin(5.0 * np.pi * xf) - ).astype(np.float32) - return [ - ("binary_sphere", sphere, 0.5), - ("dense_binary_mask", dense_mask, 0.5), - ("scalar_field", scalar, 0.0), - ] - - -def time_call(function, repeats: int, warmup: int, batches: int) -> dict[str, object]: - raw_timings = [] - batch_medians = [] - for _ in range(batches): - for _ in range(warmup): - function() - timings = [] - for _ in range(repeats): - start = perf_counter() - function() - timings.append(perf_counter() - start) - raw_timings.extend(timings) - batch_medians.append(median(timings)) +def make_workload(name: str, n: int, seed: int) -> tuple[np.ndarray, float]: + shape = (n, n, n) + if name == "binary_sphere": + # Build slice-wise so the 512^3 memory probe is not polluted by large + # temporary broadcast arrays before the measured call. + y, x = np.ogrid[:n, :n] + center = np.float32((n - 1) / 2.0) + radius = np.float32(0.28 * n) + plane_distance = (y - center) ** 2 + (x - center) ** 2 + volume = np.empty(shape, dtype=np.uint8) + for z in range(n): + volume[z] = ( + (np.float32(z) - center) ** 2 + plane_distance <= radius**2 + ) + return volume, 0.5 + if name == "dense_binary_mask": + rng = np.random.default_rng(seed + n) + return ( + (rng.random(shape, dtype=np.float32) < np.float32(0.10)).astype(np.uint8), + 0.5, + ) + if name == "scalar_field": + q = np.linspace(0.0, 1.0, n, dtype=np.float32) + z = np.sin(np.float32(4.0 * np.pi) * q)[:, None, None] + y = np.cos(np.float32(6.0 * np.pi) * q)[None, :, None] + x = np.sin(np.float32(5.0 * np.pi) * q)[None, None, :] + return np.asarray(z + y + x, dtype=np.float32, order="C"), 0.0 + raise ValueError(f"unknown workload: {name}") + + +def case_specs(args: argparse.Namespace) -> list[tuple[str, int, str]]: + workloads = WORKLOADS if args.workload == "all" else (args.workload,) + methods = METHODS if args.method == "all" else (args.method,) + if args.suite == "single": + n = shape_for(args.size)[0] + return [(workload, n, method) for workload in workloads for method in methods] + + specs = [] + if "lewiner" in methods: + for workload in workloads: + sizes = (64, 128, 192, 256) if workload == "dense_binary_mask" else (64, 128, 256, 512) + specs.extend((workload, n, "lewiner") for n in sizes) + if "lorensen" in methods: + specs.extend((workload, 128, "lorensen") for workload in workloads) + return specs + + +def call_bic(volume: np.ndarray, level: float, method: str): + return bic.mesh.marching_cubes(volume, level, method=method) + + +def call_skimage(volume: np.ndarray, level: float, method: str): + return reference_marching_cubes(volume, level, method=method) + + +def time_once(function, volume, level, method) -> float: + start = perf_counter() + result = function(volume, level, method) + elapsed = perf_counter() - start + del result + return elapsed + + +def statistics(samples: list[float], batch_medians: list[float]) -> dict[str, object]: return { - "raw_s": raw_timings, + "raw_s": samples, "batch_medians_s": batch_medians, "median_s": median(batch_medians), - "p10_s": float(np.percentile(raw_timings, 10)), - "p90_s": float(np.percentile(raw_timings, 90)), - "min_s": min(raw_timings), + "min_s": min(samples), + "q25_s": float(np.percentile(samples, 25)), + "q75_s": float(np.percentile(samples, 75)), + "p10_s": float(np.percentile(samples, 10)), + "p90_s": float(np.percentile(samples, 90)), } +def time_backends( + volume: np.ndarray, + level: float, + method: str, + backend: str, + repeats: int, + warmup: int, + batches: int, + seed: int, +) -> tuple[dict[str, object] | None, dict[str, object] | None]: + functions = {} + if backend in ("both", "bic"): + functions["bic"] = call_bic + if backend in ("both", "skimage"): + functions["skimage"] = call_skimage + samples = {name: [] for name in functions} + batch_values = {name: [] for name in functions} + rng = random.Random(seed) + gc.collect() + gc.disable() + try: + for batch in range(batches): + for _ in range(warmup): + for name, function in functions.items(): + time_once(function, volume, level, method) + current = {name: [] for name in functions} + for _ in range(repeats): + order = list(functions) + rng.shuffle(order) + for name in order: + elapsed = time_once(functions[name], volume, level, method) + samples[name].append(elapsed) + current[name].append(elapsed) + for name in functions: + batch_values[name].append(median(current[name])) + finally: + gc.enable() + return ( + None if "bic" not in functions else statistics(samples["bic"], batch_values["bic"]), + None if "skimage" not in functions else statistics(samples["skimage"], batch_values["skimage"]), + ) + + def load_baseline(path: str) -> dict[tuple[str, str, tuple[int, int, int]], float]: if not path: return {} @@ -103,131 +201,269 @@ def load_baseline(path: str) -> dict[tuple[str, str, tuple[int, int, int]], floa return { (row["workload"], row["method"], tuple(row["shape"])): row["bic_median_s"] for row in rows + if row.get("bic_median_s") is not None + } + + +def validate_preflights(specs: list[tuple[str, int, str]], seed: int) -> None: + combinations = sorted({(workload, method) for workload, _, method in specs}) + for workload, method in combinations: + volume, level = make_workload(workload, 64, seed) + actual = call_bic(volume, level, method) + reference = call_skimage(volume, level, method) + assert_mesh_matches(actual, reference) + + +def scaling_exponents(rows: list[dict[str, object]]) -> list[dict[str, object]]: + output = [] + combinations = sorted({(row["workload"], row["method"]) for row in rows}) + for workload, method in combinations: + selected = sorted( + (row for row in rows if row["workload"] == workload and row["method"] == method), + key=lambda row: row["voxel_count"], + ) + if len(selected) < 3: + continue + selected = selected[-3:] + x = np.log([row["voxel_count"] for row in selected]) + entry = {"workload": workload, "method": method, "sizes": [row["shape"][0] for row in selected]} + for backend, field in (("bic", "bic_median_s"), ("skimage", "reference_median_s")): + values = [row[field] for row in selected] + entry[f"{backend}_exponent"] = None if any(value is None for value in values) else float(np.polyfit(x, np.log(values), 1)[0]) + output.append(entry) + return output + + +def current_rss_kib() -> int: + status = "/proc/self/status" + if not os.path.exists(status): + return 0 + with open(status) as file: + for line in file: + if line.startswith("VmRSS:"): + return int(line.split()[1]) + return 0 + + +def memory_worker(args: argparse.Namespace) -> int: + if resource is None: + raise RuntimeError("memory worker requires the Unix resource module") + if args.memory_backend is None or args.memory_workload is None or args.memory_size is None: + raise SystemExit("memory worker requires backend, workload, and size") + volume, level = make_workload(args.memory_workload, args.memory_size, args.seed) + gc.collect() + rss_before = current_rss_kib() + peak_before = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + function = call_bic if args.memory_backend == "bic" else call_skimage + start = perf_counter() + result = function(volume, level, "lewiner") + elapsed = perf_counter() - start + peak_after = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + payload = { + "backend": args.memory_backend, + "workload": args.memory_workload, + "method": "lewiner", + "shape": list(volume.shape), + "elapsed_s": elapsed, + "input_nbytes": volume.nbytes, + "output_nbytes": sum(np.asarray(array).nbytes for array in result), + "vertices": len(result[0]), + "faces": len(result[1]), + "rss_before_kib": rss_before, + "rss_after_kib": current_rss_kib(), + "peak_before_kib": peak_before, + "peak_after_kib": peak_after, + "incremental_peak_kib": max(0, peak_after - peak_before), + } + print(json.dumps(payload)) + return 0 + + +def run_memory_probes(args: argparse.Namespace) -> list[dict[str, object]]: + if platform.system() != "Linux" or resource is None: + raise RuntimeError("--memory currently requires Linux /proc and ru_maxrss semantics") + backends = ("bic", "skimage") if args.backend == "both" else (args.backend,) + rows = [] + script = os.path.abspath(__file__) + for workload, size in MEMORY_CASES: + if args.workload != "all" and args.workload != workload: + continue + for backend in backends: + command = [ + sys.executable, + script, + "--memory-worker", + "--memory-backend", backend, + "--memory-workload", workload, + "--memory-size", str(size), + "--seed", str(args.seed), + ] + completed = subprocess.run(command, check=True, capture_output=True, text=True) + rows.append(json.loads(completed.stdout)) + return rows + + +def environment() -> dict[str, object]: + return { + "python": sys.version, + "platform": platform.platform(), + "numpy": np.__version__, + "scikit_image": skimage.__version__, + "cpu_count": os.cpu_count(), + "thread_environment": { + key: os.environ.get(key) + for key in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "NUMEXPR_NUM_THREADS") + }, } def main() -> int: args = parse_args() + if args.memory_worker: + return memory_worker(args) if args.repeats < 1 or args.warmup < 0 or args.batches < 1: raise SystemExit("repeats and batches must be >= 1 and warmup must be >= 0") - methods = ("lewiner", "lorensen") if args.method == "all" else (args.method,) - shape = shape_for(args.size) - n_voxels = int(np.prod(shape)) + if args.memory_only and not args.memory: + raise SystemExit("--memory-only requires --memory") + # Launch memory workers before this process has allocated any large benchmark + # volume: Linux preserves the parent's RSS high-water mark across fork/exec. + memory_rows = run_memory_probes(args) if args.memory else [] + if args.memory_only: + payload = { + "suite": args.suite, + "seed": args.seed, + "backend": args.backend, + "environment": environment(), + "results": [], + "scaling": [], + "memory": memory_rows, + } + for row in memory_rows: + print( + f"{row['workload']}/{row['backend']}/{row['shape'][0]}^3: " + f"peak={row['peak_after_kib'] / 1024:.1f} MiB " + f"increment={row['incremental_peak_kib'] / 1024:.1f} MiB" + ) + if args.json: + with open(args.json, "w") as file: + json.dump(payload, file, indent=2) + return 0 + specs = case_specs(args) baseline = load_baseline(args.baseline) + validate_preflights(specs, args.seed) + print( + f"suite={args.suite} cases={len(specs)} repeats={args.repeats} " + f"batches={args.batches} preflight=OK" + ) + print( + f"{'workload/method/size':<39} {'V':>10} {'F':>11} {'bic ms':>11} " + f"{'IQR ms':>9} {'skimage ms':>12} {'speed':>9} {'Mvox/s':>10} {'delta':>9}" + ) + print("-" * 128) rows = [] - - print(f"shape={shape} voxels={n_voxels} repeats={args.repeats} batches={args.batches}") - print(f"{'workload/method':<28} {'V':>8} {'F':>8} {'bic ms':>11} {'p10-p90 ms':>15} {'skimage ms':>12} {'speed':>9} {'Mvox/s':>10} {'delta':>9}") - print("-" * 125) - selected_workloads = workloads(shape, args.seed) - if args.workload != "all": - selected_workloads = [row for row in selected_workloads if row[0] == args.workload] - for workload, volume, level in selected_workloads: - for method in methods: - kwargs = {"method": method} - start = perf_counter() - actual = bic.mesh.marching_cubes(volume, level, **kwargs) - bic_first_call = perf_counter() - start - start = perf_counter() - reference = reference_marching_cubes(volume, level, **kwargs) - reference_first_call = perf_counter() - start - assert_mesh_matches(actual, reference) - bic_times = None - if args.backend in ("both", "bic"): - bic_times = time_call( - lambda: bic.mesh.marching_cubes(volume, level, **kwargs), - args.repeats, - args.warmup, - args.batches, - ) - reference_times = None - if args.backend in ("both", "skimage"): - reference_times = time_call( - lambda: reference_marching_cubes(volume, level, **kwargs), - args.repeats, - args.warmup, - args.batches, - ) - bic_median = None if bic_times is None else bic_times["median_s"] - reference_median = ( - None if reference_times is None else reference_times["median_s"] - ) - baseline_time = baseline.get((workload, method, shape)) - relative_change = ( - None - if baseline_time is None or bic_median is None - else bic_median / baseline_time - 1.0 - ) - speedup = ( - None - if bic_median is None or reference_median is None - else reference_median / bic_median + for case_index, (workload, size, method) in enumerate(specs): + volume, level = make_workload(workload, size, args.seed) + start = perf_counter() + actual = call_bic(volume, level, method) + bic_first_call = perf_counter() - start + start = perf_counter() + reference = call_skimage(volume, level, method) + reference_first_call = perf_counter() - start + counts_match = len(actual[0]) == len(reference[0]) and len(actual[1]) == len(reference[1]) + valid_faces = bool(np.all((actual[1] >= 0) & (actual[1] < len(actual[0])))) + finite = all(np.all(np.isfinite(array)) for array in (actual[0], actual[2], actual[3])) + if not counts_match or not valid_faces or not finite: + raise AssertionError( + f"large-case validation failed for {workload}/{method}/{size}: " + f"counts={counts_match}, faces={valid_faces}, finite={finite}" ) - row = { - "workload": workload, - "method": method, - "shape": shape, - "vertices": len(actual[0]), - "faces": len(actual[1]), - "bic_first_call_s": bic_first_call, - "reference_first_call_s": reference_first_call, - "bic_raw_s": None if bic_times is None else bic_times["raw_s"], - "bic_batch_medians_s": None if bic_times is None else bic_times["batch_medians_s"], - "bic_median_s": bic_median, - "bic_min_s": None if bic_times is None else bic_times["min_s"], - "bic_p10_s": None if bic_times is None else bic_times["p10_s"], - "bic_p90_s": None if bic_times is None else bic_times["p90_s"], - "reference_raw_s": None if reference_times is None else reference_times["raw_s"], - "reference_batch_medians_s": None if reference_times is None else reference_times["batch_medians_s"], - "reference_median_s": reference_median, - "reference_min_s": None if reference_times is None else reference_times["min_s"], - "speedup": speedup, - "mvox_per_s": None if bic_median is None else n_voxels / bic_median / 1e6, - "baseline_bic_median_s": baseline_time, - "relative_change": relative_change, - } - rows.append(row) - delta = "-" if row["relative_change"] is None else f"{100.0 * row['relative_change']:+.1f}%" - bic_text = "-" if bic_median is None else f"{bic_median * 1e3:>11.2f}" - spread_text = ( - " - " - if bic_times is None - else f"{bic_times['p10_s'] * 1e3:>6.2f}-{bic_times['p90_s'] * 1e3:>6.2f}" - ) - reference_text = ( - "-" if reference_median is None else f"{reference_median * 1e3:>12.2f}" - ) - speed_text = "-" if speedup is None else f"{speedup:>8.2f}x" - throughput_text = ( - "-" if row["mvox_per_s"] is None else f"{row['mvox_per_s']:>10.2f}" + bic_times, reference_times = time_backends( + volume, + level, + method, + args.backend, + args.repeats, + args.warmup, + args.batches, + args.seed + case_index, + ) + bic_median = None if bic_times is None else bic_times["median_s"] + reference_median = None if reference_times is None else reference_times["median_s"] + baseline_time = baseline.get((workload, method, tuple(volume.shape))) + relative_change = None if baseline_time is None or bic_median is None else bic_median / baseline_time - 1.0 + speedup = None if bic_median is None or reference_median is None else reference_median / bic_median + mvox = None if bic_median is None else volume.size / bic_median / 1e6 + row = { + "workload": workload, + "method": method, + "shape": list(volume.shape), + "voxel_count": int(volume.size), + "vertices": len(actual[0]), + "faces": len(actual[1]), + "counts_match": counts_match, + "valid_faces": valid_faces, + "finite_outputs": finite, + "bic_first_call_s": bic_first_call, + "reference_first_call_s": reference_first_call, + "bic": bic_times, + "skimage": reference_times, + "bic_median_s": bic_median, + "reference_median_s": reference_median, + "speedup": speedup, + "mvox_per_s": mvox, + "baseline_bic_median_s": baseline_time, + "relative_change": relative_change, + } + rows.append(row) + bic_text = "-" if bic_median is None else f"{bic_median * 1e3:.2f}" + iqr_text = "-" if bic_times is None else f"{(bic_times['q75_s'] - bic_times['q25_s']) * 1e3:.2f}" + reference_text = "-" if reference_median is None else f"{reference_median * 1e3:.2f}" + speed_text = "-" if speedup is None else f"{speedup:.2f}x" + throughput_text = "-" if mvox is None else f"{mvox:.2f}" + delta_text = "-" if relative_change is None else f"{relative_change * 100:+.1f}%" + name = f"{workload}/{method}/{size}^3" + print( + f"{name:<39} {len(actual[0]):>10} {len(actual[1]):>11} {bic_text:>11} " + f"{iqr_text:>9} {reference_text:>12} {speed_text:>9} " + f"{throughput_text:>10} {delta_text:>9}" + ) + del actual, reference, volume + + exponents = scaling_exponents(rows) + if exponents: + print("\nscaling exponents (time ~ voxel_count^p, largest three sizes)") + for entry in exponents: + bic_exponent = entry["bic_exponent"] + reference_exponent = entry["skimage_exponent"] + bic_text = "-" if bic_exponent is None else f"{bic_exponent:.3f}" + reference_text = "-" if reference_exponent is None else f"{reference_exponent:.3f}" + print( + f" {entry['workload']}/{entry['method']}: " + f"bic={bic_text} skimage={reference_text}" ) + if memory_rows: + print("\npeak memory") + for row in memory_rows: print( - f"{workload + '/' + method:<28} {row['vertices']:>8} {row['faces']:>8} " - f"{bic_text:>11} {spread_text:>15} {reference_text:>12} " - f"{speed_text:>9} {throughput_text:>10} {delta:>9}" + f" {row['workload']}/{row['backend']}/{row['shape'][0]}^3: " + f"peak={row['peak_after_kib'] / 1024:.1f} MiB " + f"increment={row['incremental_peak_kib'] / 1024:.1f} MiB" ) + payload = { + "suite": args.suite, + "repeats": args.repeats, + "warmup": args.warmup, + "batches": args.batches, + "seed": args.seed, + "backend": args.backend, + "environment": environment(), + "results": rows, + "scaling": exponents, + "memory": memory_rows, + } if args.json: with open(args.json, "w") as file: - json.dump( - { - "shape": shape, - "repeats": args.repeats, - "warmup": args.warmup, - "batches": args.batches, - "seed": args.seed, - "workload": args.workload, - "backend": args.backend, - "environment": { - "python": sys.version, - "platform": platform.platform(), - "numpy": np.__version__, - "scikit_image": skimage.__version__, - }, - "results": rows, - }, - file, - indent=2, - ) + json.dump(payload, file, indent=2) return 0 diff --git a/development/mesh/check_marching_cubes.py b/development/mesh/check_marching_cubes.py index b8e3a24..8550a84 100644 --- a/development/mesh/check_marching_cubes.py +++ b/development/mesh/check_marching_cubes.py @@ -25,6 +25,12 @@ def parse_args() -> argparse.Namespace: default=1024, help="number of deterministic random scalar one-cube configurations to compare (default: 1024)", ) + parser.add_argument( + "--random-volumes", + type=int, + default=1000, + help="number of deterministic small multi-cube stress cases (default: 1000)", + ) return parser.parse_args() @@ -89,10 +95,67 @@ def check_random_scalar_cubes(n_cases: int) -> None: checked += 1 +def collapsed_coordinate_faces(vertices: np.ndarray, faces: np.ndarray) -> int: + triangles = vertices[faces] + collapsed = ( + np.all(triangles[:, 0] == triangles[:, 1], axis=1) + | np.all(triangles[:, 0] == triangles[:, 2], axis=1) + | np.all(triangles[:, 1] == triangles[:, 2], axis=1) + ) + return int(np.count_nonzero(collapsed)) + + +def check_random_volumes(n_cases: int) -> int: + """Stress masks, strides, and exact-level degeneracies on small volumes. + + Robust transitive duplicate merging intentionally differs from a rare + scikit-image negative-index remapping quirk. Such differences are counted + and reported, while invalid indices or collapsed output faces always fail. + """ + rng = np.random.default_rng(293841) + levels = np.asarray([0.5, 1.0, 1.5, 2.0, 2.5]) + parity_differences = 0 + for _ in range(n_cases): + shape = tuple(int(size) for size in rng.integers(3, 9, size=3)) + volume = rng.integers(0, 4, size=shape, dtype=np.uint8) + volume.flat[0] = 0 + volume.flat[-1] = 3 + level = float(rng.choice(levels)) + kwargs: dict[str, object] = { + "method": "lewiner" if rng.random() < 0.7 else "lorensen", + "allow_degenerate": bool(rng.random() < 0.5), + "step_size": 2 if min(shape) >= 5 and rng.random() < 0.25 else 1, + } + if rng.random() < 0.25: + kwargs["mask"] = rng.random(shape) < 0.75 + + try: + actual = bic.mesh.marching_cubes(volume, level, **kwargs) + except RuntimeError: + try: + reference_marching_cubes(volume, level, **kwargs) + except RuntimeError: + continue + raise AssertionError("bioimage-cpp found no surface but scikit-image did") + reference = reference_marching_cubes(volume, level, **kwargs) + if not np.all((actual[1] >= 0) & (actual[1] < len(actual[0]))): + raise AssertionError("marching cubes returned an invalid face index") + if not kwargs["allow_degenerate"]: + if collapsed_coordinate_faces(actual[0], actual[1]) != 0: + raise AssertionError("degenerate removal left a collapsed-coordinate face") + try: + assert_mesh_matches(actual, reference) + except AssertionError: + parity_differences += 1 + else: + assert_mesh_matches(actual, reference) + return parity_differences + + def main() -> int: args = parse_args() - if args.random_cases < 0: - raise SystemExit("random-cases must be >= 0") + if args.random_cases < 0 or args.random_volumes < 0: + raise SystemExit("random-cases and random-volumes must be >= 0") failed = False print(f"{'case':<26} {'vertices':>10} {'faces':>10} {'status':>8}") print("-" * 60) @@ -112,6 +175,12 @@ def main() -> int: print(f"{'all binary one-cube cases':<26} {'254':>10} {'508':>10} {'OK':>8}") check_random_scalar_cubes(args.random_cases) print(f"{'random scalar one-cube cases':<26} {args.random_cases:>10} {2 * args.random_cases:>10} {'OK':>8}") + differences = check_random_volumes(args.random_volumes) + detail = f"{differences} intentional reference differences" + print( + f"{'random multi-cube cases':<26} {args.random_volumes:>10} " + f"{differences:>10} {'OK':>8} {detail}" + ) except Exception as error: print(f"{'randomized parity':<26} {'-':>10} {'-':>10} {'FAIL':>8} {error}") failed = True diff --git a/include/bioimage_cpp/mesh/marching_cubes.hxx b/include/bioimage_cpp/mesh/marching_cubes.hxx index 701082d..1308e89 100644 --- a/include/bioimage_cpp/mesh/marching_cubes.hxx +++ b/include/bioimage_cpp/mesh/marching_cubes.hxx @@ -13,6 +13,7 @@ #include "bioimage_cpp/array_view.hxx" #include "bioimage_cpp/detail/profile.hxx" #include "bioimage_cpp/mesh/detail/mc33_luts.hxx" +#include "bioimage_cpp/util/union_find.hxx" #include #include @@ -34,10 +35,15 @@ enum class MarchingCubesMethod { Lorensen, }; +enum class GradientDirection { + Descent, + Ascent, +}; + struct MarchingCubesResult { // All vector-valued arrays are flat C-order arrays with a trailing size-3 - // component axis. Vertices/normals use the reference kernel's x/y/z order; - // the binding converts them to NumPy z/y/x order before returning them. + // component axis. Vertices and normals use NumPy z/y/x order; faces have + // already been oriented according to GradientDirection. std::vector vertices; std::vector faces; std::vector normals; @@ -169,8 +175,8 @@ constexpr std::array, 12> kEdgeRelativeZ{{ class Cell { public: - Cell(const int nx, const int ny) - : nx_(nx), ny_(ny), + Cell(const int nx, const int ny, const GradientDirection gradient_direction) + : nx_(nx), ny_(ny), gradient_direction_(gradient_direction), face_layer1_(static_cast(nx) * static_cast(ny) * 4, -1), face_layer2_(static_cast(nx) * static_cast(ny) * 4, -1) {} @@ -219,6 +225,7 @@ public: for (int corner = 0; corner < 3; ++corner) { add_face_from_edge(lut.get2(lut_index, triangle * 3 + corner)); } + orient_last_face(); } } @@ -228,6 +235,7 @@ public: for (int corner = 0; corner < 3; ++corner) { add_face_from_edge(lut.get3(lut_index, lut_index2, triangle * 3 + corner)); } + orient_last_face(); } } @@ -258,9 +266,9 @@ private: throw std::runtime_error("marching cubes produced too many vertices for int32 faces"); } const int index = static_cast(values_.size()); - vertices_.push_back(static_cast(x)); - vertices_.push_back(static_cast(y)); vertices_.push_back(static_cast(z)); + vertices_.push_back(static_cast(y)); + vertices_.push_back(static_cast(x)); normals_.insert(normals_.end(), {0.0F, 0.0F, 0.0F}); values_.push_back(0.0F); return index; @@ -275,9 +283,15 @@ private: void add_gradient(const int vertex, const double x, const double y, const double z) { const std::size_t base = static_cast(vertex) * 3; - normals_[base] += static_cast(x); + normals_[base] += static_cast(z); normals_[base + 1] += static_cast(y); - normals_[base + 2] += static_cast(z); + normals_[base + 2] += static_cast(x); + } + + void orient_last_face() { + if (gradient_direction_ == GradientDirection::Descent) { + std::swap(faces_[faces_.size() - 3], faces_[faces_.size() - 1]); + } } void add_gradient_from_corner(const int vertex, const int corner, const double strength) { @@ -418,6 +432,7 @@ private: int nx_ = 0; int ny_ = 0; + GradientDirection gradient_direction_ = GradientDirection::Descent; int x_ = 0; int y_ = 0; int z_ = 0; @@ -473,9 +488,7 @@ inline void select_mc33_tiling(const Luts &luts, Cell &cell, int case_, int conf inline void remove_degenerate_faces(MarchingCubesResult &result) { const std::size_t n_vertices = result.values.size(); - std::vector map(n_vertices); - for (std::size_t i = 0; i < n_vertices; ++i) map[i] = static_cast(i); - std::vector keep_face(result.faces.size() / 3, true); + util::UnionFind components(n_vertices); const auto equal_vertex = [&result](const std::int32_t a, const std::int32_t b) { const std::size_t first = static_cast(a) * 3; const std::size_t second = static_cast(b) * 3; @@ -483,30 +496,57 @@ inline void remove_degenerate_faces(MarchingCubesResult &result) { && result.vertices[first + 1] == result.vertices[second + 1] && result.vertices[first + 2] == result.vertices[second + 2]; }; - for (std::size_t face = 0; face < keep_face.size(); ++face) { + const auto merge_stably = [&components](const std::int32_t first, const std::int32_t second) { + const auto first_root = components.find(static_cast(first)); + const auto second_root = components.find(static_cast(second)); + if (first_root != second_root) { + components.merge_to( + std::min(first_root, second_root), std::max(first_root, second_root) + ); + } + }; + for (std::size_t face = 0; face < result.faces.size() / 3; ++face) { const std::int32_t a = result.faces[face * 3]; const std::int32_t b = result.faces[face * 3 + 1]; const std::int32_t c = result.faces[face * 3 + 2]; - if (equal_vertex(a, b)) { map[static_cast(a)] = map[static_cast(b)] = std::min(map[static_cast(a)], map[static_cast(b)]); keep_face[face] = false; } - if (equal_vertex(a, c)) { map[static_cast(a)] = map[static_cast(c)] = std::min(map[static_cast(a)], map[static_cast(c)]); keep_face[face] = false; } - if (equal_vertex(b, c)) { map[static_cast(b)] = map[static_cast(c)] = std::min(map[static_cast(b)], map[static_cast(c)]); keep_face[face] = false; } + if (equal_vertex(a, b)) merge_stably(a, b); + if (equal_vertex(a, c)) merge_stably(a, c); + if (equal_vertex(b, c)) merge_stably(b, c); + } + + std::vector roots(n_vertices); + for (std::size_t vertex = 0; vertex < n_vertices; ++vertex) { + roots[vertex] = components.find(static_cast(vertex)); } + std::vector compact(n_vertices, -1); MarchingCubesResult output; for (std::size_t vertex = 0; vertex < n_vertices; ++vertex) { - if (map[vertex] == static_cast(vertex)) { + if (roots[vertex] == vertex) { compact[vertex] = static_cast(output.values.size()); - output.vertices.insert(output.vertices.end(), result.vertices.begin() + static_cast(vertex * 3), result.vertices.begin() + static_cast(vertex * 3 + 3)); - output.normals.insert(output.normals.end(), result.normals.begin() + static_cast(vertex * 3), result.normals.begin() + static_cast(vertex * 3 + 3)); + const auto begin = static_cast(vertex * 3); + output.vertices.insert( + output.vertices.end(), + result.vertices.begin() + begin, + result.vertices.begin() + begin + 3 + ); + output.normals.insert( + output.normals.end(), + result.normals.begin() + begin, + result.normals.begin() + begin + 3 + ); output.values.push_back(result.values[vertex]); } } - for (std::size_t face = 0; face < keep_face.size(); ++face) { - if (!keep_face[face]) continue; - for (int corner = 0; corner < 3; ++corner) { - const auto old = static_cast(result.faces[face * 3 + static_cast(corner)]); - output.faces.push_back(compact[static_cast(map[old])]); - } + + for (std::size_t face = 0; face < result.faces.size() / 3; ++face) { + const auto a = roots[static_cast(result.faces[face * 3])]; + const auto b = roots[static_cast(result.faces[face * 3 + 1])]; + const auto c = roots[static_cast(result.faces[face * 3 + 2])]; + if (a == b || a == c || b == c) continue; + output.faces.push_back(compact[static_cast(a)]); + output.faces.push_back(compact[static_cast(b)]); + output.faces.push_back(compact[static_cast(c)]); } result = std::move(output); } @@ -769,14 +809,16 @@ inline void traverse_cells( } // namespace detail::marching_cubes -// Extract an isosurface from a 3-D float32 image. The core deliberately does -// not know about Python layout conventions beyond the C-order ArrayView; the -// binding handles argument validation, output axis order, and spacing. +// Extract an isosurface from a 3-D float32 image. Vertices and normals are +// returned in NumPy (z, y, x) axis order. `gradient_direction` controls face +// winding only, matching the public Python contract; spacing remains a binding +// concern because it does not affect the float32 extraction kernel. inline MarchingCubesResult marching_cubes( const ConstArrayView &volume, const double level, const int step_size, const MarchingCubesMethod method, + const GradientDirection gradient_direction, const ConstArrayView *mask, const bool allow_degenerate ) { @@ -803,7 +845,7 @@ inline MarchingCubesResult marching_cubes( } const auto &tables = *tables_ptr; std::optional cell; - cell.emplace(nx, ny); + cell.emplace(nx, ny, gradient_direction); { BIOIMAGE_PROFILE_SCOPE(profiler, "cell_traversal"); if (mask == nullptr) { diff --git a/src/bindings/mesh.cxx b/src/bindings/mesh.cxx index 2e9000c..7511ebb 100644 --- a/src/bindings/mesh.cxx +++ b/src/bindings/mesh.cxx @@ -113,6 +113,7 @@ nb::tuple marching_cubes_float32( level, step_size, classic ? mesh::MarchingCubesMethod::Lorensen : mesh::MarchingCubesMethod::Lewiner, + descent ? mesh::GradientDirection::Descent : mesh::GradientDirection::Ascent, mask_ptr, allow_degenerate ); @@ -121,22 +122,6 @@ nb::tuple marching_cubes_float32( throw std::runtime_error("No surface found at the given iso value."); } - { - BIOIMAGE_PROFILE_SCOPE(profiler, "orient_output"); - // The core mirrors the reference kernel's x/y/z convention. Public arrays - // use NumPy's z/y/x axis order, matching skimage.measure.marching_cubes. - for (std::size_t vertex = 0; vertex < result.values.size(); ++vertex) { - const std::size_t base = vertex * 3; - std::swap(result.vertices[base], result.vertices[base + 2]); - std::swap(result.normals[base], result.normals[base + 2]); - } - if (descent) { - for (std::size_t face = 0; face < result.faces.size(); face += 3) { - std::swap(result.faces[face], result.faces[face + 2]); - } - } - } - const std::size_t n_vertices = result.values.size(); const std::size_t n_faces = result.faces.size() / 3; FloatOutput vertices; diff --git a/src/bioimage_cpp/mesh/__init__.py b/src/bioimage_cpp/mesh/__init__.py index d0c0ffd..edb9dc3 100644 --- a/src/bioimage_cpp/mesh/__init__.py +++ b/src/bioimage_cpp/mesh/__init__.py @@ -10,13 +10,20 @@ from .. import _core -def _as_spacing(spacing: Sequence[float]) -> np.ndarray: - try: - values = np.asarray(spacing, dtype=np.float64) - except (TypeError, ValueError) as error: - raise ValueError("spacing must consist of three floats") from error +def _as_spacing(spacing: float | Sequence[float]) -> np.ndarray: + if np.isscalar(spacing): + try: + value = float(spacing) + except (TypeError, ValueError) as error: + raise ValueError("spacing must be a scalar or consist of three floats") from error + values = np.full(3, value, dtype=np.float64) + else: + try: + values = np.asarray(spacing, dtype=np.float64) + except (TypeError, ValueError) as error: + raise ValueError("spacing must be a scalar or consist of three floats") from error if values.shape != (3,): - raise ValueError("spacing must consist of three floats") + raise ValueError("spacing must be a scalar or consist of three floats") if not np.all(np.isfinite(values)) or np.any(values <= 0.0): raise ValueError("spacing entries must be positive and finite") return values @@ -26,7 +33,7 @@ def marching_cubes( volume: np.ndarray, level: float | None = None, *, - spacing: Sequence[float] = (1.0, 1.0, 1.0), + spacing: float | Sequence[float] = (1.0, 1.0, 1.0), gradient_direction: str = "descent", step_size: int = 1, allow_degenerate: bool = True, @@ -50,14 +57,17 @@ def marching_cubes( level: Iso-value to extract. ``None`` uses the midpoint of the data range. spacing: - Three physical spacings in NumPy axis order ``(z, y, x)``. + One isotropic spacing or three physical spacings in NumPy axis order + ``(z, y, x)``. gradient_direction: ``"descent"`` (the default) treats objects as values larger than the exterior; ``"ascent"`` reverses triangle winding. step_size: Sampling stride in voxels. Values above one produce a coarser mesh. allow_degenerate: - If false, remove faces with repeated vertex coordinates. + If false, transitively merge repeated vertex coordinates and remove + faces that collapse after remapping. Representatives preserve the + first vertex occurrence and all returned face indices are valid. method: ``"lewiner"`` or ``"lorensen"``. mask: diff --git a/tests/mesh/test_marching_cubes.py b/tests/mesh/test_marching_cubes.py index a143e12..91cf3e3 100644 --- a/tests/mesh/test_marching_cubes.py +++ b/tests/mesh/test_marching_cubes.py @@ -51,6 +51,47 @@ def _collapsed_coordinate_faces(vertices: np.ndarray, faces: np.ndarray) -> int: return int(np.count_nonzero(first_second | first_third | second_third)) +def _trilinear_sample(volume: np.ndarray, coordinates: np.ndarray) -> np.ndarray: + """Independently sample a scalar field at NumPy-order (z, y, x) coordinates.""" + image = np.asarray(volume, dtype=np.float64) + output = np.empty(len(coordinates), dtype=np.float64) + nz, ny, nx = image.shape + for index, (z, y, x) in enumerate(coordinates): + z0, y0, x0 = int(np.floor(z)), int(np.floor(y)), int(np.floor(x)) + z1, y1, x1 = min(z0 + 1, nz - 1), min(y0 + 1, ny - 1), min(x0 + 1, nx - 1) + fz, fy, fx = z - z0, y - y0, x - x0 + value = 0.0 + for cz, wz in ((z0, 1.0 - fz), (z1, fz)): + for cy, wy in ((y0, 1.0 - fy), (y1, fy)): + for cx, wx in ((x0, 1.0 - fx), (x1, fx)): + value += wz * wy * wx * image[cz, cy, cx] + output[index] = value + return output + + +def _transitive_degenerate_regression_volume() -> np.ndarray: + """Reproduce a duplicate-vertex chain that previously yielded face index -1.""" + rng = np.random.default_rng(293841) + levels = np.asarray([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]) + for _ in range(16): + shape = tuple(int(size) for size in rng.integers(3, 9, size=3)) + volume = rng.integers(0, 4, size=shape, dtype=np.uint8) + volume.flat[0] = 0 + volume.flat[-1] = 3 + level = float(rng.choice(levels)) + method = "lewiner" if rng.random() < 0.7 else "lorensen" + allow_degenerate = bool(rng.random() < 0.5) + step_size = 2 if min(shape) >= 5 and rng.random() < 0.25 else 1 + if rng.random() < 0.25: + rng.random(shape) + assert shape == (6, 5, 8) + assert level == 1.0 + assert method == "lorensen" + assert not allow_degenerate + assert step_size == 1 + return volume + + def test_box_mesh_has_reference_dimensions_and_dtypes(): vertices, faces, normals, values = bic.mesh.marching_cubes(_interior_box(), 0.5) @@ -69,6 +110,38 @@ def test_box_mesh_has_reference_dimensions_and_dtypes(): np.testing.assert_array_equal(values, np.ones_like(values)) +@pytest.mark.parametrize("method", ["lewiner", "lorensen"]) +def test_method_output_invariants(method): + volume = _interior_ball() + vertices, faces, normals, values = bic.mesh.marching_cubes( + volume, 0.5, method=method + ) + + assert vertices.ndim == 2 and vertices.shape[1] == 3 + assert faces.ndim == 2 and faces.shape[1] == 3 + assert normals.shape == vertices.shape + assert values.shape == (len(vertices),) + assert vertices.dtype == normals.dtype == values.dtype == np.float32 + assert faces.dtype == np.int32 + assert np.all((faces >= 0) & (faces < len(vertices))) + np.testing.assert_allclose(np.linalg.norm(normals, axis=1), 1.0, atol=1e-4) + assert vertices.min() >= 0.0 + for axis, size in enumerate(volume.shape): + assert vertices[:, axis].max() <= size - 1 + + +def test_lorensen_vertices_lie_on_isosurface(): + rng = np.random.default_rng(0) + volume = rng.random((12, 13, 14), dtype=np.float32) + level = 0.5 + vertices, _, _, _ = bic.mesh.marching_cubes( + volume, level, method="lorensen" + ) + np.testing.assert_allclose( + _trilinear_sample(volume, vertices), level, atol=1e-4 + ) + + def test_output_arrays_are_contiguous_writable_and_survive_collection(): outputs = bic.mesh.marching_cubes(_interior_box(), 0.5) expected = tuple(output.copy() for output in outputs) @@ -113,6 +186,15 @@ def test_spacing_preserves_mesh_and_uses_float64_vertices(): np.testing.assert_array_equal(values, unit_values) +def test_scalar_spacing_is_broadcast_to_all_axes(): + volume = _interior_box() + scalar = bic.mesh.marching_cubes(volume, 0.5, spacing=2.0) + explicit = bic.mesh.marching_cubes(volume, 0.5, spacing=(2.0, 2.0, 2.0)) + for actual, expected in zip(scalar, explicit, strict=True): + np.testing.assert_array_equal(actual, expected) + assert scalar[0].dtype == np.float64 + + @pytest.mark.parametrize("dtype", [np.bool_, np.uint16, np.float64]) def test_numeric_dtype_conversion_matches_uint8_input(dtype): volume = _interior_box() @@ -180,6 +262,18 @@ def test_mask_and_step_size_are_honoured(): assert values.shape == (len(vertices),) +def test_mask_reduces_surface_and_step_size_coarsens_mesh(): + rng = np.random.default_rng(4) + volume = rng.random((20, 22, 24), dtype=np.float32) + mask = np.ones(volume.shape, dtype=bool) + mask[:, :, : volume.shape[2] // 2] = False + fine = bic.mesh.marching_cubes(volume, 0.5) + masked = bic.mesh.marching_cubes(volume, 0.5, mask=mask) + coarse = bic.mesh.marching_cubes(volume, 0.5, step_size=2) + assert len(masked[1]) < len(fine[1]) + assert len(coarse[1]) < len(fine[1]) + + @pytest.mark.parametrize("method", ["lewiner", "lorensen"]) @pytest.mark.parametrize("step_size", [1, 2, 3]) def test_vertex_cache_with_step_size_and_mask_is_deterministic(method, step_size): @@ -241,6 +335,29 @@ def test_allow_degenerate_false_removes_collapsed_faces(): assert removed[1].dtype == np.int32 +def test_transitive_degenerate_vertex_merges_have_valid_indices(): + volume = _transitive_degenerate_regression_volume() + first = bic.mesh.marching_cubes( + volume, + 1.0, + method="lorensen", + allow_degenerate=False, + ) + second = bic.mesh.marching_cubes( + volume, + 1.0, + method="lorensen", + allow_degenerate=False, + ) + + assert first[0].shape == (251, 3) + assert first[1].shape == (360, 3) + assert np.all((first[1] >= 0) & (first[1] < len(first[0]))) + assert _collapsed_coordinate_faces(first[0], first[1]) == 0 + for actual, expected in zip(first, second, strict=True): + np.testing.assert_array_equal(actual, expected) + + def test_default_level_is_midpoint_of_original_volume_when_padding(): volume = np.full((5, 5, 5), 2.0, dtype=np.float32) volume[1:4, 1:4, 1:4] = 4.0