feat(np_utils): add non-convex ray marching and label interface thickness

TPTBox/core/nii_wrapper.py

@@ -49,6 +49,7 @@
     np_get_connected_components_center_of_mass,
     np_is_empty,
     np_isin,
+    np_label_interface_thickness,
     np_map_labels,
     np_map_labels_based_on_majority_label_mask_overlap,
     np_point_coordinates,
@@ -3055,6 +3056,45 @@ def filter_connected_components_by_bbox_chain(
         )
         return self.set_array(arr, inplace=inplace)
 
+    def label_interface_thickness(
+        self,
+        label: int | Sequence[int],
+        other_label: int | Sequence[int],
+        max_count_component: int | None = None,
+        sigma: float = 1.0,
+        max_steps: int | None = 1000,
+        max_distance: float | None = None,
+    ) -> np.ndarray:
+        """Measures how thick a structure is where it meets another structure.
+
+        At every voxel of ``other_label`` touching ``label``, a ray is marched along the inward
+        surface normal until it exits ``label``; the distance travelled is the local thickness.
+        Distances are returned in **millimetres**, using :attr:`zoom`.
+
+        Args:
+            label (int | Sequence[int]): The structure whose thickness is measured.
+            other_label (int | Sequence[int]): The structure the measurement starts from.
+            max_count_component (int | None, optional): Keep only this many largest connected
+                components of ``label`` before measuring. Defaults to None (keep all).
+            sigma (float, optional): Smoothing applied before computing the normals. Defaults to 1.0.
+            max_steps (int | None, optional): Step limit per ray. Defaults to 1000.
+            max_distance (float | None, optional): Distance limit per ray, in voxels. Defaults to None.
+
+        Returns:
+            np.ndarray: One thickness in mm per interface voxel; ``np.nan`` where a ray hit a
+            limit without leaving ``label``. Empty when the two labels do not touch.
+        """
+        return np_label_interface_thickness(
+            self.get_seg_array(),
+            label,
+            other_label,
+            zoom=self.zoom,
+            max_count_component=max_count_component,
+            sigma=sigma,
+            max_steps=max_steps,
+            max_distance=max_distance,
+        )
+
 
 
 

TPTBox/core/np_utils.py

@@ -1064,6 +1064,145 @@ def np_find_closest_point_index(point_arr: np.ndarray, point) -> int:
     return int(np.argmin(cdist([p], arr)[0]))
 
 
