|
| 1 | +"""Reference geodesic-distance implementations for development cross-checks. |
| 2 | +
|
| 3 | +These are the oracles that ``bic.distance`` geodesic functions are validated |
| 4 | +against. There is one backend per geometry, because scikit-fmm only supports |
| 5 | +regular Cartesian grids: |
| 6 | +
|
| 7 | +- **masks** -> scikit-fmm (``skfmm``): fast marching / Eikonal solver. |
| 8 | +- **meshes** -> pygeodesic: exact Mitchell-Mount-Papadimitriou (MMP) geodesics. |
| 9 | +
|
| 10 | +Not part of the pytest suite; requires ``scikit-fmm`` and/or ``pygeodesic``. |
| 11 | +Importing this module is cheap (numpy only); each function imports its backend |
| 12 | +lazily and raises a clear ``ImportError`` if the backend is missing. |
| 13 | +
|
| 14 | +Note on the scikit-fmm point-source idiom: to get the distance to a set of |
| 15 | +seed voxels, we set ``phi = +1`` everywhere and ``phi = -1`` at the seeds and |
| 16 | +take ``|skfmm.distance(phi)|``. The zero contour then sits roughly half a cell |
| 17 | +away from each seed, so the reference is offset from a "distance is exactly 0 at |
| 18 | +the seed" convention by ~0.5 * dx near the sources. Compare with a tolerance, |
| 19 | +or exclude the immediate neighbourhood of the seeds, when checking against an |
| 20 | +implementation that initialises the seeds to exactly 0. |
| 21 | +""" |
| 22 | + |
| 23 | +from __future__ import annotations |
| 24 | + |
| 25 | +import numpy as np |
| 26 | + |
| 27 | + |
| 28 | +# --------------------------------------------------------------------------- # |
| 29 | +# masks: scikit-fmm |
| 30 | +# --------------------------------------------------------------------------- # |
| 31 | + |
| 32 | + |
| 33 | +def _skfmm(): |
| 34 | + try: |
| 35 | + import skfmm |
| 36 | + except ImportError as error: # pragma: no cover - dev script |
| 37 | + raise ImportError( |
| 38 | + "scikit-fmm is required for the mask geodesic reference " |
| 39 | + "(`pip install scikit-fmm`)" |
| 40 | + ) from error |
| 41 | + return skfmm |
| 42 | + |
| 43 | + |
| 44 | +def _normalize_dx(sampling, ndim): |
| 45 | + if sampling is None: |
| 46 | + return 1.0 |
| 47 | + if np.isscalar(sampling): |
| 48 | + return float(sampling) |
| 49 | + dx = [float(s) for s in sampling] |
| 50 | + if len(dx) != ndim: |
| 51 | + raise ValueError(f"sampling must have length {ndim}, got {len(dx)}") |
| 52 | + return dx |
| 53 | + |
| 54 | + |
| 55 | +def reference_geodesic_field_mask(mask, sources, sampling=None, speed=None): |
| 56 | + """Geodesic distance field within a mask from a set of source coordinates. |
| 57 | +
|
| 58 | + Mirrors ``bic.distance.geodesic_distance_field``: for every voxel, the |
| 59 | + shortest-path distance to the nearest source, constrained to the nonzero |
| 60 | + region of ``mask``. Voxels outside the mask (and voxels unreachable from any |
| 61 | + source) are returned as ``+inf``. |
| 62 | + """ |
| 63 | + skfmm = _skfmm() |
| 64 | + mask = np.ascontiguousarray(mask) != 0 |
| 65 | + sources = np.atleast_2d(np.asarray(sources, dtype=np.int64)) |
| 66 | + dx = _normalize_dx(sampling, mask.ndim) |
| 67 | + |
| 68 | + phi = np.ones(mask.shape, dtype=np.float64) |
| 69 | + phi[tuple(sources.T)] = -1.0 |
| 70 | + phi = np.ma.MaskedArray(phi, mask=~mask) |
| 71 | + |
| 72 | + # order=1 matches bioimage-cpp's first-order Godunov scheme (so agreement is |
| 73 | + # limited only by the ~0.5-cell seed offset documented above, not by the |
| 74 | + # stencil order). |
| 75 | + if speed is None: |
| 76 | + field = np.ma.abs(skfmm.distance(phi, dx=dx, order=1)) |
| 77 | + else: |
| 78 | + speed = np.ascontiguousarray(speed, dtype=np.float64) |
| 79 | + field = skfmm.travel_time(phi, speed, dx=dx, order=1) |
| 80 | + |
| 81 | + out = np.full(mask.shape, np.inf, dtype=np.float64) |
| 82 | + field_data = np.ma.getdata(field) |
| 83 | + valid = ~np.ma.getmaskarray(field) |
| 84 | + out[valid] = field_data[valid] |
| 85 | + return out |
| 86 | + |
| 87 | + |
| 88 | +def reference_geodesic_distances_mask(mask, points, sampling=None, speed=None): |
| 89 | + """Full pairwise geodesic distance matrix between points within a mask. |
| 90 | +
|
| 91 | + Runs the field solve once per point and reads the field at the other |
| 92 | + points. The result is symmetrized (per-source solves can differ by a tiny |
| 93 | + numerical amount) with a zero diagonal. |
| 94 | + """ |
| 95 | + points = np.atleast_2d(np.asarray(points, dtype=np.int64)) |
| 96 | + n = len(points) |
| 97 | + out = np.full((n, n), np.inf, dtype=np.float64) |
| 98 | + index = tuple(points.T) |
| 99 | + for i in range(n): |
| 100 | + field = reference_geodesic_field_mask(mask, points[i][None, :], sampling, speed) |
| 101 | + out[i] = field[index] |
| 102 | + out = 0.5 * (out + out.T) |
| 103 | + np.fill_diagonal(out, 0.0) |
| 104 | + return out |
| 105 | + |
| 106 | + |
| 107 | +# --------------------------------------------------------------------------- # |
| 108 | +# meshes: pygeodesic (exact MMP) |
| 109 | +# --------------------------------------------------------------------------- # |
| 110 | + |
| 111 | + |
| 112 | +def _pygeodesic(): |
| 113 | + try: |
| 114 | + import pygeodesic.geodesic as geodesic |
| 115 | + except ImportError as error: # pragma: no cover - dev script |
| 116 | + raise ImportError( |
| 117 | + "pygeodesic is required for the mesh geodesic reference " |
| 118 | + "(`pip install pygeodesic`)" |
| 119 | + ) from error |
| 120 | + return geodesic |
| 121 | + |
| 122 | + |
| 123 | +def _mesh_algorithm(vertices, faces): |
| 124 | + geodesic = _pygeodesic() |
| 125 | + vertices = np.ascontiguousarray(vertices, dtype=np.float64) |
| 126 | + faces = np.ascontiguousarray(faces, dtype=np.int32) |
| 127 | + return geodesic.PyGeodesicAlgorithmExact(vertices, faces) |
| 128 | + |
| 129 | + |
| 130 | +def reference_geodesic_field_mesh(vertices, faces, sources): |
| 131 | + """Geodesic distance field on a triangle mesh from a set of source vertices. |
| 132 | +
|
| 133 | + Mirrors ``bic.distance.geodesic_distance_field_mesh``. Returns a |
| 134 | + ``(n_vertices,)`` array of surface geodesic distance to the nearest source. |
| 135 | + """ |
| 136 | + algo = _mesh_algorithm(vertices, faces) |
| 137 | + sources = np.atleast_1d(np.asarray(sources, dtype=np.int32)) |
| 138 | + distances, _ = algo.geodesicDistances(sources, None) |
| 139 | + return np.asarray(distances, dtype=np.float64) |
| 140 | + |
| 141 | + |
| 142 | +def reference_geodesic_distances_mesh(vertices, faces, points): |
| 143 | + """Full pairwise geodesic distance matrix between mesh vertices.""" |
| 144 | + algo = _mesh_algorithm(vertices, faces) |
| 145 | + points = np.atleast_1d(np.asarray(points, dtype=np.int32)) |
| 146 | + n = len(points) |
| 147 | + out = np.full((n, n), np.inf, dtype=np.float64) |
| 148 | + for i in range(n): |
| 149 | + # Read the full single-source field and index the targets; passing an |
| 150 | + # explicit target array to pygeodesic can overflow on some meshes. |
| 151 | + distances, _ = algo.geodesicDistances(points[i : i + 1], None) |
| 152 | + out[i] = np.asarray(distances, dtype=np.float64)[points] |
| 153 | + out = 0.5 * (out + out.T) |
| 154 | + np.fill_diagonal(out, 0.0) |
| 155 | + return out |
0 commit comments