diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index f4e7403..878a075 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -2188,6 +2188,98 @@ Important differences: pairwise distance matrix; threshold the distance map first to keep the candidate count modest. +## scikit-fmm + +`scikit-fmm` computes geodesic distances on regular grids with the fast +marching method. `bioimage-cpp` groups geodesic distances under `bic.distance` +with two operations — a **distance field** to a set of sources, and the full +**pairwise distance matrix** between a set of points — for two geometry types: +regular-grid **masks** (the scikit-fmm equivalent) and triangle **meshes** +(which scikit-fmm does not support; the reference for those is the exact MMP +algorithm in `pygeodesic`). + +Both are solved with a first-order fast marching method (Godunov/Eikonal on the +grid, Kimmel–Sethian on triangles). Masks match scikit-fmm's `order=1` scheme; +mesh distances are first-order approximations of the exact surface geodesics. + +scikit-fmm: + +```python +import numpy as np +import skfmm + +# Distance from a set of seed voxels, constrained to a mask (domain). +phi = np.ones(mask.shape) +phi[tuple(sources.T)] = -1 # zero contour marks the seeds +phi = np.ma.MaskedArray(phi, ~mask) # obstacles / outside-domain +field = np.abs(skfmm.distance(phi, dx=sampling)) +# Weighted (travel-time) variant: +tt = skfmm.travel_time(phi, speed, dx=sampling) +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +# masks (regular grid); sources/points are (n, ndim) int64 voxel coordinates +field = bic.distance.geodesic_distance_field(mask, sources, sampling=None, + speed=None) # -> mask.shape, float64 +matrix = bic.distance.geodesic_distances(mask, points) # -> (N, N) float64 + +# optional per-axis gradient of the field (like vector_difference_transform) +field, grad = bic.distance.geodesic_distance_field(mask, sources, + return_gradient=True) +# grad: (*mask.shape, ndim) float32, d(field)/d(axis), points away from source +gradient = bic.distance.geodesic_gradient_field(mask, sources) # just the gradient + +# surfaces (triangle mesh); sources/points are 1-D int64 vertex indices +field = bic.distance.geodesic_distance_field_mesh(vertices, faces, sources, + speed=None) # -> (n_vertices,) float64 +matrix = bic.distance.geodesic_distances_mesh(vertices, faces, points) # -> (N, N) float64 +``` + +Name mapping: + +| scikit-fmm / pygeodesic | bioimage-cpp name | +| --- | --- | +| `skfmm.distance` / `skfmm.travel_time` (from seed voxels, masked) | `geodesic_distance_field` | +| gradient of the distance field (per-axis) | `geodesic_distance_field(..., return_gradient=True)` / `geodesic_gradient_field` | +| pairwise via repeated `skfmm.distance` | `geodesic_distances` | +| `pygeodesic … geodesicDistances(sources, None)` | `geodesic_distance_field_mesh` | +| `pygeodesic … geodesicDistances` (pairwise) | `geodesic_distances_mesh` | + +Important differences: + +- Explicit, geometry-specific functions rather than a level-set encoding: the + domain is passed directly as a `mask` (nonzero = inside) or as a + `(vertices, faces)` mesh, and the sources are given as coordinates / vertex + indices instead of being baked into a signed `phi`. +- Outputs are `float64`. Voxels outside the mask and points/vertices + unreachable from any source are `+inf`; pairwise matrices are symmetric with + a zero diagonal. +- `sampling` (per-axis voxel spacing, scalar or per-axis) maps to scikit-fmm's + `dx` and applies to masks only — meshes carry real vertex coordinates. +- `speed` is optional (`None` = unit-speed geodesic distance). When supplied it + generalizes to a weighted travel time, matching `skfmm.travel_time`; for + masks it has the mask's shape, for meshes it is per-vertex. +- `number_of_threads` follows the `bic.distance` convention (`1` default, + `0` = hardware concurrency); the pairwise solves parallelize over sources. +- `geodesic_distance_field(..., return_gradient=True)` also returns the per-axis + gradient `∇T` of the field (or use `geodesic_gradient_field` for the gradient + alone), analogous to `vector_difference_transform`. It is `float32` with shape + `(*mask.shape, ndim)` (channel-last, NumPy axis order); component `i` is + `d(field)/d(axis_i)`, pointing **away** from the nearest source with + `‖∇T‖ ≈ 1/speed`. Negate it to trace back toward the source — e.g. feed + `-grad` (transposed to channel-first) to `bic.flow.compute_flow_density`. It + is zero at sources, background, and unreachable voxels. Masks only. +- Mesh geodesics are surface (2-manifold) distances from the Kimmel–Sethian + triangle fast-marching method — a first-order approximation (a few % mean + error vs the exact `pygeodesic` MMP algorithm, larger on very obtuse + triangles, since obtuse-angle unfolding is not done yet). Like the mask + solver it slightly overestimates. See `development/distance/` for the + reference oracles and `benchmark_geodesic.py` for timings. + ## I/O and Build Dependencies `bioimage-cpp` intentionally does not replace nifty or affogato I/O helpers. diff --git a/development/distance/PERFORMANCE_NOTES.md b/development/distance/PERFORMANCE_NOTES.md new file mode 100644 index 0000000..3367569 --- /dev/null +++ b/development/distance/PERFORMANCE_NOTES.md @@ -0,0 +1,103 @@ +# Geodesic FMM performance notes + +Optimization of the two geodesic Fast Marching Method solvers, with the +before/after measurements that justify the changes. Both changes are pure +reorganizations of identical arithmetic — the solved fields are **bitwise +identical** to the pre-optimization solvers on every code path tested. + +## What changed + +- **Lead A — mask/grid solver** (`distance/detail/fast_marching.hxx`): the FMM + march no longer calls `detail::valid_offset_target` (one integer division + + modulo *per axis, per call*, ~42 calls / popped voxel in 3D). Each popped voxel + is decoded **once** with the new `detail::coords_from_index` (divide-and-subtract, + no modulo); axis neighbours are reached with a single `p ± strides_[d]` add and an + integer bounds compare, and the neighbour's coordinates are handed to + `solve_eikonal` so it never re-decodes. The `offsets_` `vector` is gone. +- **Lead B — mesh solver** (`distance/detail/mesh_fast_marching.hxx`): when a + vertex `v` freezes, each incident face relaxes its other corners from *that face + alone* instead of calling `update_vertex`, which rescanned every incident face of + every touched vertex. `dist_[w]` is a running minimum that already reflects the + other faces' contributions (applied when their corners froze), so the single-face + update reaches the same value. `update_vertex` was removed. Operands are passed to + `triangle_update` in face-storage order so the acute update is bit-identical, not + just equal to within rounding. + +## Measurement setup + +8-core machine, `powersave` governor (frequency-scaling noise present; `perf` +unavailable, `perf_event_paranoid=4`). Mitigations: 1 warmup call + 20 timed +repeats, **min** used as the headline estimator (least-interference), spread +reported. Single-source **field** solve, `--threads 1` — the single-source solve is +the whole kernel (pairwise = N × field). Baseline and optimized were built +separately; the **`mesh/field` case is a null control for Lead A** (its code is +untouched between the baseline and A-only builds) and moved <1%, bounding machine +drift below the measured mask gains. Lead A's mask numbers also reproduced within +~1% across two independent builds (A-only and A+B). + +Reproduce (from `development/distance/`): + +``` +python benchmark_geodesic.py --large --only field --threads 1 --repeats 20 --no-ref --json base.json +python benchmark_geodesic.py --xlarge --only field --threads 1 --repeats 20 --no-ref --json opt.json +python /path/to/compare.py base.json opt.json # or diff the JSON directly +``` + +## Wall-clock: baseline → optimized (single-thread field, min of 20) + +| case | baseline min | optimized min | speedup | +|---|---|---|---| +| mask2d/field (1024²) | 294.0 ms | 208.9 ms | **1.41×** | +| mask2d/field (1536²) | 685.6 ms | 492.7 ms | **1.39×** | +| mask3d/field (64·128²) | 745.1 ms | 370.2 ms | **2.01×** | +| mask3d/field (128³) | 1811.6 ms | 912.4 ms | **1.99×** | +| mesh/field (V=20000) | 23.2 ms | 7.9 ms | **2.94×** | +| mesh/field (V=40000) | 49.1 ms | 17.5 ms | **2.81×** | + +Lead A: ~2.0× in 3D (matches the div/mod arithmetic dominating 3 axes), ~1.4× in +2D. Lead B: ~2.8–2.9× on the mesh. Both exceed the pre-work predictions (A: 1.5–2× +3D; B: ~2×). + +## Mechanism confirmation (BIOIMAGE_PROFILE=ON, one field solve) + +Coarse phases in `solve()`: `pop` (heap pop, untouched by either change) vs `relax` +(the neighbour/face loop where the removed work lived). + +| solve | phase | baseline | optimized | change | +|---|---|---|---|---| +| mask 128³ | pop | 0.319 s | 0.316 s | unchanged | +| mask 128³ | relax | 1.700 s | 0.762 s | **2.23× smaller** | +| mesh 40k | pop | 0.0035 s | 0.0034 s | unchanged | +| mesh 40k | relax | 0.0500 s | 0.0163 s | **3.07× smaller** | + +The heap phase is invariant; the targeted `relax` phase shrank by 2.2× (mask) and +3.1× (mesh). `pop`'s *share* rose (15.8%→29.3% mask; 6.5%→17.1% mesh) only because +`relax` shrank around it — the gain came from exactly the phase each change targeted. + +## Versus reference (single-thread field, min of 10) + +| case | bioimage-cpp | reference | speedup | +|---|---|---|---| +| mask2d/field (1024²) | 210.5 ms | scikit-fmm 331.8 ms | **1.58×** | +| mask3d/field (64·128²) | 368.9 ms | scikit-fmm 761.4 ms | **2.06×** | +| mesh/field (V=20000) | 8.4 ms | pygeodesic 1745 ms | 208× (apples-to-oranges: pygeodesic is *exact* MMP, ours is first-order) | + +The mask solver went from parity with scikit-fmm's compiled C (pre-work ~1.0–1.1×) +to a clear lead (1.58× in 2D, 2.06× in 3D). + +## Correctness (all green on the shipping build) + +- `pytest tests/distance/ -q` — 64 passed. The tight external cross-checks ran + (not skipped): `test_mask_matches_scikit_fmm` (residual < 1e-6), + `test_mesh_matches_pygeodesic` (mean rel < 0.06, p95 < 0.12). +- `python check_geodesic_distance.py` — all 6 cases OK, exit 0 (mask residuals + ~1e-12 vs scikit-fmm; mesh rel mean 0.027). +- **Bitwise equivalence**: 13 solver outputs (2D/3D field, anisotropic `sampling`, + `speed`, gradient, mask + mesh pairwise, mesh `speed`) are exactly equal + (`np.array_equal`) between the baseline and optimized builds. Determinism verified + (repeat solves identical). This is the dev-time check behind the "bitwise + identical" claim; the pytest suite guards correctness against the external + references on every run. + +The profiler instrumentation (`BIOIMAGE_PROFILE_SCOPE` "pop"/"relax") is left in both +solvers; it is a no-op unless built with `-C cmake.define.BIOIMAGE_PROFILE=ON`. diff --git a/development/distance/_geodesic_reference.py b/development/distance/_geodesic_reference.py new file mode 100644 index 0000000..a539d69 --- /dev/null +++ b/development/distance/_geodesic_reference.py @@ -0,0 +1,155 @@ +"""Reference geodesic-distance implementations for development cross-checks. + +These are the oracles that ``bic.distance`` geodesic functions are validated +against. There is one backend per geometry, because scikit-fmm only supports +regular Cartesian grids: + +- **masks** -> scikit-fmm (``skfmm``): fast marching / Eikonal solver. +- **meshes** -> pygeodesic: exact Mitchell-Mount-Papadimitriou (MMP) geodesics. + +Not part of the pytest suite; requires ``scikit-fmm`` and/or ``pygeodesic``. +Importing this module is cheap (numpy only); each function imports its backend +lazily and raises a clear ``ImportError`` if the backend is missing. + +Note on the scikit-fmm point-source idiom: to get the distance to a set of +seed voxels, we set ``phi = +1`` everywhere and ``phi = -1`` at the seeds and +take ``|skfmm.distance(phi)|``. The zero contour then sits roughly half a cell +away from each seed, so the reference is offset from a "distance is exactly 0 at +the seed" convention by ~0.5 * dx near the sources. Compare with a tolerance, +or exclude the immediate neighbourhood of the seeds, when checking against an +implementation that initialises the seeds to exactly 0. +""" + +from __future__ import annotations + +import numpy as np + + +# --------------------------------------------------------------------------- # +# masks: scikit-fmm +# --------------------------------------------------------------------------- # + + +def _skfmm(): + try: + import skfmm + except ImportError as error: # pragma: no cover - dev script + raise ImportError( + "scikit-fmm is required for the mask geodesic reference " + "(`pip install scikit-fmm`)" + ) from error + return skfmm + + +def _normalize_dx(sampling, ndim): + if sampling is None: + return 1.0 + if np.isscalar(sampling): + return float(sampling) + dx = [float(s) for s in sampling] + if len(dx) != ndim: + raise ValueError(f"sampling must have length {ndim}, got {len(dx)}") + return dx + + +def reference_geodesic_field_mask(mask, sources, sampling=None, speed=None): + """Geodesic distance field within a mask from a set of source coordinates. + + Mirrors ``bic.distance.geodesic_distance_field``: for every voxel, the + shortest-path distance to the nearest source, constrained to the nonzero + region of ``mask``. Voxels outside the mask (and voxels unreachable from any + source) are returned as ``+inf``. + """ + skfmm = _skfmm() + mask = np.ascontiguousarray(mask) != 0 + sources = np.atleast_2d(np.asarray(sources, dtype=np.int64)) + dx = _normalize_dx(sampling, mask.ndim) + + phi = np.ones(mask.shape, dtype=np.float64) + phi[tuple(sources.T)] = -1.0 + phi = np.ma.MaskedArray(phi, mask=~mask) + + # order=1 matches bioimage-cpp's first-order Godunov scheme (so agreement is + # limited only by the ~0.5-cell seed offset documented above, not by the + # stencil order). + if speed is None: + field = np.ma.abs(skfmm.distance(phi, dx=dx, order=1)) + else: + speed = np.ascontiguousarray(speed, dtype=np.float64) + field = skfmm.travel_time(phi, speed, dx=dx, order=1) + + out = np.full(mask.shape, np.inf, dtype=np.float64) + field_data = np.ma.getdata(field) + valid = ~np.ma.getmaskarray(field) + out[valid] = field_data[valid] + return out + + +def reference_geodesic_distances_mask(mask, points, sampling=None, speed=None): + """Full pairwise geodesic distance matrix between points within a mask. + + Runs the field solve once per point and reads the field at the other + points. The result is symmetrized (per-source solves can differ by a tiny + numerical amount) with a zero diagonal. + """ + points = np.atleast_2d(np.asarray(points, dtype=np.int64)) + n = len(points) + out = np.full((n, n), np.inf, dtype=np.float64) + index = tuple(points.T) + for i in range(n): + field = reference_geodesic_field_mask(mask, points[i][None, :], sampling, speed) + out[i] = field[index] + out = 0.5 * (out + out.T) + np.fill_diagonal(out, 0.0) + return out + + +# --------------------------------------------------------------------------- # +# meshes: pygeodesic (exact MMP) +# --------------------------------------------------------------------------- # + + +def _pygeodesic(): + try: + import pygeodesic.geodesic as geodesic + except ImportError as error: # pragma: no cover - dev script + raise ImportError( + "pygeodesic is required for the mesh geodesic reference " + "(`pip install pygeodesic`)" + ) from error + return geodesic + + +def _mesh_algorithm(vertices, faces): + geodesic = _pygeodesic() + vertices = np.ascontiguousarray(vertices, dtype=np.float64) + faces = np.ascontiguousarray(faces, dtype=np.int32) + return geodesic.PyGeodesicAlgorithmExact(vertices, faces) + + +def reference_geodesic_field_mesh(vertices, faces, sources): + """Geodesic distance field on a triangle mesh from a set of source vertices. + + Mirrors ``bic.distance.geodesic_distance_field_mesh``. Returns a + ``(n_vertices,)`` array of surface geodesic distance to the nearest source. + """ + algo = _mesh_algorithm(vertices, faces) + sources = np.atleast_1d(np.asarray(sources, dtype=np.int32)) + distances, _ = algo.geodesicDistances(sources, None) + return np.asarray(distances, dtype=np.float64) + + +def reference_geodesic_distances_mesh(vertices, faces, points): + """Full pairwise geodesic distance matrix between mesh vertices.""" + algo = _mesh_algorithm(vertices, faces) + points = np.atleast_1d(np.asarray(points, dtype=np.int32)) + n = len(points) + out = np.full((n, n), np.inf, dtype=np.float64) + for i in range(n): + # Read the full single-source field and index the targets; passing an + # explicit target array to pygeodesic can overflow on some meshes. + distances, _ = algo.geodesicDistances(points[i : i + 1], None) + out[i] = np.asarray(distances, dtype=np.float64)[points] + out = 0.5 * (out + out.T) + np.fill_diagonal(out, 0.0) + return out diff --git a/development/distance/benchmark_geodesic.py b/development/distance/benchmark_geodesic.py new file mode 100644 index 0000000..bbcf105 --- /dev/null +++ b/development/distance/benchmark_geodesic.py @@ -0,0 +1,243 @@ +"""Benchmark geodesic distances: bioimage-cpp vs scikit-fmm / pygeodesic. + +Times the ``bic.distance`` geodesic entry points against their reference oracles +(scikit-fmm for masks, pygeodesic exact for meshes) on a range of mask grids and +sphere meshes, reporting per-call statistics and (optionally) the speedup versus +the reference. Reference columns are skipped when the backend is not installed or +when ``--no-ref`` is passed. + +For each case it reports the min, median and spread (p75-p25) of the timed calls +plus ``ns/unit`` — nanoseconds per element processed (voxels for a field solve, +sources x voxels for a pairwise matrix), a size-independent throughput number that +is directly comparable across sizes and dimensionalities and is the right metric +for before/after optimization comparisons. + +For rigorous before/after measurement use a high ``--repeats`` (>=15), pin +``--threads 1`` on the field cases, run on a quiet machine, and dump ``--json`` for +each build so the two runs can be diffed. Interleave baseline/optimized builds to +cancel slow thermal/scheduler drift. + +Not part of the pytest suite; requires scikit-fmm, pygeodesic and scipy (the last +only for the sphere-mesh construction; ``--no-ref`` still needs scipy for meshes). + +Run:: + + python benchmark_geodesic.py --small + python benchmark_geodesic.py --large --threads 0 + python benchmark_geodesic.py --xlarge --only field --threads 1 --repeats 20 \\ + --no-ref --json /tmp/baseline.json +""" + +from __future__ import annotations + +import argparse +import importlib.util +import json +import sys +from statistics import median +from time import perf_counter + +import numpy as np + +import bioimage_cpp as bic + +from _geodesic_reference import ( + reference_geodesic_distances_mask, + reference_geodesic_distances_mesh, + reference_geodesic_field_mask, + reference_geodesic_field_mesh, +) + + +def available(name: str) -> bool: + return importlib.util.find_spec(name) is not None + + +def time_call(fn, repeats: int, warmup: int = 1) -> list[float]: + """Return the list of per-call wall-clock times (seconds), after warmup.""" + for _ in range(warmup): + fn() + timings = [] + for _ in range(repeats): + start = perf_counter() + fn() + timings.append(perf_counter() - start) + return timings + + +def stats(timings: list[float]) -> dict: + ordered = sorted(timings) + n = len(ordered) + p25 = ordered[max(0, (n - 1) // 4)] + p75 = ordered[min(n - 1, (3 * (n - 1)) // 4)] + return { + "min": ordered[0], + "median": median(ordered), + "max": ordered[-1], + "spread": p75 - p25, + } + + +def make_sphere(n_points: int, radius: float = 5.0, seed: int = 0): + from scipy.spatial import ConvexHull + + rng = np.random.default_rng(seed) + pts = rng.standard_normal((n_points, 3)) + pts /= np.linalg.norm(pts, axis=1, keepdims=True) + pts *= radius + verts = np.ascontiguousarray(pts, dtype=np.float64) + faces = np.ascontiguousarray(ConvexHull(pts).simplices, dtype=np.int64) + return verts, faces + + +def scattered_points(n_axis_points: int, ndim: int, extent: int) -> np.ndarray: + coords = np.linspace(1, extent - 2, num=n_axis_points, dtype=np.int64) + grids = np.meshgrid(*([coords] * ndim), indexing="ij") + return np.stack([g.ravel() for g in grids], axis=1).astype(np.int64) + + +def build_cases(args): + """Return a list of (name, bic_fn, ref_fn, n_units) tuples for the size. + + ``n_units`` is the number of elements processed by one call: the voxel/vertex + count for a field solve, or ``n_sources * count`` for a pairwise matrix. It is + the denominator for the ns/unit throughput metric. + """ + if args.small: + mask2d_shape, mask3d_shape, mesh_n = (256, 256), (24, 48, 48), 1500 + elif args.large: + mask2d_shape, mask3d_shape, mesh_n = (1024, 1024), (64, 128, 128), 20000 + elif args.xlarge: + mask2d_shape, mask3d_shape, mesh_n = (1536, 1536), (128, 128, 128), 40000 + else: + mask2d_shape, mask3d_shape, mesh_n = (512, 512), (32, 96, 96), 5000 + + def prod(shape): + n = 1 + for s in shape: + n *= s + return n + + cases = [] + + mask2d = np.ones(mask2d_shape, np.uint8) + src2d = np.array([[mask2d_shape[0] // 2, mask2d_shape[1] // 2]], np.int64) + pts2d = scattered_points(args.n_pairwise, 2, min(mask2d_shape)) + cases.append(( + f"mask2d/field {mask2d_shape}", + lambda: bic.distance.geodesic_distance_field(mask2d, src2d, number_of_threads=args.threads), + lambda: reference_geodesic_field_mask(mask2d, src2d), + prod(mask2d_shape), + )) + cases.append(( + f"mask2d/pairwise N={len(pts2d)}", + lambda: bic.distance.geodesic_distances(mask2d, pts2d, number_of_threads=args.threads), + lambda: reference_geodesic_distances_mask(mask2d, pts2d), + len(pts2d) * prod(mask2d_shape), + )) + + mask3d = np.ones(mask3d_shape, np.uint8) + src3d = np.array([[s // 2 for s in mask3d_shape]], np.int64) + pts3d = scattered_points(args.n_pairwise, 3, min(mask3d_shape)) + cases.append(( + f"mask3d/field {mask3d_shape}", + lambda: bic.distance.geodesic_distance_field(mask3d, src3d, number_of_threads=args.threads), + lambda: reference_geodesic_field_mask(mask3d, src3d), + prod(mask3d_shape), + )) + cases.append(( + f"mask3d/pairwise N={len(pts3d)}", + lambda: bic.distance.geodesic_distances(mask3d, pts3d, number_of_threads=args.threads), + lambda: reference_geodesic_distances_mask(mask3d, pts3d), + len(pts3d) * prod(mask3d_shape), + )) + + verts, faces = make_sphere(mesh_n, seed=args.seed) + src_v = np.array([0], np.int64) + pts_v = np.linspace(0, len(verts) - 1, num=args.n_pairwise, dtype=np.int64) + cases.append(( + f"mesh/field V={len(verts)}", + lambda: bic.distance.geodesic_distance_field_mesh(verts, faces, src_v, number_of_threads=args.threads), + lambda: reference_geodesic_field_mesh(verts, faces, src_v), + len(verts), + )) + cases.append(( + f"mesh/pairwise N={len(pts_v)}", + lambda: bic.distance.geodesic_distances_mesh(verts, faces, pts_v, number_of_threads=args.threads), + lambda: reference_geodesic_distances_mesh(verts, faces, pts_v), + len(pts_v) * len(verts), + )) + + if args.only: + cases = [c for c in cases if args.only in c[0]] + return cases + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--repeats", type=int, default=7) + parser.add_argument("--small", action="store_true") + parser.add_argument("--large", action="store_true") + parser.add_argument("--xlarge", action="store_true") + parser.add_argument("--threads", type=int, default=1, help="0 = hardware concurrency") + parser.add_argument("--n-pairwise", type=int, default=8, dest="n_pairwise", + help="points per axis (mask) / total points (mesh) for pairwise") + parser.add_argument("--seed", type=int, default=0) + parser.add_argument("--no-ref", action="store_true", help="skip reference timings") + parser.add_argument("--only", type=str, default="", + help="substring filter on case name, e.g. 'field' or 'mask3d'") + parser.add_argument("--json", type=str, default="", help="dump results to this JSON path") + args = parser.parse_args() + n_sizes = sum(bool(x) for x in (args.small, args.large, args.xlarge)) + if n_sizes > 1: + print("choose at most one of --small / --large / --xlarge", file=sys.stderr) + return 2 + + have_mask_ref = available("skfmm") and not args.no_ref + have_mesh_ref = available("pygeodesic") and not args.no_ref + if not args.no_ref: + if not available("skfmm"): + print("scikit-fmm not installed: mask reference column skipped", file=sys.stderr) + if not available("pygeodesic"): + print("pygeodesic not installed: mesh reference column skipped", file=sys.stderr) + + cases = build_cases(args) + + header = (f"{'case':>26} {'bic_min':>10} {'bic_med':>10} {'spread':>9} " + f"{'ns/unit':>9} {'ref_med':>10} {'speedup':>9}") + print(f"threads={args.threads} repeats={args.repeats}") + print(header) + print("-" * len(header)) + results = [] + for name, bic_fn, ref_fn, n_units in cases: + bic_t = time_call(bic_fn, args.repeats) + bs = stats(bic_t) + ns_per_unit = bs["median"] / n_units * 1e9 + is_mesh = name.startswith("mesh") + has_ref = have_mesh_ref if is_mesh else have_mask_ref + entry = {"case": name, "n_units": n_units, "threads": args.threads, + "bic_s": bs, "ns_per_unit": ns_per_unit} + if has_ref: + ref_t = time_call(ref_fn, args.repeats) + rs = stats(ref_t) + speed = rs["median"] / bs["median"] + entry["ref_s"] = rs + entry["speedup"] = speed + ref_col, speed_col = f"{rs['median']*1e3:10.2f}", f"{speed:6.2f}x" + else: + ref_col, speed_col = f"{'n/a':>10}", f"{'n/a':>9}" + print(f"{name:>26} {bs['min']*1e3:10.2f} {bs['median']*1e3:10.2f} " + f"{bs['spread']*1e3:9.2f} {ns_per_unit:9.1f} {ref_col} {speed_col:>9}") + results.append(entry) + + if args.json: + with open(args.json, "w") as fh: + json.dump({"threads": args.threads, "repeats": args.repeats, + "results": results}, fh, indent=2) + print(f"wrote {args.json}", file=sys.stderr) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/development/distance/check_geodesic_distance.py b/development/distance/check_geodesic_distance.py new file mode 100644 index 0000000..dc85a84 --- /dev/null +++ b/development/distance/check_geodesic_distance.py @@ -0,0 +1,224 @@ +"""Cross-check bioimage-cpp geodesic distances against scikit-fmm / pygeodesic. + +Builds small mask (grid) and triangle-mesh test geometries, computes the +reference geodesic fields and pairwise matrices with the oracles in +``_geodesic_reference`` (scikit-fmm for masks, pygeodesic for meshes), and +compares them against ``bic.distance``: + + geodesic_distance_field <-> skfmm (mask field) + geodesic_distances <-> skfmm (mask pairwise) + geodesic_distance_field_mesh <-> pygeodesic exact (mesh field) + geodesic_distances_mesh <-> pygeodesic exact (mesh pairwise) + +The bioimage-cpp solvers are not implemented yet, so the bic calls currently +raise ``RuntimeError("... not yet implemented")``. Those cases are reported as +PENDING (not a failure); the script is the acceptance check that turns green +once the algorithms land. Cases whose reference backend is not installed are +reported as NO-REF and also skipped. + +Not part of the pytest suite; requires scikit-fmm, pygeodesic and scipy. +""" + +from __future__ import annotations + +import argparse +import sys + +import numpy as np + +import bioimage_cpp as bic + +from _geodesic_reference import ( + reference_geodesic_distances_mask, + reference_geodesic_distances_mesh, + reference_geodesic_field_mask, + reference_geodesic_field_mesh, +) + + +# --------------------------------------------------------------------------- # +# test geometries +# --------------------------------------------------------------------------- # + + +def make_mask_2d(): + """40x40 domain with a horizontal wall + gap; geodesic paths must detour.""" + mask = np.ones((40, 40), dtype=np.uint8) + mask[20, :34] = 0 # wall across most of the row, gap near the right edge + sources = np.array([[5, 5]], dtype=np.int64) + points = np.array([[5, 5], [35, 5], [35, 35], [5, 35]], dtype=np.int64) + return mask, sources, points + + +def make_mask_3d(): + """Solid ball on a 30^3 grid.""" + n, c, r = 30, 15, 12 + zz, yy, xx = np.ogrid[:n, :n, :n] + mask = ((zz - c) ** 2 + (yy - c) ** 2 + (xx - c) ** 2 <= r * r).astype(np.uint8) + sources = np.array([[c - r + 1, c, c]], dtype=np.int64) + points = np.array( + [[c - r + 1, c, c], [c + r - 1, c, c], [c, c - r + 1, c], [c, c, c + r - 1]], + dtype=np.int64, + ) + return mask, sources, points + + +def make_mesh(n_points=600, radius=5.0, seed=0): + """A closed sphere mesh from the convex hull of points on the sphere. + + scipy's ConvexHull gives a clean manifold triangulation (pygeodesic can + overflow on the degenerate faces marching_cubes produces). + """ + from scipy.spatial import ConvexHull + + rng = np.random.default_rng(seed) + pts = rng.standard_normal((n_points, 3)) + pts /= np.linalg.norm(pts, axis=1, keepdims=True) + pts *= radius + verts = np.ascontiguousarray(pts, dtype=np.float64) + faces = np.ascontiguousarray(ConvexHull(pts).simplices, dtype=np.int64) + n_verts = len(verts) + sources = np.array([0], dtype=np.int64) + points = np.linspace(0, n_verts - 1, num=6, dtype=np.int64) + return verts, faces, sources, points + + +# --------------------------------------------------------------------------- # +# comparison harness +# --------------------------------------------------------------------------- # + + +def compare(case, atol, rtol, mesh_rtol): + """Compare one case; ``case`` is ``(name, kind, reference_fn, bic_fn)``. + + Masks use the same first-order scheme as scikit-fmm, so after removing the + ~0.5-cell level-set seed offset (see ``_geodesic_reference``) the residual is + tiny. Meshes compare our first-order FMM to the exact pygeodesic distances + via relative error (looser: first-order + no obtuse unfolding yet). + """ + name, kind, reference_fn, bic_fn = case + try: + ref = np.asarray(reference_fn(), dtype=np.float64) + except ImportError as error: # backend missing + return {"name": name, "status": "NO-REF", "metric": str(error).split("(")[0].strip()} + + row = {"name": name, "ref_shape": ref.shape} + finite = np.isfinite(ref) + if finite.any(): + row["ref_range"] = (float(ref[finite].min()), float(ref[finite].max())) + + try: + got = np.asarray(bic_fn(), dtype=np.float64) + except RuntimeError as error: + if "not yet implemented" in str(error): + row["status"] = "PENDING" + return row + raise + + inf_match = np.array_equal(~finite, ~np.isfinite(got)) + both = finite & np.isfinite(got) + if not both.any(): + row["status"] = "OK" if inf_match else "FAIL" + row["metric"] = "no finite overlap" + return row + + if kind == "mask": + away = both & (ref > 1.0) # skip the seed neighbourhood / matrix diagonal + sel = away if away.any() else both + diff = got[sel] - ref[sel] + offset = float(np.median(diff)) + max_resid = float(np.abs(diff - offset).max()) + scale = float(ref[sel].max()) + ok = inf_match and max_resid < (rtol * scale + atol) + row["metric"] = f"offset={offset:+.3f} resid_max={max_resid:.2e}" + else: # mesh, vs exact pygeodesic + reach = both & (ref > 1e-9) # exclude the sources (ref == 0) + rel = np.abs(got[reach] - ref[reach]) / ref[reach] + mean_rel, p95, mx = float(rel.mean()), float(np.percentile(rel, 95)), float(rel.max()) + ok = inf_match and mean_rel < 0.06 and p95 < mesh_rtol + row["metric"] = f"rel mean={mean_rel:.3f} p95={p95:.3f} max={mx:.3f}" + + row["status"] = "OK" if ok else "FAIL" + return row + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--seed", type=int, default=0) + parser.add_argument("--repeats", type=int, default=1) + parser.add_argument("--atol", type=float, default=1e-2, help="mask absolute residual floor") + parser.add_argument("--rtol", type=float, default=2e-2, help="mask residual vs max distance") + parser.add_argument("--mesh-rtol", type=float, default=0.12, help="mesh p95 relative error") + args = parser.parse_args() + + mask2d, src2d, pts2d = make_mask_2d() + mask3d, src3d, pts3d = make_mask_3d() + try: + verts, faces, src_v, pts_v = make_mesh() + mesh_ok = True + except ImportError as error: + sys.stderr.write(f"scipy not installed, skipping mesh cases: {error}\n") + mesh_ok = False + + cases = [ + ( + "mask2d/field", "mask", + lambda: reference_geodesic_field_mask(mask2d, src2d), + lambda: bic.distance.geodesic_distance_field(mask2d, src2d), + ), + ( + "mask2d/pairwise", "mask", + lambda: reference_geodesic_distances_mask(mask2d, pts2d), + lambda: bic.distance.geodesic_distances(mask2d, pts2d), + ), + ( + "mask3d/field", "mask", + lambda: reference_geodesic_field_mask(mask3d, src3d), + lambda: bic.distance.geodesic_distance_field(mask3d, src3d), + ), + ( + "mask3d/pairwise", "mask", + lambda: reference_geodesic_distances_mask(mask3d, pts3d), + lambda: bic.distance.geodesic_distances(mask3d, pts3d), + ), + ] + if mesh_ok: + cases += [ + ( + "mesh/field", "mesh", + lambda: reference_geodesic_field_mesh(verts, faces, src_v), + lambda: bic.distance.geodesic_distance_field_mesh(verts, faces, src_v), + ), + ( + "mesh/pairwise", "mesh", + lambda: reference_geodesic_distances_mesh(verts, faces, pts_v), + lambda: bic.distance.geodesic_distances_mesh(verts, faces, pts_v), + ), + ] + + header = f"{'case':>16} {'status':>8} {'ref_shape':>14} {'ref_range':>20} {'agreement':>34}" + print(header) + print("-" * len(header)) + any_fail = False + for case in cases: + r = compare(case, args.atol, args.rtol, args.mesh_rtol) + rng = r.get("ref_range") + rng_s = f"[{rng[0]:.3f}, {rng[1]:.3f}]" if rng else "-" + shape_s = str(r.get("ref_shape", "-")) + metric_s = r.get("metric", "-") + print(f"{r['name']:>16} {r['status']:>8} {shape_s:>14} {rng_s:>20} {metric_s:>34}") + if r["status"] == "FAIL": + any_fail = True + + print() + if any_fail: + print("FAIL: at least one case disagrees with the reference.", file=sys.stderr) + sys.exit(1) + print( + "No mismatches. (PENDING = bic solver not implemented yet; " + "NO-REF = reference backend not installed.)" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/distance/complex-shape.tif b/examples/distance/complex-shape.tif new file mode 100644 index 0000000..cd9350e Binary files /dev/null and b/examples/distance/complex-shape.tif differ diff --git a/examples/distance/geodesic_2d.py b/examples/distance/geodesic_2d.py new file mode 100644 index 0000000..33b579b --- /dev/null +++ b/examples/distance/geodesic_2d.py @@ -0,0 +1,21 @@ +import imageio.v3 as imageio +import napari +import numpy as np + +from bioimage_cpp.distance import geodesic_distance_field, distance_transform + +mask = (imageio.imread("./complex-shape.tif") == 2).astype("uint8") +# compute the most central point --- distance transform maximum +max_point = np.argmax(distance_transform(mask).ravel()) +max_point = np.unravel_index(max_point, mask.shape) + +sources = np.array([max_point]) +distance_field, gradients = geodesic_distance_field(mask, sources, return_gradient=True) +gradients = gradients.transpose((2, 0, 1)) + +v = napari.Viewer() +v.add_image(distance_field) +v.add_image(gradients) +v.add_labels(mask) +v.add_points([max_point]) +napari.run() diff --git a/examples/distance/geodesic_3d.py b/examples/distance/geodesic_3d.py new file mode 100644 index 0000000..8154590 --- /dev/null +++ b/examples/distance/geodesic_3d.py @@ -0,0 +1,20 @@ +import napari +import numpy as np + +from bioimage_cpp.distance import geodesic_distance_field, distance_transform +from shape import make_shape + +mask = make_shape() +max_point = np.argmax(distance_transform(mask).ravel()) +max_point = np.unravel_index(max_point, mask.shape) + +sources = np.array([max_point]) +distance_field, gradients = geodesic_distance_field(mask, sources, return_gradient=True) +gradients = gradients.transpose((3, 0, 1, 2)) + +v = napari.Viewer() +v.add_image(distance_field) +v.add_image(gradients) +v.add_labels(mask) +v.add_points([max_point]) +napari.run() diff --git a/examples/distance/geodesic_mesh.py b/examples/distance/geodesic_mesh.py new file mode 100644 index 0000000..ea819bf --- /dev/null +++ b/examples/distance/geodesic_mesh.py @@ -0,0 +1,36 @@ +import napari +import numpy as np +from scipy.spatial import cKDTree +from skimage.measure import marching_cubes + +from bioimage_cpp.distance import geodesic_distance_field_mesh +from shape import make_shape + +# Create the mesh based on our sample shape. +mask = make_shape() +verts, faces, _, _ = marching_cubes(mask, level=0.5) + +# Create some points in the volume as sources for the distances. +c = mask.shape[0] / 2.0 +angles = np.linspace(0.0, 2.0 * np.pi, num=5, endpoint=False) +points = np.array( + [[c, c + 33.0 * np.sin(a), c + 33.0 * np.cos(a)] for a in angles] +) + +# Snap the points to the closest mesh vertex. +_, idx = cKDTree(verts).query(points) +sources = np.unique(idx) + +# Compute the distance field on the mesh. +field = geodesic_distance_field_mesh(verts, faces, sources) + +# Alternatively, you can also compute the pairwise geodesic distance +# between all the points like this: +# pairwise_distances = geodesic_distances_mesh(verts, faces, sources, number_of_threads=1) + +# Display the mesh and points in napari. +# The mesh gets the distance field as values. +v = napari.Viewer() +v.add_surface((verts, faces, field)) +v.add_points(points) +napari.run() diff --git a/examples/distance/shape.py b/examples/distance/shape.py new file mode 100644 index 0000000..3991958 --- /dev/null +++ b/examples/distance/shape.py @@ -0,0 +1,10 @@ +import numpy as np + + +def make_shape(size=80): + R, r = 24.0, 9.0 + c = size / 2.0 + zz, yy, xx = np.ogrid[:size, :size, :size] + q = np.sqrt((yy - c) ** 2 + (xx - c) ** 2) - R + mask = (q ** 2 + (zz - c) ** 2 <= r ** 2).astype("uint8") + return mask diff --git a/include/bioimage_cpp/detail/grid.hxx b/include/bioimage_cpp/detail/grid.hxx index d6afc76..ce89e9c 100644 --- a/include/bioimage_cpp/detail/grid.hxx +++ b/include/bioimage_cpp/detail/grid.hxx @@ -28,6 +28,26 @@ inline std::vector c_order_strides(const std::vector &strides, + std::size_t ndim, + std::ptrdiff_t *coords_out +) { + for (std::size_t axis = 0; axis < ndim; ++axis) { + const auto stride = static_cast(strides[axis]); + const auto coord = node / stride; + coords_out[axis] = static_cast(coord); + node -= coord * stride; + } +} + // Translate a flat node index by a per-axis offset on a row-major grid. // // Returns true when the offset keeps the result inside the grid, in which case diff --git a/include/bioimage_cpp/distance/detail/fast_marching.hxx b/include/bioimage_cpp/distance/detail/fast_marching.hxx new file mode 100644 index 0000000..fd86eaf --- /dev/null +++ b/include/bioimage_cpp/distance/detail/fast_marching.hxx @@ -0,0 +1,264 @@ +#pragma once + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/indexed_heap.hxx" +#include "bioimage_cpp/detail/profile.hxx" + +#include +#include +#include +#include +#include +#include +#include + +// First-order Godunov fast marching method (Eikonal solver) on a regular grid. +// This is the mask backend for the geodesic-distance functions; it matches the +// scheme used by scikit-fmm (order=1). See geodesic_mask.hxx for the public API. + +namespace bioimage_cpp::distance::detail { + +inline constexpr double kFmmInfinity = std::numeric_limits::infinity(); + +enum class FmmState : std::uint8_t { Far = 0, Narrow = 1, Frozen = 2 }; + +// Reusable workspace for repeated single-/multi-source solves on the same grid. +// `dist`/`state` and the narrow-band heap are allocated once and reset between +// solves so the pairwise driver can reuse one workspace per thread. +class GridFastMarching { +public: + GridFastMarching( + const ConstArrayView &mask, + const std::vector &sampling, + const ConstArrayView *speed + ) + : shape_(mask.shape), + strides_(bioimage_cpp::detail::c_order_strides(mask.shape)), + n_(bioimage_cpp::detail::number_of_elements(mask.shape)), + ndim_(mask.shape.size()), + mask_(mask.data), + speed_(speed != nullptr ? speed->data : nullptr), + dist_(n_, kFmmInfinity), + state_(n_, FmmState::Far), + heap_(n_) { + h_ = sampling; + inv_h2_.resize(ndim_); + for (std::size_t d = 0; d < ndim_; ++d) { + const double h = sampling[d]; + inv_h2_[d] = 1.0 / (h * h); + } + // Axis neighbours are exactly +-1 along one axis, so the flat-index + // delta to a neighbour is +-strides_[d] (a single add), and its bounds + // reduce to a per-axis compare against the popped voxel's coordinate. + // We therefore decode the popped voxel once (coords_) instead of calling + // valid_offset_target (div + mod per axis) ~2*ndim + (2*ndim)^2 times. + coords_.resize(ndim_); + nb_coords_.resize(ndim_); + terms_.reserve(ndim_); + } + + [[nodiscard]] std::size_t size() const { return n_; } + [[nodiscard]] double distance(std::size_t index) const { return dist_[index]; } + [[nodiscard]] const std::vector &distances() const { return dist_; } + + // Run the marching front from the given in-mask source voxels (linear + // indices). Sources outside the mask are ignored. Results are left in + // `dist_`; unreachable / out-of-domain voxels keep +inf. + void solve(const std::vector &sources) { + reset(); + BIOIMAGE_PROFILE_INIT(prof); + for (const auto source : sources) { + if (mask_[source] == 0) { + continue; // a source outside the domain cannot seed the front + } + dist_[source] = 0.0; + state_[source] = FmmState::Narrow; + heap_.push_or_change(source, 0.0); + } + + while (!heap_.empty()) { + std::size_t p = 0; + { + BIOIMAGE_PROFILE_SCOPE(prof, "pop"); + const auto top = heap_.pop(); + p = top.key; + } + state_[p] = FmmState::Frozen; + + BIOIMAGE_PROFILE_SCOPE(prof, "relax"); + // Decode p once; a neighbour's coords differ from p's in one axis + // only, so nb_coords_ is p's coords with that single entry adjusted. + bioimage_cpp::detail::coords_from_index(p, strides_, ndim_, coords_.data()); + for (std::size_t d = 0; d < ndim_; ++d) { + nb_coords_[d] = coords_[d]; + } + // Visit +axis then -axis for each axis in turn: the same neighbour + // order (and thus heap-insertion order) as the previous offsets_ loop. + for (std::size_t d = 0; d < ndim_; ++d) { + const std::ptrdiff_t cd = coords_[d]; + for (int side = 0; side < 2; ++side) { + const std::ptrdiff_t ncd = (side == 0) ? cd + 1 : cd - 1; + if (ncd < 0 || ncd >= shape_[d]) { + continue; + } + const std::uint64_t nb = (side == 0) + ? p + static_cast(strides_[d]) + : p - static_cast(strides_[d]); + if (mask_[nb] == 0 || state_[nb] == FmmState::Frozen) { + continue; + } + nb_coords_[d] = ncd; + const double t = solve_eikonal(nb, nb_coords_.data()); + nb_coords_[d] = cd; + if (t < dist_[nb]) { + dist_[nb] = t; + state_[nb] = FmmState::Narrow; + heap_.push_or_change(nb, t); + } + } + } + } + BIOIMAGE_PROFILE_REPORT(prof); + } + + // Write the first-order upwind gradient of the solved field into `grad` + // (shape (*shape, ndim), row-major float32). Component i is dT/dx_i, pointing + // toward increasing distance (away from the nearest source); it is 0 at + // sources / local minima, background, and unreachable voxels. Following + // -grad traces the geodesic back toward the source; ||grad|| ~= slowness. + void write_gradient(ArrayView &grad) const { + std::vector coords(ndim_); + for (std::size_t p = 0; p < n_; ++p) { + float *out = grad.data + p * ndim_; + const double tp = dist_[p]; + if (mask_[p] == 0 || !std::isfinite(tp)) { + for (std::size_t d = 0; d < ndim_; ++d) { + out[d] = 0.0f; + } + continue; + } + bioimage_cpp::detail::coords_from_index(p, strides_, ndim_, coords.data()); + for (std::size_t d = 0; d < ndim_; ++d) { + // +strides_[d] is the plus side, -strides_[d] the minus side. + // Keep only strictly-upwind (smaller-valued, in-mask) neighbours. + double t_plus = kFmmInfinity; + double t_minus = kFmmInfinity; + if (coords[d] + 1 < shape_[d]) { + const std::uint64_t nb = p + static_cast(strides_[d]); + if (mask_[nb] != 0 && dist_[nb] < tp) { + t_plus = dist_[nb]; + } + } + if (coords[d] > 0) { + const std::uint64_t nb = p - static_cast(strides_[d]); + if (mask_[nb] != 0 && dist_[nb] < tp) { + t_minus = dist_[nb]; + } + } + + double g = 0.0; + if (t_minus == kFmmInfinity && t_plus == kFmmInfinity) { + g = 0.0; // source / local minimum: no upwind neighbour + } else if (t_minus <= t_plus) { + g = (tp - t_minus) / h_[d]; // front from -i => dT/dx_i >= 0 + } else { + g = (t_plus - tp) / h_[d]; // front from +i => dT/dx_i <= 0 + } + out[d] = static_cast(g); + } + } + } + +private: + struct Term { + double u; // smaller frozen neighbour value along an axis + double inv_h2; // 1 / spacing^2 for that axis + }; + + void reset() { + std::fill(dist_.begin(), dist_.end(), kFmmInfinity); + std::fill(state_.begin(), state_.end(), FmmState::Far); + heap_.clear(); + } + + // First-order Godunov update at voxel p (assumed in-mask, not frozen). + // `coords` are p's per-axis coordinates (p differs from the popped voxel in + // one axis only, so the caller supplies them without a fresh decode). + double solve_eikonal(std::uint64_t p, const std::ptrdiff_t *coords) { + const double f = (speed_ != nullptr) ? speed_[p] : 1.0; + if (!(f > 0.0)) { + return kFmmInfinity; // zero/negative speed is impassable + } + const double slowness = 1.0 / f; + + terms_.clear(); + for (std::size_t d = 0; d < ndim_; ++d) { + double best = kFmmInfinity; + const std::ptrdiff_t cd = coords[d]; + if (cd + 1 < shape_[d]) { + const std::uint64_t nb = p + static_cast(strides_[d]); + if (state_[nb] == FmmState::Frozen) { + best = std::min(best, dist_[nb]); + } + } + if (cd > 0) { + const std::uint64_t nb = p - static_cast(strides_[d]); + if (state_[nb] == FmmState::Frozen) { + best = std::min(best, dist_[nb]); + } + } + if (best != kFmmInfinity) { + terms_.push_back({best, inv_h2_[d]}); + } + } + if (terms_.empty()) { + return kFmmInfinity; + } + + std::sort(terms_.begin(), terms_.end(), [](const Term &a, const Term &b) { + return a.u < b.u; + }); + + // Incrementally include axis terms (ascending u) and solve the quadratic + // sum_w t^2 - 2 sum_wu t + (sum_wu2 - slowness^2) = 0 + // stopping once the candidate no longer exceeds the next unused u. + const double s2 = slowness * slowness; + double sum_w = 0.0, sum_wu = 0.0, sum_wu2 = 0.0; + double t = kFmmInfinity; + for (std::size_t k = 0; k < terms_.size(); ++k) { + sum_w += terms_[k].inv_h2; + sum_wu += terms_[k].inv_h2 * terms_[k].u; + sum_wu2 += terms_[k].inv_h2 * terms_[k].u * terms_[k].u; + + const double disc = sum_wu * sum_wu - sum_w * (sum_wu2 - s2); + if (disc < 0.0) { + break; // keep the candidate from fewer terms + } + const double candidate = (sum_wu + std::sqrt(disc)) / sum_w; + t = candidate; + if (k + 1 == terms_.size() || candidate <= terms_[k + 1].u) { + break; + } + } + return t; + } + + std::vector shape_; + std::vector strides_; + std::size_t n_; + std::size_t ndim_; + const std::uint8_t *mask_; + const double *speed_; + std::vector inv_h2_; + std::vector h_; + std::vector coords_; // scratch: popped voxel's coords + std::vector nb_coords_; // scratch: current neighbour's coords + + std::vector dist_; + std::vector state_; + bioimage_cpp::detail::DenseIndexedHeap> heap_; + std::vector terms_; +}; + +} // namespace bioimage_cpp::distance::detail diff --git a/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx b/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx new file mode 100644 index 0000000..5262611 --- /dev/null +++ b/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx @@ -0,0 +1,246 @@ +#pragma once + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/indexed_heap.hxx" +#include "bioimage_cpp/detail/profile.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// First-order fast marching method on a triangle mesh (Kimmel & Sethian, PNAS +// 1998): the mesh backend for the geodesic-distance functions. See +// geodesic_mesh.hxx for the public API and development/distance/ for the exact +// (pygeodesic) reference oracle. +// +// The per-triangle update solves for a linear arrival-time field with gradient +// magnitude equal to the local slowness, and accepts it only when the +// characteristic falls inside the triangle wedge (otherwise it falls back to +// the two edge updates). Obtuse-angle "unfolding" is not done yet, so accuracy +// is first-order and degrades on very obtuse triangulations. + +namespace bioimage_cpp::distance::detail { + +inline constexpr double kMeshFmmInfinity = std::numeric_limits::infinity(); + +class MeshFastMarching { +public: + MeshFastMarching( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const ConstArrayView *speed + ) + : vertices_(vertices.data), + n_vertices_(static_cast(vertices.shape[0])), + faces_(faces.data), + n_faces_(static_cast(faces.shape[0])), + speed_(speed != nullptr ? speed->data : nullptr), + dist_(n_vertices_, kMeshFmmInfinity), + state_(n_vertices_, MeshState::Far), + heap_(n_vertices_) { + build_vertex_face_incidence(); + } + + [[nodiscard]] std::size_t size() const { return n_vertices_; } + [[nodiscard]] double distance(std::size_t vertex) const { return dist_[vertex]; } + [[nodiscard]] const std::vector &distances() const { return dist_; } + + // Run the marching front from the given source vertex indices. Results are + // left in `dist_`; vertices in a different connected component keep +inf. + void solve(const std::vector &sources) { + reset(); + BIOIMAGE_PROFILE_INIT(prof); + for (const auto source : sources) { + dist_[source] = 0.0; + state_[source] = MeshState::Narrow; + heap_.push_or_change(source, 0.0); + } + + while (!heap_.empty()) { + std::size_t v = 0; + { + BIOIMAGE_PROFILE_SCOPE(prof, "pop"); + const auto top = heap_.pop(); + v = top.key; + } + state_[v] = MeshState::Frozen; + + BIOIMAGE_PROFILE_SCOPE(prof, "relax"); + // Only the faces incident to the just-frozen v carry new frozen + // information, so relax each of their other corners from that face + // alone (using v, now frozen, as the/one frozen source) rather than + // rescanning every incident face of w. This is exactly equivalent to + // the previous update_vertex rescan: dist_[w] is a running minimum + // and already reflects every other incident face's contribution from + // the earlier freeze of that face's corner, so min(dist_[w], full + // rescan) collapses to min(dist_[w], this-face candidate). + for (std::size_t k = face_offsets_[v]; k < face_offsets_[v + 1]; ++k) { + const std::size_t f = face_ids_[k]; + for (std::size_t corner = 0; corner < 3; ++corner) { + const std::size_t w = static_cast(faces_[f * 3 + corner]); + if (w == v || state_[w] == MeshState::Frozen) { + continue; + } + // The two corners of f other than w, in face-storage order, + // matching the previous update_vertex so triangle_update sees + // identical operands (its acute update is only symmetric up + // to rounding, so operand order is load-bearing for + // bit-identical output). + std::array others{}; + std::size_t count = 0; + for (std::size_t c2 = 0; c2 < 3; ++c2) { + const auto u = static_cast(faces_[f * 3 + c2]); + if (u != w) { + others[count++] = u; + } + } + const std::size_t a = others[0]; + const std::size_t b = others[1]; + const bool fa = state_[a] == MeshState::Frozen; + const bool fb = state_[b] == MeshState::Frozen; + double t = kMeshFmmInfinity; + if (fa && fb) { + t = triangle_update(w, a, b); + } else if (fa) { + t = dist_[a] + slowness(w) * edge_length(w, a); + } else if (fb) { + t = dist_[b] + slowness(w) * edge_length(w, b); + } + if (t < dist_[w]) { + dist_[w] = t; + state_[w] = MeshState::Narrow; + heap_.push_or_change(w, t); + } + } + } + } + BIOIMAGE_PROFILE_REPORT(prof); + } + +private: + enum class MeshState : std::uint8_t { Far = 0, Narrow = 1, Frozen = 2 }; + + void reset() { + std::fill(dist_.begin(), dist_.end(), kMeshFmmInfinity); + std::fill(state_.begin(), state_.end(), MeshState::Far); + heap_.clear(); + } + + // Build the vertex -> incident-faces CSR (offsets + face ids), mirroring the + // count-then-scatter pattern in detail::mesh_smoothing::build_adjacency. + void build_vertex_face_incidence() { + face_offsets_.assign(n_vertices_ + 1, 0); + for (std::size_t f = 0; f < n_faces_; ++f) { + for (std::size_t corner = 0; corner < 3; ++corner) { + const std::int64_t v = faces_[f * 3 + corner]; + if (v < 0 || static_cast(v) >= n_vertices_) { + throw std::invalid_argument( + "faces contains vertex index " + std::to_string(v) + + " out of range [0, " + std::to_string(n_vertices_) + ")" + ); + } + ++face_offsets_[static_cast(v) + 1]; + } + } + for (std::size_t i = 1; i < face_offsets_.size(); ++i) { + face_offsets_[i] += face_offsets_[i - 1]; + } + face_ids_.resize(face_offsets_.back()); + std::vector insert_pos(face_offsets_.begin(), face_offsets_.end() - 1); + for (std::size_t f = 0; f < n_faces_; ++f) { + for (std::size_t corner = 0; corner < 3; ++corner) { + const auto v = static_cast(faces_[f * 3 + corner]); + face_ids_[insert_pos[v]++] = f; + } + } + } + + [[nodiscard]] double slowness(std::size_t vertex) const { + if (speed_ == nullptr) { + return 1.0; + } + const double f = speed_[vertex]; + return (f > 0.0) ? 1.0 / f : kMeshFmmInfinity; + } + + [[nodiscard]] double edge_length(std::size_t u, std::size_t v) const { + const double *pu = vertices_ + u * 3; + const double *pv = vertices_ + v * 3; + const double dx = pu[0] - pv[0]; + const double dy = pu[1] - pv[1]; + const double dz = pu[2] - pv[2]; + return std::sqrt(dx * dx + dy * dy + dz * dz); + } + + // Kimmel-Sethian planar-wavefront update of vertex c from the triangle + // (c, a, b) with a, b frozen. Returns the accepted arrival time, or the + // better of the two edge updates when the acute update is not causal. + double triangle_update(std::size_t c, std::size_t a, std::size_t b) { + const double s = slowness(c); + const double len_a = edge_length(c, b); // |CB| pairs with T_B + const double len_b = edge_length(c, a); // |CA| pairs with T_A + const double edge_cand = + std::min(dist_[a] + s * len_b, dist_[b] + s * len_a); + if (!(len_a > 0.0) || !(len_b > 0.0) || std::isinf(s)) { + return edge_cand; + } + + const double *pc = vertices_ + c * 3; + const double *pa = vertices_ + a * 3; + const double *pb = vertices_ + b * 3; + double dot = 0.0; + for (std::size_t d = 0; d < 3; ++d) { + dot += (pa[d] - pc[d]) * (pb[d] - pc[d]); + } + const double cos_theta = dot / (len_a * len_b); + const double sin2 = std::max(0.0, 1.0 - cos_theta * cos_theta); + + const double ta = dist_[a]; + const double tb = dist_[b]; + const double aa = len_a; // |CB| + const double bb = len_b; // |CA| + + // F t^2 - 2 K t + L = 0 for the linear arrival-time field over the + // triangle with |grad T| = s (see header note for the derivation). + const double F = aa * aa + bb * bb - 2.0 * aa * bb * cos_theta; + const double K = aa * aa * ta - aa * bb * cos_theta * (ta + tb) + bb * bb * tb; + const double L = aa * aa * ta * ta - 2.0 * aa * bb * cos_theta * ta * tb + + bb * bb * tb * tb - s * s * aa * aa * bb * bb * sin2; + const double disc = K * K - F * L; + if (F > 0.0 && disc >= 0.0) { + const double t = (K + std::sqrt(disc)) / F; + if (t >= ta && t >= tb) { + // Causality: the characteristic -grad T must point into the + // wedge, i.e. both barycentric weights are non-negative. + const double lam = aa * aa * (t - ta) - aa * bb * cos_theta * (t - tb); + const double mu = bb * bb * (t - tb) - aa * bb * cos_theta * (t - ta); + if (lam >= 0.0 && mu >= 0.0) { + return std::min(edge_cand, t); + } + } + } + return edge_cand; + } + + const double *vertices_; + std::size_t n_vertices_; + const std::int64_t *faces_; + std::size_t n_faces_; + const double *speed_; + + std::vector face_offsets_; + std::vector face_ids_; + + std::vector dist_; + std::vector state_; + bioimage_cpp::detail::DenseIndexedHeap> heap_; +}; + +} // namespace bioimage_cpp::distance::detail diff --git a/include/bioimage_cpp/distance_transform.hxx b/include/bioimage_cpp/distance/distance_transform.hxx similarity index 100% rename from include/bioimage_cpp/distance_transform.hxx rename to include/bioimage_cpp/distance/distance_transform.hxx diff --git a/include/bioimage_cpp/distance/geodesic_mask.hxx b/include/bioimage_cpp/distance/geodesic_mask.hxx new file mode 100644 index 0000000..d6a8998 --- /dev/null +++ b/include/bioimage_cpp/distance/geodesic_mask.hxx @@ -0,0 +1,140 @@ +#pragma once + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/threading.hxx" +#include "bioimage_cpp/distance/detail/fast_marching.hxx" + +#include +#include +#include +#include +#include + +// Geodesic distances within a mask (regular Cartesian grid). +// +// A "mask" defines the geometry: distances are shortest-path lengths that stay +// inside the nonzero (foreground) region and never cross background voxels. +// Solved with the first-order Godunov fast marching method (see +// detail/fast_marching.hxx), matching the scikit-fmm scheme. The reference +// oracle lives in development/distance/. + +namespace bioimage_cpp::distance { + +namespace detail { + +// Convert an (m, ndim) array of int64 voxel coordinates (NumPy axis order) to +// flat row-major indices, validating that every coordinate is in bounds. +inline std::vector linear_indices_from_coords( + const ConstArrayView &coords, + const std::vector &shape, + const std::vector &strides, + const char *name +) { + const auto m = static_cast(coords.shape[0]); + const auto ndim = shape.size(); + std::vector indices(m); + for (std::size_t i = 0; i < m; ++i) { + std::size_t linear = 0; + for (std::size_t d = 0; d < ndim; ++d) { + const std::int64_t c = coords.data[i * ndim + d]; + if (c < 0 || c >= shape[d]) { + throw std::invalid_argument( + std::string(name) + " coordinate out of bounds: " + name + "[" + + std::to_string(i) + ", " + std::to_string(d) + "]=" + + std::to_string(c) + " not in [0, " + std::to_string(shape[d]) + ")" + ); + } + linear += static_cast(c) * static_cast(strides[d]); + } + indices[i] = linear; + } + return indices; +} + +} // namespace detail + +// Geodesic distance field within a mask from a set of source coordinates. +// +// mask nonzero = inside domain; distances propagate only through nonzero voxels. +// sources: (n_sources, ndim) int64 voxel coords (NumPy axis order). A single +// source is a 1-row array. sampling: per-axis spacing (length ndim). speed: +// optional, same shape as mask (nullptr => unit speed). distances (out): +// float64, mask.shape; +inf where unreachable / outside the domain. gradient: +// optional (nullptr to skip); when given, the first-order upwind gradient of +// the field is written into it as float32 with shape (*mask.shape, ndim) — see +// GridFastMarching::write_gradient. +inline void geodesic_distance_field( + const ConstArrayView &mask, + const ConstArrayView &sources, + const std::vector &sampling, + const ConstArrayView *speed, + ArrayView &distances, + std::size_t /*n_threads*/ = 1, + ArrayView *gradient = nullptr +) { + const auto strides = bioimage_cpp::detail::c_order_strides(mask.shape); + const auto source_indices = + detail::linear_indices_from_coords(sources, mask.shape, strides, "sources"); + + detail::GridFastMarching solver(mask, sampling, speed); + solver.solve(source_indices); + + const auto &dist = solver.distances(); + std::copy(dist.begin(), dist.end(), distances.data); + + if (gradient != nullptr) { + solver.write_gradient(*gradient); + } +} + +// Full pairwise geodesic distance matrix between points within a mask. +// +// points: (n_points, ndim) int64. distances (out): (n_points, n_points) float64, +// symmetric with a zero diagonal; +inf where two points are not connected inside +// the domain. +inline void geodesic_distances( + const ConstArrayView &mask, + const ConstArrayView &points, + const std::vector &sampling, + const ConstArrayView *speed, + ArrayView &distances, + std::size_t n_threads = 1 +) { + const auto strides = bioimage_cpp::detail::c_order_strides(mask.shape); + const auto point_indices = + detail::linear_indices_from_coords(points, mask.shape, strides, "points"); + const std::size_t n = point_indices.size(); + double *out = distances.data; + + const std::size_t threads = + bioimage_cpp::detail::normalize_thread_count(n_threads, n); + + // One independent single-source solve per point (row of the matrix). The + // mask/speed inputs are read-only, so the per-thread solves are safe. + bioimage_cpp::detail::parallel_for_chunks( + threads, n, + [&](std::size_t /*thread_id*/, std::size_t begin, std::size_t end) { + detail::GridFastMarching solver(mask, sampling, speed); + for (std::size_t i = begin; i < end; ++i) { + solver.solve({point_indices[i]}); + for (std::size_t j = 0; j < n; ++j) { + out[i * n + j] = solver.distance(point_indices[j]); + } + } + } + ); + + // Symmetrize (per-source solves may differ by a tiny numerical amount) and + // pin the diagonal to exactly zero. + for (std::size_t i = 0; i < n; ++i) { + out[i * n + i] = 0.0; + for (std::size_t j = i + 1; j < n; ++j) { + const double avg = 0.5 * (out[i * n + j] + out[j * n + i]); + out[i * n + j] = avg; + out[j * n + i] = avg; + } + } +} + +} // namespace bioimage_cpp::distance diff --git a/include/bioimage_cpp/distance/geodesic_mesh.hxx b/include/bioimage_cpp/distance/geodesic_mesh.hxx new file mode 100644 index 0000000..fe0b923 --- /dev/null +++ b/include/bioimage_cpp/distance/geodesic_mesh.hxx @@ -0,0 +1,121 @@ +#pragma once + +#include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/threading.hxx" +#include "bioimage_cpp/distance/detail/mesh_fast_marching.hxx" + +#include +#include +#include +#include +#include +#include + +// Geodesic distances on a triangle-mesh surface. +// +// A triangle mesh defines the geometry: distances are shortest-path lengths +// measured across the 2-manifold surface (not through the embedding space). +// Solved with the first-order Kimmel-Sethian fast marching method (see +// detail/mesh_fast_marching.hxx). The exact reference oracle (pygeodesic) lives +// in development/distance/. Vertices carry real coordinates, so there is no +// sampling argument. + +namespace bioimage_cpp::distance { + +namespace detail { + +// Validate 1-D int64 vertex indices against the vertex count and return them as +// linear indices. +inline std::vector vertex_indices( + const ConstArrayView &indices, + std::size_t n_vertices, + const char *name +) { + const auto m = static_cast(indices.shape[0]); + std::vector out(m); + for (std::size_t i = 0; i < m; ++i) { + const std::int64_t v = indices.data[i]; + if (v < 0 || static_cast(v) >= n_vertices) { + throw std::invalid_argument( + std::string(name) + " contains vertex index " + std::to_string(v) + + " out of range [0, " + std::to_string(n_vertices) + ")" + ); + } + out[i] = static_cast(v); + } + return out; +} + +} // namespace detail + +// Geodesic distance field on a triangle-mesh surface from a set of sources. +// +// vertices: (n_vertices, 3) float64. faces: (n_faces, 3) int64 triangle +// indices. sources: (n_sources,) int64 vertex indices. speed: optional +// (n_vertices,) per-vertex speed (nullptr => unit speed). distances (out): +// (n_vertices,) float64; +inf for vertices unreachable from any source. +inline void geodesic_distance_field_mesh( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const ConstArrayView &sources, + const ConstArrayView *speed, + ArrayView &distances, + std::size_t /*n_threads*/ = 1 +) { + detail::MeshFastMarching solver(vertices, faces, speed); + const auto source_indices = + detail::vertex_indices(sources, solver.size(), "sources"); + solver.solve(source_indices); + + const auto &dist = solver.distances(); + std::copy(dist.begin(), dist.end(), distances.data); +} + +// Full pairwise geodesic distance matrix between mesh vertices. +// +// points: (n_points,) int64 vertex indices. distances (out): (n_points, +// n_points) float64, symmetric with a zero diagonal; +inf when two points lie +// in different connected components. +inline void geodesic_distances_mesh( + const ConstArrayView &vertices, + const ConstArrayView &faces, + const ConstArrayView &points, + const ConstArrayView *speed, + ArrayView &distances, + std::size_t n_threads = 1 +) { + // Build one solver up front to validate the mesh and the point indices. + detail::MeshFastMarching probe(vertices, faces, speed); + const auto point_indices = detail::vertex_indices(points, probe.size(), "points"); + const std::size_t n = point_indices.size(); + double *out = distances.data; + + const std::size_t threads = + bioimage_cpp::detail::normalize_thread_count(n_threads, n); + + // One independent single-source solve per point (row of the matrix); the + // mesh/speed inputs are read-only, so the per-thread solves are safe. + bioimage_cpp::detail::parallel_for_chunks( + threads, n, + [&](std::size_t /*thread_id*/, std::size_t begin, std::size_t end) { + detail::MeshFastMarching solver(vertices, faces, speed); + for (std::size_t i = begin; i < end; ++i) { + solver.solve({point_indices[i]}); + for (std::size_t j = 0; j < n; ++j) { + out[i * n + j] = solver.distance(point_indices[j]); + } + } + } + ); + + for (std::size_t i = 0; i < n; ++i) { + out[i * n + i] = 0.0; + for (std::size_t j = i + 1; j < n; ++j) { + const double avg = 0.5 * (out[i * n + j] + out[j * n + i]); + out[i * n + j] = avg; + out[j * n + i] = avg; + } + } +} + +} // namespace bioimage_cpp::distance diff --git a/include/bioimage_cpp/non_maximum_distance_suppression.hxx b/include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx similarity index 100% rename from include/bioimage_cpp/non_maximum_distance_suppression.hxx rename to include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx diff --git a/src/bindings/distance.cxx b/src/bindings/distance.cxx index b05bdb6..087fdde 100644 --- a/src/bindings/distance.cxx +++ b/src/bindings/distance.cxx @@ -1,8 +1,10 @@ #include "distance.hxx" #include "bioimage_cpp/array_view.hxx" -#include "bioimage_cpp/distance_transform.hxx" -#include "bioimage_cpp/non_maximum_distance_suppression.hxx" +#include "bioimage_cpp/distance/distance_transform.hxx" +#include "bioimage_cpp/distance/non_maximum_distance_suppression.hxx" +#include "bioimage_cpp/distance/geodesic_mask.hxx" +#include "bioimage_cpp/distance/geodesic_mesh.hxx" #include #include @@ -25,6 +27,9 @@ namespace { using UInt8Input = nb::ndarray; using FloatArray = nb::ndarray; using Int32Array = nb::ndarray; +using Int64Input = nb::ndarray; +using DoubleInput = nb::ndarray; +using DoubleArray = nb::ndarray; template std::vector ndarray_shape(const Array &array) { @@ -265,6 +270,241 @@ nb::ndarray non_maximum_distance_suppression_im return output; } +// --------------------------------------------------------------------------- +// Geodesic distances (masks + meshes). Interface only for now: the core +// functions throw "not yet implemented" until the solvers land, but the full +// Python -> binding -> core path is wired and validated here. +// --------------------------------------------------------------------------- + +// Build a ConstArrayView over an optional speed ndarray, checking its +// shape against `expected`. Returns nullptr when no speed was supplied. +struct OptionalSpeed { + ConstArrayView view; + const ConstArrayView *ptr = nullptr; +}; + +OptionalSpeed make_optional_speed( + const std::optional &speed, + const std::vector &shape +) { + OptionalSpeed result; + if (speed.has_value()) { + check_buffer_shape("speed", ndarray_shape(*speed), size_t_shape(shape)); + result.view = ConstArrayView{speed->data(), shape, {}}; + result.ptr = &result.view; + } + return result; +} + +nb::tuple geodesic_distance_field_mask( + UInt8Input mask, + Int64Input sources, + const std::vector &sampling, + std::optional speed, + const bool return_gradient, + const std::size_t n_threads +) { + if (mask.ndim() == 0) { + throw std::invalid_argument("mask must have ndim >= 1, got ndim=0"); + } + if (sampling.size() != mask.ndim()) { + throw std::invalid_argument( + "sampling must have length matching mask ndim, got ndim=" + + std::to_string(mask.ndim()) + ", sampling length=" + + std::to_string(sampling.size()) + ); + } + if (sources.ndim() != 2) { + throw std::invalid_argument( + "sources must have ndim == 2, got ndim=" + std::to_string(sources.ndim()) + ); + } + if (static_cast(sources.shape(1)) != mask.ndim()) { + throw std::invalid_argument( + "sources.shape[1] must match mask ndim, got sources.shape[1]=" + + std::to_string(sources.shape(1)) + ", mask.ndim()=" + + std::to_string(mask.ndim()) + ); + } + + const auto shape = ndarray_shape(mask); + const auto ndim = shape.size(); + auto distances_array = make_array(size_t_shape(shape)); + + ConstArrayView mask_view{mask.data(), shape, {}}; + ConstArrayView sources_view{sources.data(), ndarray_shape(sources), {}}; + ArrayView distances_view{distances_array.data(), shape, {}}; + const OptionalSpeed speed_opt = make_optional_speed(speed, shape); + + // Optional per-axis gradient output, shape (*mask.shape, ndim), float32. + FloatArray gradient_array; + ArrayView gradient_view; + ArrayView *gradient_ptr = nullptr; + if (return_gradient) { + auto gradient_shape = size_t_shape(shape); + gradient_shape.push_back(ndim); + gradient_array = make_array(gradient_shape); + auto gradient_view_shape = shape; + gradient_view_shape.push_back(static_cast(ndim)); + gradient_view = ArrayView{gradient_array.data(), gradient_view_shape, {}}; + gradient_ptr = &gradient_view; + } + + { + nb::gil_scoped_release release; + distance::geodesic_distance_field( + mask_view, sources_view, sampling, speed_opt.ptr, distances_view, n_threads, + gradient_ptr + ); + } + + auto gradient_result = + return_gradient ? nb::cast(gradient_array) : nb::object(nb::none()); + return nb::make_tuple(nb::cast(distances_array), gradient_result); +} + +DoubleArray geodesic_distances_mask( + UInt8Input mask, + Int64Input points, + const std::vector &sampling, + std::optional speed, + const std::size_t n_threads +) { + if (mask.ndim() == 0) { + throw std::invalid_argument("mask must have ndim >= 1, got ndim=0"); + } + if (sampling.size() != mask.ndim()) { + throw std::invalid_argument( + "sampling must have length matching mask ndim, got ndim=" + + std::to_string(mask.ndim()) + ", sampling length=" + + std::to_string(sampling.size()) + ); + } + if (points.ndim() != 2) { + throw std::invalid_argument( + "points must have ndim == 2, got ndim=" + std::to_string(points.ndim()) + ); + } + if (static_cast(points.shape(1)) != mask.ndim()) { + throw std::invalid_argument( + "points.shape[1] must match mask ndim, got points.shape[1]=" + + std::to_string(points.shape(1)) + ", mask.ndim()=" + + std::to_string(mask.ndim()) + ); + } + + const auto shape = ndarray_shape(mask); + const auto n_points = static_cast(points.shape(0)); + auto distances_array = make_array({n_points, n_points}); + + ConstArrayView mask_view{mask.data(), shape, {}}; + ConstArrayView points_view{points.data(), ndarray_shape(points), {}}; + ArrayView distances_view{ + distances_array.data(), + {static_cast(n_points), static_cast(n_points)}, + {}, + }; + const OptionalSpeed speed_opt = make_optional_speed(speed, shape); + + { + nb::gil_scoped_release release; + distance::geodesic_distances( + mask_view, points_view, sampling, speed_opt.ptr, distances_view, n_threads + ); + } + return distances_array; +} + +// Validate a triangle mesh (vertices (V, 3) float64, faces (F, 3) int64) and +// return the vertex count. +std::size_t check_mesh(const DoubleInput &vertices, const Int64Input &faces) { + if (vertices.ndim() != 2 || vertices.shape(1) != 3) { + throw std::invalid_argument( + "vertices must have shape (n_vertices, 3), got ndim=" + + std::to_string(vertices.ndim()) + ); + } + if (faces.ndim() != 2 || faces.shape(1) != 3) { + throw std::invalid_argument( + "faces must have shape (n_faces, 3), got ndim=" + std::to_string(faces.ndim()) + ); + } + return static_cast(vertices.shape(0)); +} + +DoubleArray geodesic_distance_field_mesh( + DoubleInput vertices, + Int64Input faces, + Int64Input sources, + std::optional speed, + const std::size_t n_threads +) { + const auto n_vertices = check_mesh(vertices, faces); + if (sources.ndim() != 1) { + throw std::invalid_argument( + "sources must have ndim == 1 (vertex indices), got ndim=" + + std::to_string(sources.ndim()) + ); + } + + auto distances_array = make_array({n_vertices}); + + ConstArrayView vertices_view{vertices.data(), ndarray_shape(vertices), {}}; + ConstArrayView faces_view{faces.data(), ndarray_shape(faces), {}}; + ConstArrayView sources_view{sources.data(), ndarray_shape(sources), {}}; + ArrayView distances_view{ + distances_array.data(), {static_cast(n_vertices)}, {} + }; + const OptionalSpeed speed_opt = + make_optional_speed(speed, {static_cast(n_vertices)}); + + { + nb::gil_scoped_release release; + distance::geodesic_distance_field_mesh( + vertices_view, faces_view, sources_view, speed_opt.ptr, distances_view, n_threads + ); + } + return distances_array; +} + +DoubleArray geodesic_distances_mesh( + DoubleInput vertices, + Int64Input faces, + Int64Input points, + std::optional speed, + const std::size_t n_threads +) { + const auto n_vertices = check_mesh(vertices, faces); + if (points.ndim() != 1) { + throw std::invalid_argument( + "points must have ndim == 1 (vertex indices), got ndim=" + + std::to_string(points.ndim()) + ); + } + + const auto n_points = static_cast(points.shape(0)); + auto distances_array = make_array({n_points, n_points}); + + ConstArrayView vertices_view{vertices.data(), ndarray_shape(vertices), {}}; + ConstArrayView faces_view{faces.data(), ndarray_shape(faces), {}}; + ConstArrayView points_view{points.data(), ndarray_shape(points), {}}; + ArrayView distances_view{ + distances_array.data(), + {static_cast(n_points), static_cast(n_points)}, + {}, + }; + const OptionalSpeed speed_opt = + make_optional_speed(speed, {static_cast(n_vertices)}); + + { + nb::gil_scoped_release release; + distance::geodesic_distances_mesh( + vertices_view, faces_view, points_view, speed_opt.ptr, distances_view, n_threads + ); + } + return distances_array; +} + } // namespace void bind_distance(nb::module_ &m) { @@ -311,6 +551,45 @@ void bind_distance(nb::module_ &m) { &non_maximum_distance_suppression_impl, nb::arg("distance_map"), nb::arg("points"), nb::arg("n_threads"), nms_doc ); + + m.def( + "_geodesic_distance_field_mask", + &geodesic_distance_field_mask, + nb::arg("mask"), nb::arg("sources"), nb::arg("sampling"), + nb::arg("speed").none(), nb::arg("return_gradient"), nb::arg("n_threads"), + "Geodesic distance field within a mask from a set of source coordinates.\n" + "mask nonzero = inside the domain. sources is (n_sources, ndim) int64.\n" + "Returns (field, gradient): field is float64 of mask.shape (unreachable\n" + "voxels +inf); gradient is float32 (*mask.shape, ndim) when\n" + "return_gradient else None." + ); + m.def( + "_geodesic_distances_mask", + &geodesic_distances_mask, + nb::arg("mask"), nb::arg("points"), nb::arg("sampling"), + nb::arg("speed").none(), nb::arg("n_threads"), + "Full pairwise geodesic distance matrix between points within a mask.\n" + "points is (n_points, ndim) int64. Returns a symmetric (n_points,\n" + "n_points) float64 matrix; +inf where two points are not connected." + ); + m.def( + "_geodesic_distance_field_mesh", + &geodesic_distance_field_mesh, + nb::arg("vertices"), nb::arg("faces"), nb::arg("sources"), + nb::arg("speed").none(), nb::arg("n_threads"), + "Geodesic distance field on a triangle mesh from a set of source\n" + "vertices. vertices (n_vertices, 3) float64, faces (n_faces, 3) int64,\n" + "sources (n_sources,) int64 vertex indices. Returns (n_vertices,) float64." + ); + m.def( + "_geodesic_distances_mesh", + &geodesic_distances_mesh, + nb::arg("vertices"), nb::arg("faces"), nb::arg("points"), + nb::arg("speed").none(), nb::arg("n_threads"), + "Full pairwise geodesic distance matrix between mesh vertices.\n" + "points is (n_points,) int64 vertex indices. Returns a symmetric\n" + "(n_points, n_points) float64 matrix." + ); } } // namespace bioimage_cpp::bindings diff --git a/src/bioimage_cpp/_version.py b/src/bioimage_cpp/_version.py index 3d18726..906d362 100644 --- a/src/bioimage_cpp/_version.py +++ b/src/bioimage_cpp/_version.py @@ -1 +1 @@ -__version__ = "0.5.0" +__version__ = "0.6.0" diff --git a/src/bioimage_cpp/distance/__init__.py b/src/bioimage_cpp/distance/__init__.py index ac0e836..22fd436 100644 --- a/src/bioimage_cpp/distance/__init__.py +++ b/src/bioimage_cpp/distance/__init__.py @@ -1,13 +1,25 @@ -"""Distance transforms.""" +"""Distance transforms and geodesic distances.""" from ._distance import ( distance_transform, non_maximum_distance_suppression, vector_difference_transform, ) +from ._geodesic import ( + geodesic_distance_field, + geodesic_distance_field_mesh, + geodesic_distances, + geodesic_distances_mesh, + geodesic_gradient_field, +) __all__ = [ "distance_transform", "non_maximum_distance_suppression", "vector_difference_transform", + "geodesic_distance_field", + "geodesic_gradient_field", + "geodesic_distances", + "geodesic_distance_field_mesh", + "geodesic_distances_mesh", ] diff --git a/src/bioimage_cpp/distance/_geodesic.py b/src/bioimage_cpp/distance/_geodesic.py new file mode 100644 index 0000000..3100d59 --- /dev/null +++ b/src/bioimage_cpp/distance/_geodesic.py @@ -0,0 +1,310 @@ +"""Python wrappers for the geodesic-distance bindings. + +Geodesic distance is the shortest-path length constrained to a geometry: + +- **masks** (regular grid): distances stay inside the nonzero region and never + cross background voxels (fast-marching / Eikonal formulation, à la + scikit-fmm). +- **surfaces** (triangle meshes): distances are measured across the mesh + surface (exact MMP geodesics, à la pygeodesic). + +.. note:: + + The C++ solvers are not implemented yet. These wrappers validate their + arguments fully and dispatch to the ``_core`` bindings, which currently + raise ``RuntimeError("... not yet implemented")``. The reference behaviour + lives in ``development/distance/``. +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np + +from .. import _core +from ._distance import _as_binary_input, _normalize_sampling, _normalize_threads + + +def _as_coordinates(points: np.ndarray, ndim: int, function: str, name: str) -> np.ndarray: + """Coerce a point set to a C-contiguous ``(n, ndim)`` int64 array.""" + array = np.ascontiguousarray(points) + if array.ndim == 1 and array.size == ndim: + # A single coordinate may be given as a flat (ndim,) vector. + array = array.reshape(1, ndim) + if array.ndim != 2: + raise ValueError( + f"{function}: {name} must have shape (n, {ndim}), got ndim={array.ndim}" + ) + if array.shape[1] != ndim: + raise ValueError( + f"{function}: {name}.shape[1] must equal the mask ndim ({ndim}), " + f"got {name}.shape[1]={array.shape[1]}" + ) + return np.ascontiguousarray(array, dtype=np.int64) + + +def _as_vertex_indices(indices: np.ndarray, function: str, name: str) -> np.ndarray: + """Coerce vertex indices to a C-contiguous 1-D int64 array.""" + array = np.atleast_1d(np.ascontiguousarray(indices)) + if array.ndim != 1: + raise ValueError( + f"{function}: {name} must be 1-D vertex indices, got ndim={array.ndim}" + ) + return np.ascontiguousarray(array, dtype=np.int64) + + +def _as_mesh(vertices: np.ndarray, faces: np.ndarray, function: str): + """Coerce a triangle mesh to ``(vertices float64 (V,3), faces int64 (F,3))``.""" + vertices = np.ascontiguousarray(vertices, dtype=np.float64) + if vertices.ndim != 2 or vertices.shape[1] != 3: + raise ValueError( + f"{function}: vertices must have shape (n_vertices, 3), " + f"got shape={vertices.shape}" + ) + faces = np.ascontiguousarray(faces, dtype=np.int64) + if faces.ndim != 2 or faces.shape[1] != 3: + raise ValueError( + f"{function}: faces must have shape (n_faces, 3), got shape={faces.shape}" + ) + return vertices, faces + + +def _as_speed( + speed: np.ndarray | None, expected_shape: tuple[int, ...], function: str +) -> np.ndarray | None: + """Coerce an optional speed field to a C-contiguous float64 array.""" + if speed is None: + return None + array = np.ascontiguousarray(speed, dtype=np.float64) + if array.shape != tuple(expected_shape): + raise ValueError( + f"{function}: speed must have shape {tuple(expected_shape)}, " + f"got shape={array.shape}" + ) + return array + + +def geodesic_distance_field( + mask: np.ndarray, + sources: np.ndarray, + sampling: float | Sequence[float] | None = None, + speed: np.ndarray | None = None, + return_gradient: bool = False, + number_of_threads: int = 1, +): + """Geodesic distance field within a mask from a set of source coordinates. + + Computes, for every foreground voxel, the shortest-path distance to the + nearest source, where paths must stay inside the mask (never crossing + background). A single source is just a one-row ``sources`` array. + + Parameters + ---------- + mask + Binary/label array of any ndim. Nonzero is inside the domain; distances + propagate only through nonzero voxels. Coerced to C-contiguous + ``uint8``. + sources + Integer array of shape ``(n_sources, ndim)`` (or a flat ``(ndim,)`` + vector for a single source); each row is a voxel coordinate in NumPy + axis order. + sampling + Per-axis voxel spacing. Scalar or per-axis sequence; default 1.0. + speed + Optional per-voxel speed, same shape as ``mask``. ``None`` gives + unit-speed geodesic distance; otherwise the result is the weighted + travel time. + return_gradient + If ``True``, also return the per-axis gradient of the field (see + Returns). Analogous to :func:`vector_difference_transform`. + number_of_threads + ``0`` uses ``hardware_concurrency``; a positive value pins the thread + count. Default ``1``. + + Returns + ------- + np.ndarray or (np.ndarray, np.ndarray) + The distance field, a ``float64`` array of shape ``mask.shape``. + Background voxels and voxels unreachable from any source are ``+inf``. + If ``return_gradient`` is ``True``, returns ``(field, gradient)`` where + ``gradient`` is a ``float32`` array of shape ``(*mask.shape, ndim)`` + holding the first-order upwind gradient ``d(field)/d(axis)`` (channel + last, NumPy axis order). The gradient points **away from the source** + (direction of increasing distance) with ``norm ~= 1/speed``; negate it + to trace back toward the source. It is zero at sources, background, and + unreachable voxels. + """ + function = "geodesic_distance_field" + binary = _as_binary_input(mask, function) + ndim = binary.ndim + sources_arr = _as_coordinates(sources, ndim, function, "sources") + sampling_values = _normalize_sampling(sampling, ndim, function) + speed_arr = _as_speed(speed, tuple(binary.shape), function) + n_threads = _normalize_threads(number_of_threads, function) + field, gradient = _core._geodesic_distance_field_mask( + binary, sources_arr, sampling_values, speed_arr, bool(return_gradient), n_threads + ) + if return_gradient: + return field, gradient + return field + + +def geodesic_gradient_field( + mask: np.ndarray, + sources: np.ndarray, + sampling: float | Sequence[float] | None = None, + speed: np.ndarray | None = None, + number_of_threads: int = 1, +) -> np.ndarray: + """Return the per-axis gradient of the geodesic distance field within a mask. + + Thin wrapper around :func:`geodesic_distance_field` with + ``return_gradient=True`` that returns only the gradient. Output has shape + ``mask.shape + (ndim,)`` and dtype ``float32``; the trailing vector axis + follows NumPy axis order. Each component is ``d(field)/d(axis)`` pointing + away from the nearest source (``norm ~= 1/speed``); negate to trace toward + it (e.g. to feed :func:`bioimage_cpp.flow.compute_flow_density`). + """ + _, gradient = geodesic_distance_field( + mask, + sources, + sampling=sampling, + speed=speed, + return_gradient=True, + number_of_threads=number_of_threads, + ) + return gradient + + +def geodesic_distances( + mask: np.ndarray, + points: np.ndarray, + sampling: float | Sequence[float] | None = None, + speed: np.ndarray | None = None, + number_of_threads: int = 1, +) -> np.ndarray: + """Full pairwise geodesic distance matrix between points within a mask. + + Parameters + ---------- + mask + Binary/label array of any ndim; nonzero is inside the domain. Coerced + to C-contiguous ``uint8``. + points + Integer array of shape ``(n_points, ndim)``; each row is a voxel + coordinate in NumPy axis order. + sampling + Per-axis voxel spacing. Scalar or per-axis sequence; default 1.0. + speed + Optional per-voxel speed, same shape as ``mask``. ``None`` gives + unit-speed distances. + number_of_threads + ``0`` uses ``hardware_concurrency``; a positive value pins the thread + count. Default ``1``. + + Returns + ------- + np.ndarray + Symmetric ``float64`` matrix of shape ``(n_points, n_points)``. Entry + ``(i, j)`` is the geodesic distance from ``points[i]`` to ``points[j]`` + within the mask, ``+inf`` when they are not connected inside the domain, + and ``0`` on the diagonal. + """ + function = "geodesic_distances" + binary = _as_binary_input(mask, function) + ndim = binary.ndim + points_arr = _as_coordinates(points, ndim, function, "points") + sampling_values = _normalize_sampling(sampling, ndim, function) + speed_arr = _as_speed(speed, tuple(binary.shape), function) + n_threads = _normalize_threads(number_of_threads, function) + return _core._geodesic_distances_mask( + binary, points_arr, sampling_values, speed_arr, n_threads + ) + + +def geodesic_distance_field_mesh( + vertices: np.ndarray, + faces: np.ndarray, + sources: np.ndarray, + speed: np.ndarray | None = None, + number_of_threads: int = 1, +) -> np.ndarray: + """Geodesic distance field on a triangle-mesh surface from source vertices. + + Parameters + ---------- + vertices + Vertex positions of shape ``(n_vertices, 3)``. Coerced to + C-contiguous ``float64``. + faces + Triangle vertex indices of shape ``(n_faces, 3)``. Coerced to + C-contiguous ``int64``. + sources + 1-D array of source vertex indices (a scalar/one-element array for a + single source). + speed + Optional per-vertex speed of shape ``(n_vertices,)``. ``None`` gives + unit-speed geodesic distance. + number_of_threads + ``0`` uses ``hardware_concurrency``; a positive value pins the thread + count. Default ``1``. + + Returns + ------- + np.ndarray + ``float64`` array of shape ``(n_vertices,)``. Vertices unreachable from + any source (a disconnected component) are ``+inf``. + """ + function = "geodesic_distance_field_mesh" + vertices_arr, faces_arr = _as_mesh(vertices, faces, function) + sources_arr = _as_vertex_indices(sources, function, "sources") + speed_arr = _as_speed(speed, (vertices_arr.shape[0],), function) + n_threads = _normalize_threads(number_of_threads, function) + return _core._geodesic_distance_field_mesh( + vertices_arr, faces_arr, sources_arr, speed_arr, n_threads + ) + + +def geodesic_distances_mesh( + vertices: np.ndarray, + faces: np.ndarray, + points: np.ndarray, + speed: np.ndarray | None = None, + number_of_threads: int = 1, +) -> np.ndarray: + """Full pairwise geodesic distance matrix between mesh vertices. + + Parameters + ---------- + vertices + Vertex positions of shape ``(n_vertices, 3)``. Coerced to + C-contiguous ``float64``. + faces + Triangle vertex indices of shape ``(n_faces, 3)``. Coerced to + C-contiguous ``int64``. + points + 1-D array of vertex indices. + speed + Optional per-vertex speed of shape ``(n_vertices,)``. ``None`` gives + unit-speed distances. + number_of_threads + ``0`` uses ``hardware_concurrency``; a positive value pins the thread + count. Default ``1``. + + Returns + ------- + np.ndarray + Symmetric ``float64`` matrix of shape ``(n_points, n_points)``. Entry + ``(i, j)`` is the surface geodesic distance from ``points[i]`` to + ``points[j]``, ``+inf`` when they lie in different connected components, + and ``0`` on the diagonal. + """ + function = "geodesic_distances_mesh" + vertices_arr, faces_arr = _as_mesh(vertices, faces, function) + points_arr = _as_vertex_indices(points, function, "points") + speed_arr = _as_speed(speed, (vertices_arr.shape[0],), function) + n_threads = _normalize_threads(number_of_threads, function) + return _core._geodesic_distances_mesh( + vertices_arr, faces_arr, points_arr, speed_arr, n_threads + ) diff --git a/tests/distance/test_geodesic_distance.py b/tests/distance/test_geodesic_distance.py new file mode 100644 index 0000000..3441212 --- /dev/null +++ b/tests/distance/test_geodesic_distance.py @@ -0,0 +1,438 @@ +"""Tests for the geodesic-distance solvers (masks + triangle meshes). + +The solvers use a first-order fast marching method, so these tests assert the +properties a correct geodesic solver must satisfy (source distance 0, +monotonicity, ``geodesic >= euclidean``, exact behaviour along grid axes, exact +speed scaling, symmetry of pairwise matrices, ``+inf`` for unreachable regions) +rather than tight agreement with an exact Euclidean field (which first-order FMM +overestimates by up to ~20% near diagonals). Tight numeric agreement with the +scikit-fmm / pygeodesic references lives in ``development/distance/`` and in the +optional, dependency-guarded cross-checks at the end of this file. +""" + +import numpy as np +import pytest + +import bioimage_cpp as bic + +INF = np.inf + + +# --------------------------------------------------------------------------- # +# helpers +# --------------------------------------------------------------------------- # + + +def euclidean_from(shape, source, sampling=None): + """Analytic Euclidean distance from a single source voxel over a full grid.""" + sampling = np.ones(len(shape)) if sampling is None else np.asarray(sampling, float) + coords = np.indices(shape, dtype=np.float64) + total = np.zeros(shape, dtype=np.float64) + for axis, s in enumerate(source): + total += ((coords[axis] - s) * sampling[axis]) ** 2 + return np.sqrt(total) + + +def flat_grid_mesh(n): + """A planar (z=0) triangulated ``n x n`` grid; geodesic == 2D Euclidean.""" + xs, ys = np.meshgrid(np.arange(n), np.arange(n), indexing="ij") + verts = np.stack([xs.ravel(), ys.ravel(), np.zeros(n * n)], axis=1).astype(np.float64) + + def vid(i, j): + return i * n + j + + faces = [] + for i in range(n - 1): + for j in range(n - 1): + faces.append([vid(i, j), vid(i + 1, j), vid(i, j + 1)]) + faces.append([vid(i + 1, j), vid(i + 1, j + 1), vid(i, j + 1)]) + return verts, np.array(faces, np.int64), vid + + +def sphere_mesh(n_points=300, radius=5.0, seed=0): + """A closed sphere mesh from the convex hull of points on the sphere.""" + hull = pytest.importorskip("scipy.spatial").ConvexHull + rng = np.random.default_rng(seed) + pts = rng.standard_normal((n_points, 3)) + pts /= np.linalg.norm(pts, axis=1, keepdims=True) + pts *= radius + faces = hull(pts).simplices.astype(np.int64) + return np.ascontiguousarray(pts, np.float64), np.ascontiguousarray(faces) + + +# --------------------------------------------------------------------------- # +# mask: field +# --------------------------------------------------------------------------- # + + +@pytest.mark.parametrize("shape", [(25, 25), (9, 10, 11)]) +def test_mask_field_basic_properties(shape): + mask = np.ones(shape, np.uint8) + source = tuple(0 for _ in shape) + field = bic.distance.geodesic_distance_field(mask, np.array([source], np.int64)) + + assert field.shape == shape + assert field.dtype == np.float64 + assert field[source] == 0.0 + assert np.all(np.isfinite(field)) + # geodesic in an obstacle-free domain is >= the straight-line distance + eucl = euclidean_from(shape, source) + assert np.all(field >= eucl - 1e-9) + + +def test_mask_field_axis_aligned_exact(): + mask = np.ones((30, 30), np.uint8) + field = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + np.testing.assert_allclose(field[0, :], np.arange(30), atol=1e-9) + np.testing.assert_allclose(field[:, 0], np.arange(30), atol=1e-9) + + +def test_mask_field_farfield_close_to_euclidean(): + mask = np.ones((60, 60), np.uint8) + field = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + eucl = euclidean_from((60, 60), (0, 0)) + far = eucl > 15 + rel = np.abs(field[far] - eucl[far]) / eucl[far] + assert rel.max() < 0.1 # first-order FMM error is modest in the far field + + +def test_mask_anisotropic_sampling_axis_exact(): + mask = np.ones((20, 20), np.uint8) + field = bic.distance.geodesic_distance_field( + mask, np.array([[0, 0]], np.int64), sampling=(2.0, 0.5) + ) + np.testing.assert_allclose(field[:, 0], 2.0 * np.arange(20), atol=1e-9) + np.testing.assert_allclose(field[0, :], 0.5 * np.arange(20), atol=1e-9) + + +def test_mask_speed_scaling_is_exact(): + mask = np.ones((32, 32), np.uint8) + src = np.array([[3, 5]], np.int64) + base = bic.distance.geodesic_distance_field(mask, src) + scaled = bic.distance.geodesic_distance_field(mask, src, speed=2.0 * np.ones((32, 32))) + np.testing.assert_allclose(scaled, base / 2.0, atol=1e-9) + + +def test_mask_multi_source_monotone(): + mask = np.ones((30, 30), np.uint8) + a = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + b = bic.distance.geodesic_distance_field(mask, np.array([[29, 29]], np.int64)) + both = bic.distance.geodesic_distance_field( + mask, np.array([[0, 0], [29, 29]], np.int64) + ) + # adding a source can only lower arrival times, so the two-source field is + # bounded above by the pointwise min of the single-source fields; it dips + # slightly below near the medial axis where the update stencil mixes fronts. + assert np.all(both <= np.minimum(a, b) + 1e-9) + np.testing.assert_allclose(both, np.minimum(a, b), atol=0.5) + + +def test_mask_single_source_as_flat_vector(): + mask = np.ones((10, 10), np.uint8) + field = bic.distance.geodesic_distance_field(mask, np.array([2, 3], np.int64)) + assert field[2, 3] == 0.0 + + +def test_mask_obstacle_detour_and_unreachable(): + # Wall across the middle with a gap on the right; source top-left. + mask = np.ones((21, 21), np.uint8) + mask[10, :18] = 0 # gap at columns 18, 19, 20 + field = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + + # wall voxels are background -> +inf + assert np.isinf(field[10, 0]) + # target directly across the wall must detour around the gap + target = (20, 0) + straight = np.hypot(20, 0) + assert np.isfinite(field[target]) + assert field[target] > straight + 1.0 + # the detour is bounded below by the best two-leg path through a gap cell + gaps = [(10, c) for c in range(18, 21)] + ref = min(np.hypot(*np.subtract((0, 0), g)) + np.hypot(*np.subtract(g, target)) for g in gaps) + assert 0.95 * ref <= field[target] <= 1.25 * ref + + +def test_mask_isolated_region_is_inf(): + mask = np.zeros((10, 10), np.uint8) + mask[0:3, 0:3] = 1 # region containing the source + mask[7:10, 7:10] = 1 # disconnected region + field = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + assert np.all(np.isinf(field[7:10, 7:10])) + assert np.all(np.isfinite(field[0:3, 0:3])) + + +# --------------------------------------------------------------------------- # +# mask: pairwise +# --------------------------------------------------------------------------- # + + +def test_mask_pairwise_symmetric_zero_diagonal(): + mask = np.ones((20, 20), np.uint8) + points = np.array([[0, 0], [19, 0], [0, 19], [19, 19], [10, 10]], np.int64) + D = bic.distance.geodesic_distances(mask, points) + assert D.shape == (5, 5) + assert D.dtype == np.float64 + np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12) + np.testing.assert_allclose(D, D.T, atol=1e-9) + assert np.all(D >= 0.0) + + +def test_mask_pairwise_matches_field_rows(): + mask = np.ones((18, 18), np.uint8) + points = np.array([[0, 0], [17, 5], [3, 12]], np.int64) + D = bic.distance.geodesic_distances(mask, points) + for i, p in enumerate(points): + field = bic.distance.geodesic_distance_field(mask, p[None, :]) + row = np.array([field[tuple(q)] for q in points]) + np.testing.assert_allclose(D[i], row, atol=1e-9) + + +# --------------------------------------------------------------------------- # +# mask: validation +# --------------------------------------------------------------------------- # + + +def test_mask_invalid_arguments(): + mask = np.ones((8, 8), np.uint8) + with pytest.raises(ValueError, match="out of bounds"): + bic.distance.geodesic_distance_field(mask, np.array([[8, 0]], np.int64)) + with pytest.raises(ValueError, match="out of bounds"): + bic.distance.geodesic_distance_field(mask, np.array([[-1, 0]], np.int64)) + with pytest.raises(ValueError, match="speed"): + bic.distance.geodesic_distance_field( + mask, np.array([[0, 0]], np.int64), speed=np.ones((3, 3)) + ) + with pytest.raises(ValueError, match="sources.shape"): + bic.distance.geodesic_distance_field(mask, np.zeros((1, 3), np.int64)) + + +# --------------------------------------------------------------------------- # +# mask: gradient of the field +# --------------------------------------------------------------------------- # + + +def _interior_mask(shape, source, margin=3, min_radius=8): + """Voxels away from the source and the array border (smooth-field region).""" + coords = np.indices(shape, dtype=np.float64) + eucl = np.sqrt(sum((coords[a] - source[a]) ** 2 for a in range(len(shape)))) + sel = eucl > min_radius + for a, size in enumerate(shape): + sel &= (coords[a] > margin) & (coords[a] < size - 1 - margin) + return sel + + +def test_mask_gradient_shape_dtype_and_backcompat(): + mask = np.ones((20, 20), np.uint8) + src = np.array([[0, 0]], np.int64) + + field_only = bic.distance.geodesic_distance_field(mask, src) + assert isinstance(field_only, np.ndarray) + assert field_only.dtype == np.float64 + + field, grad = bic.distance.geodesic_distance_field(mask, src, return_gradient=True) + assert np.array_equal(field, field_only) # field is unchanged by the option + assert grad.shape == mask.shape + (mask.ndim,) + assert grad.dtype == np.float32 + + +@pytest.mark.parametrize("shape", [(60, 60), (18, 30, 30)]) +def test_mask_gradient_eikonal_norm(shape): + mask = np.ones(shape, np.uint8) + source = tuple(0 for _ in shape) + _, grad = bic.distance.geodesic_distance_field( + mask, np.array([source], np.int64), return_gradient=True + ) + norm = np.linalg.norm(grad, axis=-1) + interior = _interior_mask(shape, source) + # |grad(T)| = slowness = 1 for unit speed (the Eikonal equation) + np.testing.assert_allclose(norm[interior], 1.0, atol=0.02) + + # speed = 2 -> |grad| = 0.5 + _, grad2 = bic.distance.geodesic_distance_field( + mask, np.array([source], np.int64), speed=2.0 * np.ones(shape), return_gradient=True + ) + np.testing.assert_allclose(np.linalg.norm(grad2, axis=-1)[interior], 0.5, atol=0.02) + + +def test_mask_gradient_points_away_from_source(): + shape = (60, 60) + source = (0, 0) + mask = np.ones(shape, np.uint8) + _, grad = bic.distance.geodesic_distance_field( + mask, np.array([source], np.int64), return_gradient=True + ) + coords = np.indices(shape, dtype=np.float64) + # dot(grad, position - source) > 0 => grad points away from the source + dot = grad[..., 0] * coords[0] + grad[..., 1] * coords[1] + interior = _interior_mask(shape, source) + assert np.all(dot[interior] > 0) + # zero at the source itself (a local minimum) + np.testing.assert_allclose(grad[source], 0.0, atol=1e-6) + + +def test_mask_gradient_anisotropic_norm(): + shape = (50, 50) + source = (0, 0) + mask = np.ones(shape, np.uint8) + _, grad = bic.distance.geodesic_distance_field( + mask, np.array([source], np.int64), sampling=(2.0, 0.5), return_gradient=True + ) + interior = _interior_mask(shape, source) + np.testing.assert_allclose(np.linalg.norm(grad, axis=-1)[interior], 1.0, atol=0.02) + + +def test_mask_gradient_field_wrapper_matches(): + mask = np.ones((30, 30), np.uint8) + src = np.array([[5, 7]], np.int64) + _, grad = bic.distance.geodesic_distance_field(mask, src, return_gradient=True) + only = bic.distance.geodesic_gradient_field(mask, src) + assert np.array_equal(only, grad) + assert only.dtype == np.float32 + + +def test_mask_gradient_obstacle_and_background_zero(): + mask = np.ones((21, 21), np.uint8) + mask[10, :18] = 0 # wall with a gap on the right + field, grad = bic.distance.geodesic_distance_field( + mask, np.array([[0, 0]], np.int64), return_gradient=True + ) + # background (wall) voxels have a zero gradient + assert np.all(grad[mask == 0] == 0.0) + # gradient is finite everywhere and unit-norm across most of the reachable + # interior (it may deviate only on the medial axis where fronts meet) + assert np.all(np.isfinite(grad)) + reachable = np.isfinite(field) & (mask != 0) + norm = np.linalg.norm(grad, axis=-1) + interior = reachable & (field > 5) + frac_unit = np.mean(np.abs(norm[interior] - 1.0) < 0.05) + assert frac_unit > 0.95 + + +# --------------------------------------------------------------------------- # +# mesh: field +# --------------------------------------------------------------------------- # + + +def test_mesh_flat_grid_matches_euclidean(): + verts, faces, vid = flat_grid_mesh(21) + field = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([vid(0, 0)], np.int64)) + eucl = np.hypot(verts[:, 0], verts[:, 1]) + + assert field[vid(0, 0)] == 0.0 + assert np.all(np.isfinite(field)) + assert np.all(field >= eucl - 1e-9) + # exact along a mesh edge direction + axis = np.array([field[vid(i, 0)] for i in range(21)]) + np.testing.assert_allclose(axis, np.arange(21), atol=1e-9) + # close to Euclidean in the far field + far = eucl > 8 + rel = np.abs(field[far] - eucl[far]) / eucl[far] + assert rel.max() < 0.12 + + +def test_mesh_field_ge_chord_on_sphere(): + verts, faces = sphere_mesh() + field = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([0], np.int64)) + assert field[0] == 0.0 + assert np.all(np.isfinite(field)) + # surface geodesic distance is >= the straight-line (chord) distance + chord = np.linalg.norm(verts - verts[0], axis=1) + assert np.all(field >= chord - 1e-6) + + +def test_mesh_speed_scaling_is_exact(): + verts, faces = sphere_mesh() + src = np.array([0], np.int64) + base = bic.distance.geodesic_distance_field_mesh(verts, faces, src) + scaled = bic.distance.geodesic_distance_field_mesh( + verts, faces, src, speed=2.0 * np.ones(len(verts)) + ) + np.testing.assert_allclose(scaled, base / 2.0, rtol=1e-9, atol=1e-9) + + +def test_mesh_multi_source_monotone(): + verts, faces = sphere_mesh() + a = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([0], np.int64)) + b = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([50], np.int64)) + both = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([0, 50], np.int64)) + assert np.all(both <= np.minimum(a, b) + 1e-9) + np.testing.assert_allclose(both, np.minimum(a, b), atol=0.3) + + +def test_mesh_disconnected_components_are_inf(): + # two disjoint triangles (no shared vertices) + verts = np.array( + [[0, 0, 0], [1, 0, 0], [0, 1, 0], [10, 0, 0], [11, 0, 0], [10, 1, 0]], + np.float64, + ) + faces = np.array([[0, 1, 2], [3, 4, 5]], np.int64) + field = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([0], np.int64)) + assert np.all(np.isfinite(field[:3])) + assert np.all(np.isinf(field[3:])) + + +# --------------------------------------------------------------------------- # +# mesh: pairwise + validation +# --------------------------------------------------------------------------- # + + +def test_mesh_pairwise_symmetric_and_matches_field(): + verts, faces = sphere_mesh() + points = np.array([0, 40, 120, 200], np.int64) + D = bic.distance.geodesic_distances_mesh(verts, faces, points) + assert D.shape == (4, 4) + np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12) + # the matrix is symmetrized; single-source FMM on the irregular mesh is + # slightly direction-dependent, so rows match the raw field only up to that + # small asymmetry. + np.testing.assert_allclose(D, D.T, atol=1e-9) + for i, p in enumerate(points): + field = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([p], np.int64)) + np.testing.assert_allclose(D[i], field[points], rtol=0.03, atol=0.2) + + +def test_mesh_invalid_arguments(): + verts, faces = sphere_mesh(n_points=60) + with pytest.raises(ValueError, match="out of range"): + bic.distance.geodesic_distance_field_mesh( + verts, faces, np.array([len(verts)], np.int64) + ) + with pytest.raises(ValueError, match="vertices"): + bic.distance.geodesic_distance_field_mesh( + np.zeros((5, 2)), faces, np.array([0], np.int64) + ) + + +# --------------------------------------------------------------------------- # +# optional tight cross-checks against the external references +# --------------------------------------------------------------------------- # + + +def test_mask_matches_scikit_fmm(): + skfmm = pytest.importorskip("skfmm") + mask = np.ones((50, 50), np.uint8) + field = bic.distance.geodesic_distance_field(mask, np.array([[0, 0]], np.int64)) + phi = np.ones((50, 50)) + phi[0, 0] = -1 + ref = np.abs(skfmm.distance(phi, order=1)) + eucl = euclidean_from((50, 50), (0, 0)) + interior = eucl > 3 + # scikit-fmm's level-set seed idiom is offset ~0.5 cell; remove it, then the + # two first-order schemes must agree to machine precision. + offset = np.median((field - ref)[interior]) + residual = np.abs((field - ref) - offset)[interior] + assert residual.max() < 1e-6 + + +def test_mesh_matches_pygeodesic(): + geodesic = pytest.importorskip("pygeodesic.geodesic") + verts, faces = sphere_mesh(n_points=400) + ours = bic.distance.geodesic_distance_field_mesh(verts, faces, np.array([0], np.int64)) + algo = geodesic.PyGeodesicAlgorithmExact(verts, np.ascontiguousarray(faces, np.int32)) + ref, _ = algo.geodesicDistances(np.array([0], np.int32), None) + ref = np.asarray(ref) + reachable = ref > 1e-9 + rel = np.abs(ours[reachable] - ref[reachable]) / ref[reachable] + # first-order FMM vs exact MMP: modest error, no obtuse unfolding yet + assert np.mean(rel) < 0.06 + assert np.percentile(rel, 95) < 0.12