+def np_raymarch_until_background(
+    arr: np.ndarray,
+    start_coord: Sequence[float] | np.ndarray,
+    direction_vector: np.ndarray,
+    step_size: float = 0.0625,
+    max_steps: int | None = 1000,
+    max_distance: float | None = None,
+    threshold: float = 0.5,
+    interpolator=None,
+) -> np.ndarray | None:
+    """Marches a ray from a point until it leaves a mask.
+
+    Takes fixed-size steps, sampling the mask with trilinear interpolation at each one, and stops
+    as soon as the interpolated value drops below ``threshold`` or the ray leaves the volume.
+
+    Unlike :func:`TPTBox.core.poi_fun.ray_casting.max_distance_ray_cast_convex`, which bisects and
+    therefore assumes the region is convex, this walks the ray step by step and so handles
+    concave and multi-lobed structures correctly -- at the cost of being slower.
+
+    Args:
+        arr (np.ndarray): The mask to march through. Non-zero is inside.
+        start_coord (Sequence[float] | np.ndarray): Voxel coordinate the ray starts at.
+        direction_vector (np.ndarray): Direction of the ray; normalised internally.
+        step_size (float, optional): Step length in voxels. Defaults to 0.0625 (1/16 of a voxel).
+        max_steps (int | None, optional): Give up after this many steps. Defaults to 1000.
+        max_distance (float | None, optional): Give up once this distance in voxels is covered.
+            Defaults to None. At least one of ``max_steps``/``max_distance`` must be set.
+        threshold (float, optional): Interpolated value below which the ray is considered to have
+            left the mask. Defaults to 0.5.
+        interpolator (RegularGridInterpolator | None, optional): Prebuilt interpolator over
+            ``arr``. Pass one when marching many rays through the same mask -- building it is
+            far more expensive than the march itself. Defaults to None (built internally).
+
+    Returns:
+        np.ndarray | None: The coordinate where the ray left the mask, or None if it was still
+        inside when the step/distance limit was reached.
+    """
+    assert max_steps is not None or max_distance is not None, "at least one of max_steps or max_distance must be set"
+    from scipy.interpolate import RegularGridInterpolator
+
+    start = np.asarray(start_coord, dtype=np.float64)
+    direction = np.asarray(direction_vector, dtype=np.float64)
+    direction = direction / (np.linalg.norm(direction) + 1e-10)
+    step_vector = direction * step_size
+
+    if interpolator is None:
+        interpolator = RegularGridInterpolator(
+            [np.arange(s, dtype=np.float64) for s in arr.shape], arr.astype(np.float64), bounds_error=False, fill_value=0.0
+        )
+
+    pos = start.copy()
+    steps = int(max_steps) if max_steps is not None else int(1e8)
+    for _ in range(steps):
+        if max_distance is not None and np.linalg.norm(pos - start) >= max_distance:
+            return None
+        if np.any(pos < 0) or np.any(pos >= arr.shape):
+            return pos
+        if interpolator(pos) < threshold:
+            return pos
+        pos = pos + step_vector
+    return None
+
+
+def np_label_interface_thickness(
+    arr: UINTARRAY,
+    label: int | Sequence[int],
+    other_label: int | Sequence[int],
+    zoom: Sequence[float] | None = None,
+    max_count_component: int | None = None,
+    sigma: float = 1.0,
+    step_size: float | None = None,
+    max_steps: int | None = 1000,
+    max_distance: float | None = None,
+) -> np.ndarray:
+    """Measures how thick a structure is where it meets another structure.
+
+    At every voxel of ``other_label`` that touches ``label``, a ray is marched from that voxel
+    along the inward surface normal until it exits ``label``. The distance travelled is the local
+    thickness of ``label`` at that point on the interface.
+
+    Args:
+        arr (UINTARRAY): 3-dimensional label array.
+        label (int | Sequence[int]): The structure whose thickness is measured.
+        other_label (int | Sequence[int]): The structure the measurement starts from.
+        zoom (Sequence[float] | None, optional): Voxel spacing in mm. When given, distances are
+            returned in millimetres and the default step size is derived from it; otherwise
+            distances are in voxels. Defaults to None.
+        max_count_component (int | None, optional): Keep only this many largest connected
+            components of ``label`` before measuring. Defaults to None (keep all).
+        sigma (float, optional): Smoothing applied before computing the normals. Defaults to 1.0.
+        step_size (float | None, optional): March step in voxels. Defaults to ``min(zoom)/16``
+            when ``zoom`` is given, else 1/16 of a voxel.
+        max_steps (int | None, optional): Step limit per ray. Defaults to 1000.
+        max_distance (float | None, optional): Distance limit per ray, in voxels. Defaults to None.
+
+    Returns:
+        np.ndarray: 1-D array with one thickness per interface voxel. Rays that hit a limit
+        without leaving ``label`` contribute ``np.nan``. Empty when the labels do not touch.
+    """
+    from scipy.interpolate import RegularGridInterpolator
+
+    assert arr.ndim == 3, arr.ndim
+    if step_size is None:
+        step_size = (min(zoom) / 16) if zoom is not None else 1 / 16
+
+    work = np.where(np_isin(arr, label) | np_isin(arr, other_label), arr, 0)
+    if max_count_component is not None:
+        kept = np_filter_connected_components(np_isin(work, label), largest_k_components=max_count_component, connectivity=1)
+        work = np.where(np_isin(work, other_label) | (kept != 0), work, 0)
+
+    coords, normals = np_compute_boundary_normals(work, other_label, label, sigma=sigma)
+    if len(coords) == 0:
+        return np.zeros(0, dtype=float)
+
+    # March through the measured structure plus the interface voxels themselves, so a ray starting
+    # on the interface begins inside the mask rather than immediately outside it.
+    march_mask = np_isin(work, label).astype(np.float64)
+    march_mask[tuple(coords.T)] = 1.0
+    interpolator = RegularGridInterpolator(
+        [np.arange(s, dtype=np.float64) for s in march_mask.shape], march_mask, bounds_error=False, fill_value=0.0
+    )
+
+    spacing = np.asarray(zoom, dtype=float) if zoom is not None else np.ones(3)
+    out = np.empty(len(coords), dtype=float)
+    for i, (coord, normal) in enumerate(zip(coords, normals)):
+        # np_compute_boundary_normals points into `other_label`; flip it to march into `label`.
+        end = np_raymarch_until_background(
+            march_mask,
+            coord,
+            -normal,
+            step_size=step_size,
+            max_steps=max_steps,
+            max_distance=max_distance,
+            interpolator=interpolator,
+        )
+        out[i] = np.nan if end is None else float(np.linalg.norm((end - coord) * spacing))
+    return out
+
+
 def np_compute_boundary_normals(
     arr: UINTARRAY,
     label: int | Sequence[int],

TPTBox/core/poi_fun/ray_casting.py

@@ -8,7 +8,7 @@
 from sklearn.decomposition import PCA
 
 from TPTBox import NII, POI, Print_Logger, Vertebra_Instance
-from TPTBox.core.np_utils import np_unit_vector
+from TPTBox.core.np_utils import np_raymarch_until_background, np_unit_vector
 from TPTBox.core.poi_fun._help import sacrum_w_o_arcus, to_local_np
 from TPTBox.core.poi_fun.pixel_based_point_finder import get_direction
 from TPTBox.core.vert_constants import COORDINATE, DIRECTIONS, Location
@@ -183,6 +183,49 @@ def is_inside(distance):
     return start_point_np + normal_vector * ((min_v + max_v) / 2)
 
 
+def max_distance_ray_cast_non_convex(
+    region: NII,
+    start_coord: COORDINATE | np.ndarray,
+    direction_vector: np.ndarray,
+    step_size: float | None = None,
+    max_steps: int | None = 1000,
+    max_distance: float | None = None,
+    threshold: float = 0.5,
+) -> np.ndarray | None:
+    """Find the exit point of a ray inside an arbitrary NII region.
+
+    The non-convex counterpart to :func:`max_distance_ray_cast_convex`. That one bisects, which
+    silently skips interior gaps -- on a mask with a hole it reports the far outer wall rather
+    than the near edge of the hole. This walks the ray in fixed steps instead, so concave and
+    multi-lobed regions are handled correctly, at the cost of being slower.
+
+    Args:
+        region: ``NII`` object whose nonzero voxels define the region.
+        start_coord: Starting coordinate ``(x, y, z)`` of the ray in voxel space.
+        direction_vector: Direction of the ray; normalised internally.
+        step_size: Step length in voxels. Defaults to ``min(region.zoom) / 16``.
+        max_steps: Give up after this many steps. Defaults to 1000.
+        max_distance: Give up once this distance in voxels is covered. Defaults to None.
+        threshold: Interpolated mask value below which the ray has left the region.
+            Defaults to 0.5.
+
+    Returns:
+        3-element numpy array with the exit coordinate, or ``None`` if the ray was still inside
+        the region when the step or distance limit was reached.
+    """
+    if step_size is None:
+        step_size = min(region.zoom) / 16
+    return np_raymarch_until_background(
+        region.get_array(),
+        start_coord,
+        direction_vector,
+        step_size=step_size,
+        max_steps=max_steps,
+        max_distance=max_distance,
+        threshold=threshold,
+    )
+
+
 def max_distance_ray_cast_convex(
     region: NII,
     start_coord: COORDINATE | np.ndarray,

unit_tests/test_nii_extended.py

@@ -586,5 +586,25 @@ def test_is_symmetric(self):
             self.assertEqual(np_slices_overlap(a, b), np_slices_overlap(b, a))
 
 
+class Test_label_interface_thickness_nii(unittest.TestCase):
+    def _two_slabs(self, zoom=(1.0, 1.0, 1.0)):
+        seg = np.zeros((40, 20, 20), dtype=np.uint8)
+        seg[10:18] = 1
+        seg[18:22] = 2
+        return _make_nii(seg, zoom=zoom)
+
+    def test_returns_millimetres_using_zoom(self):
+        one = self._two_slabs().label_interface_thickness(1, 2)
+        two = self._two_slabs(zoom=(2.0, 1.0, 1.0)).label_interface_thickness(1, 2)
+        self.assertAlmostEqual(float(np.nanmean(two)) / float(np.nanmean(one)), 2.0, places=5)
+
+    def test_matches_the_array_level_function(self):
+        from TPTBox.core.np_utils import np_label_interface_thickness
+
+        nii = self._two_slabs(zoom=(1.5, 1.0, 1.0))
+        expected = np_label_interface_thickness(nii.get_seg_array(), 1, 2, zoom=nii.zoom)
+        self.assertTrue(np.allclose(nii.label_interface_thickness(1, 2), expected, equal_nan=True))
+
+
 if __name__ == "__main__":
     unittest.main()

unit_tests/test_nputils.py

@@ -661,6 +661,92 @@ def test_inputs_are_not_mutated(self):
         self.assertTrue(np.array_equal(b, b_before))
 
 
+class Test_raymarch_until_background(unittest.TestCase):
+    def test_exit_point_on_a_slab(self):
+        arr = np.zeros((40, 10, 10))
+        arr[10:25] = 1
+        end = np_utils.np_raymarch_until_background(arr, [10, 5, 5], [1, 0, 0])
+        self.assertIsNotNone(end)
+        # the 0.5 crossing sits at the slab's far edge, within one step
+        self.assertAlmostEqual(float(end[0]), 24.5, delta=0.1)
+
+    def test_stops_at_an_interior_gap(self):
+        """The whole point of the non-convex marcher: a bisection search skips this gap."""
+        arr = np.zeros((40, 10, 10))
+        arr[5:35] = 1
+        arr[15:20] = 0
+        end = np_utils.np_raymarch_until_background(arr, [6, 5, 5], [1, 0, 0])
+        self.assertAlmostEqual(float(end[0]), 14.5, delta=0.1)
+
+    def test_convex_bisection_disagrees_on_the_same_mask(self):
+        import nibabel as nib
+
+        from TPTBox import NII
+        from TPTBox.core.poi_fun.ray_casting import max_distance_ray_cast_convex
+
+        arr = np.zeros((40, 10, 10), dtype=np.uint8)
+        arr[5:35] = 1
+        arr[15:20] = 0
+        nii = NII(nib.Nifti1Image(arr, np.eye(4)), seg=True)
+        convex = max_distance_ray_cast_convex(nii, np.array([6, 5, 5]), np.array([1.0, 0.0, 0.0]))
+        self.assertGreater(float(convex[0]), 30.0)  # skipped the gap entirely
+
+    def test_returns_none_when_it_never_leaves(self):
+        arr = np.ones((40, 10, 10))
+        self.assertIsNone(np_utils.np_raymarch_until_background(arr, [1, 5, 5], [1, 0, 0], max_steps=5))
+
+    def test_max_distance_limit(self):
+        arr = np.ones((40, 10, 10))
+        self.assertIsNone(np_utils.np_raymarch_until_background(arr, [1, 5, 5], [1, 0, 0], max_steps=None, max_distance=2.0))
+
+    def test_requires_a_limit(self):
+        with self.assertRaises(AssertionError):
+            np_utils.np_raymarch_until_background(np.ones((5, 5, 5)), [1, 1, 1], [1, 0, 0], max_steps=None)
+
+    def test_leaving_the_volume_returns_a_position(self):
+        arr = np.ones((10, 10, 10))
+        end = np_utils.np_raymarch_until_background(arr, [8, 5, 5], [1, 0, 0])
+        self.assertIsNotNone(end)
+        self.assertGreaterEqual(float(end[0]), 9.0)
+
+
+class Test_label_interface_thickness(unittest.TestCase):
+    def _two_slabs(self):
+        seg = np.zeros((40, 20, 20), dtype=np.uint8)
+        seg[10:18] = 1  # measured structure
+        seg[18:22] = 2  # reference structure
+        return seg
+
+    def test_thickness_of_a_known_slab(self):
+        t = np_utils.np_label_interface_thickness(self._two_slabs(), label=1, other_label=2)
+        self.assertEqual(len(t), 20 * 20)
+        # measured from the interface voxel center to the 0.5 isosurface
+        self.assertAlmostEqual(float(np.nanmean(t)), 8.5, delta=0.15)
+        self.assertAlmostEqual(float(np.nanstd(t)), 0.0, delta=1e-6)
+
+    def test_zoom_scales_the_result(self):
+        seg = self._two_slabs()
+        voxels = np_utils.np_label_interface_thickness(seg, 1, 2)
+        mm = np_utils.np_label_interface_thickness(seg, 1, 2, zoom=(2.0, 1.0, 1.0))
+        self.assertAlmostEqual(float(np.nanmean(mm)) / float(np.nanmean(voxels)), 2.0, places=5)
+
+    def test_no_contact_returns_empty(self):
+        seg = np.zeros((20, 20, 20), dtype=np.uint8)
+        seg[2:5] = 1
+        seg[15:18] = 2
+        self.assertEqual(np_utils.np_label_interface_thickness(seg, 1, 2).shape, (0,))
+
+    def test_max_count_component_drops_smaller_components(self):
+        seg = np.zeros((40, 20, 20), dtype=np.uint8)
+        seg[10:18, 0:12, :] = 1  # large component
+        seg[10:18, 15:18, :] = 1  # smaller, separate component
+        seg[18:22] = 2
+        every = np_utils.np_label_interface_thickness(seg, 1, 2)
+        largest = np_utils.np_label_interface_thickness(seg, 1, 2, max_count_component=1)
+        self.assertLess(len(largest), len(every))
+        self.assertEqual(len(np_utils.np_label_interface_thickness(seg, 1, 2, max_count_component=2)), len(every))
+
+
 if __name__ == "__main__":
     unittest.main()