Implement scaffold for geodesic distance functions

Author: constantinpape

MIGRATION_GUIDE.md

@@ -2188,6 +2188,81 @@ 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`).
+
+> **Status:** the C++ solvers are not implemented yet. The functions below are
+> present and validate their arguments, but calling a well-formed input
+> currently raises `RuntimeError("... not yet implemented")`. This section
+> documents the intended API and behaviour.
+
+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
+
+# 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` |
+| 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.
+- Mesh geodesics are true surface (2-manifold) distances, validated against the
+  exact MMP algorithm in `pygeodesic`, not the grid fast-marching used for
+  masks. See `development/distance/` for the reference oracles.
+
 ## I/O and Build Dependencies
 
 `bioimage-cpp` intentionally does not replace nifty or affogato I/O helpers.

development/distance/_geodesic_reference.py

@@ -0,0 +1,150 @@
+"""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)
+
+    if speed is None:
+        field = np.ma.abs(skfmm.distance(phi, dx=dx))
+    else:
+        speed = np.ascontiguousarray(speed, dtype=np.float64)
+        field = skfmm.travel_time(phi, speed, dx=dx)
+
+    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):
+        distances, _ = algo.geodesicDistances(points[i : i + 1], points)
+        out[i] = np.asarray(distances, dtype=np.float64)
+    out = 0.5 * (out + out.T)
+    np.fill_diagonal(out, 0.0)
+    return out