feat(core): add point-search, boundary-normal and vertebra-angle helpers

TPTBox/core/np_utils.py

@@ -20,10 +20,13 @@
 from scipy.ndimage import (
     binary_erosion,
     center_of_mass,
+    convolve,
     distance_transform_edt,
     gaussian_filter,
     generate_binary_structure,
+    sobel,
 )
+from scipy.spatial.distance import cdist
 from skimage.measure import euler_number as _euler_number
 from skimage.measure import label as _label
 
@@ -893,6 +896,131 @@ def np_point_coordinates(
     return surface_points
 
 
+def np_unit_vector(vector: np.ndarray) -> np.ndarray:
+    """Returns the unit vector of the input vector.
+
+    Args:
+        vector (np.ndarray): Any non-zero numeric array.
+
+    Returns:
+        np.ndarray: Array with the same direction as ``vector`` but unit length.
+    """
+    return vector / np.linalg.norm(vector)
+
+
+def np_angle_between(v1, v2, degrees: bool = False) -> float:
+    """Calculates the angle between two vectors.
+
+    Args:
+        v1: The first vector.
+        v2: The second vector.
+        degrees (bool, optional): Return the angle in degrees instead of radians. Defaults to False.
+
+    Returns:
+        float: The angle between ``v1`` and ``v2``.
+
+    Examples:
+        >>> np_angle_between((1, 0, 0), (0, 1, 0))
+        1.5707963267948966
+        >>> np_angle_between((1, 0, 0), (0, 1, 0), degrees=True)
+        90.0
+        >>> np_angle_between((1, 0, 0), (-1, 0, 0))
+        3.141592653589793
+    """
+    rad = np.arccos(np.clip(np.dot(np_unit_vector(v1), np_unit_vector(v2)), -1.0, 1.0))
+    return float(np.degrees(rad)) if degrees else float(rad)
+
+
+def np_index(arr: np.ndarray, entry) -> np.ndarray:
+    """Finds the rows of a point array that exactly equal a given point.
+
+    Args:
+        arr (np.ndarray): Array of shape ``(N, D)`` holding ``N`` points.
+        entry: The point to look for, broadcastable to shape ``(D,)``.
+
+    Returns:
+        np.ndarray: 1-D integer array with the indices of every matching row.
+        Empty if the point is not present.
+    """
+    arr = np.asarray(arr)
+    assert arr.ndim == 2, f"expected a (N, D) point array, got shape {arr.shape}"
+    return np.flatnonzero((arr == np.asarray(entry)).all(axis=1))
+
+
+def np_find_closest_point_index(point_arr: np.ndarray, point) -> int:
+    """Finds the index of the point in ``point_arr`` closest to ``point``.
+
+    For integer point arrays (i.e. voxel coordinates) an exact match is tried first, which is
+    the common case when the query already sits on the grid; only if that fails is the full
+    distance search run. ``scipy``'s compiled ``cdist`` is used for the search -- a hand-rolled
+    ``einsum`` over ``point_arr - point`` is roughly 2-6x slower because it materialises the
+    difference array.
+
+    Args:
+        point_arr (np.ndarray): Array of shape ``(N, D)`` holding the candidate points.
+        point: The query point, broadcastable to shape ``(D,)``.
+
+    Returns:
+        int: Index into ``point_arr`` of the nearest point. When several points are equally
+        close, the lowest index is returned.
+    """
+    arr = np.asarray(point_arr)
+    assert arr.ndim == 2, f"expected a (N, D) point array, got shape {arr.shape}"
+    assert len(arr) != 0, "cannot search an empty point array"
+    p = np.asarray(point)
+    if np.issubdtype(arr.dtype, np.integer):
+        # Only meaningful for grids; rounding a float query would silently snap it to integers.
+        matches = np_index(arr, np.round(p).astype(arr.dtype))
+        if len(matches) != 0:
+            return int(matches[0])
+    return int(np.argmin(cdist([p], arr)[0]))
+
+
+def np_compute_boundary_normals(
+    arr: UINTARRAY,
+    label: int | Sequence[int],
+    other_label: int | Sequence[int],
+    sigma: float = 1.0,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Computes per-voxel surface normals where one label touches another.
+
+    Voxels of ``label`` that share a face with ``other_label`` form the interface. The normal at
+    each of them is the gradient of a Gaussian-smoothed ``label`` mask, so it points *into*
+    ``label`` (the direction of increasing mask density); negate it to point outwards.
+
+    Unlike :func:`TPTBox.core.poi_fun.ray_casting.calculate_pca_normal_np`, which returns a single
+    principal axis for a whole segmentation, this returns a normal per interface voxel.
+
+    Args:
+        arr (UINTARRAY): 3-dimensional label array.
+        label (int | Sequence[int]): Label(s) whose interface voxels are returned.
+        other_label (int | Sequence[int]): Label(s) that ``label`` must touch to count as interface.
+        sigma (float, optional): Standard deviation of the Gaussian applied before the gradient.
+            Larger values give smoother, less noisy normals. Defaults to 1.0.
+
+    Returns:
+        tuple[np.ndarray, np.ndarray]: ``(coords, normals)``, both of shape ``(N, 3)``. ``coords``
+        are the integer voxel coordinates of the interface, ``normals`` the matching unit vectors.
+        Both are empty when the two labels do not touch.
+    """
+    assert arr.ndim == 3, arr.ndim
+    mask = np_isin(arr, label)
+    other = np_isin(arr, other_label)
+    # 6-connectivity kernel: the six face neighbours, so a purely diagonal contact does not count.
+    kernel = np.zeros((3, 3, 3), dtype=np.uint8)
+    kernel[1, 1, 0] = kernel[1, 1, 2] = 1
+    kernel[1, 0, 1] = kernel[1, 2, 1] = 1
+    kernel[0, 1, 1] = kernel[2, 1, 1] = 1
+    touching = mask & (convolve(other.astype(np.uint8), kernel, mode="constant") > 0)
+    coords = np.argwhere(touching)
+    if len(coords) == 0:
+        return coords, np.zeros((0, 3), dtype=float)
+    smoothed = gaussian_filter(mask.astype(float), sigma=sigma)
+    gradient = np.stack([sobel(smoothed, axis=a) for a in range(3)], axis=-1)
+    gradient /= np.linalg.norm(gradient, axis=-1, keepdims=True) + 1e-8
+    return coords, gradient[touching]
+
+
 def np_connected_components(
     arr: UINTARRAY,
     label_ref: LABEL_REFERENCE | None = None,

TPTBox/core/poi_fun/pixel_based_point_finder.py

@@ -7,6 +7,7 @@
 from scipy.spatial.distance import cdist
 
 from TPTBox import NII, POI, Logger_Interface, Print_Logger
+from TPTBox.core.np_utils import np_find_closest_point_index
 from TPTBox.core.poi_fun._help import to_local_np
 from TPTBox.core.poi_fun.vertebra_direction import _get_sub_array_by_direction, get_direction, get_vert_direction_matrix
 from TPTBox.core.vert_constants import COORDINATE, DIRECTIONS, Location
@@ -32,13 +33,9 @@ def get_nearest_neighbor(
         1-D integer array ``(x, y, z)`` of the voxel in ``sr_msk`` with label
         ``region_label`` that minimises the Euclidean distance to ``p``.
     """
-    if len(p.shape) == 1:
-        p = np.expand_dims(p, 1)
-    locs = np.where(sr_msk == region_label)
-    locs_array = np.array(list(locs)).T
-    distances = cdist(p.T, locs_array)
-
-    return locs_array[distances.argmin()]
+    p = np.asarray(p).reshape(-1)
+    locs_array = np.array(list(np.where(sr_msk == region_label))).T
+    return locs_array[np_find_closest_point_index(locs_array, p)]
 
 
 def max_distance_ray_cast_pixel_level(
@@ -248,9 +245,7 @@ def project_pois_onto_set_of_points(poi: POI, point_set: list[COORDINATE]) -> PO
     point_arr = np.asarray(point_set)
 
     for r, s, c in poi.items():
-        distance_to_point = cdist_to_point(c, point_arr)
-        new_coord = point_arr[np.argmin(distance_to_point)]
-        poi_n[r, s] = new_coord
+        poi_n[r, s] = point_arr[np_find_closest_point_index(point_arr, c)]
 
     return poi_n
 

TPTBox/core/poi_fun/ray_casting.py

@@ -8,6 +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.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
@@ -16,9 +17,7 @@
 _log = Print_Logger()
 
 
-def unit_vector(vector: np.ndarray) -> np.ndarray:
-    """Returns the unit vector of the vector."""
-    return vector / np.linalg.norm(vector)
+unit_vector = np_unit_vector
 
 
 # @njit(fastmath=True)

TPTBox/core/poi_fun/vertebra_direction.py

@@ -8,6 +8,7 @@
 from numpy.linalg import norm
 
 from TPTBox import NII, POI, Print_Logger, calc_poi_from_subreg_vert
+from TPTBox.core.np_utils import np_angle_between
 from TPTBox.core.poi_fun._help import make_spine_plot, sacrum_w_o_direction
 from TPTBox.core.vert_constants import DIRECTIONS, Location, Vertebra_Instance, _plane_dict, never_called
 
@@ -409,6 +410,32 @@ def get_vert_direction_matrix(poi: POI, vert_id: int, to_pir: bool = False) -> t
     return to_vert_orient, from_vert_orient
 
 
+def get_vert_direction_angles(poi: POI, vert_id: int, to_pir: bool = False, degrees: bool = True) -> tuple[float, float, float]:
+    """Return how far a vertebra's local frame is tilted from the global PIR axes.
+
+    Each of the vertebra's Posterior/Inferior/Right direction vectors is compared with the
+    corresponding global axis, giving one angle per axis. Useful as a quality check: a healthy
+    frame stays within a few degrees of orthogonal to its neighbours, whereas an implausibly
+    large angle (e.g. a "posterior" direction more than 90 degrees from global posterior)
+    indicates that the direction landmarks are wrong.
+
+    Args:
+        poi: ``POI`` object with pre-computed vertebra direction landmarks.
+        vert_id: Vertebra identifier (integer label).
+        to_pir: Whether to convert the POI to isotropic PIR space before computing.
+            Defaults to ``False``.
+        degrees: Return the angles in degrees rather than radians. Defaults to ``True``.
+
+    Returns:
+        Tuple of three angles ``(posterior, inferior, right)`` between the vertebra's direction
+        vectors and the global PIR axes.
+    """
+    directions = get_vert_direction_PIR(poi, vert_id=vert_id, to_pir=to_pir)
+    global_pir = (np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1]))
+    a, b, c = (np_angle_between(v, g, degrees=degrees) for v, g in zip(directions, global_pir))
+    return a, b, c
+
+
 def calc_center_spinal_cord(
     poi: POI,
     subreg: NII,

TPTBox/spine/spinestats/angles.py

@@ -10,6 +10,7 @@
 from TPTBox import POI, Image_Reference
 from TPTBox.core.compat import zip_strict
 from TPTBox.core.nii_wrapper import to_nii
+from TPTBox.core.np_utils import np_angle_between, np_unit_vector
 from TPTBox.core.vert_constants import DIRECTIONS, Location, Vertebra_Instance
 from TPTBox.spine.snapshot2D.snapshot_modular import Snapshot_Frame, create_snapshot
 
@@ -158,39 +159,9 @@ def get_stop_vert(self, poi) -> Vertebra_Instance:
 }
 
 
-def unit_vector(vector: np.ndarray) -> np.ndarray:
-    """Return the unit vector of the input vector.
-
-    Args:
-        vector: Any non-zero numeric array.
-
-    Returns:
-        Array with the same direction as ``vector`` but unit length.
-    """
-    return vector / np.linalg.norm(vector)
-
-
-def angle_between(v1, v2) -> float:
-    """Calculates the angle in radians between two vectors.
-
-    Args:
-        v1 (tuple): The first vector.
-        v2 (tuple): The second vector.
-
-    Returns:
-        float: The angle in radians between vectors 'v1' and 'v2'.
-
-    Examples:
-        >>> angle_between((1, 0, 0), (0, 1, 0))
-        1.5707963267948966
-        >>> angle_between((1, 0, 0), (1, 0, 0))
-        0.0
-        >>> angle_between((1, 0, 0), (-1, 0, 0))
-        3.141592653589793
-    """
-    v1_u = unit_vector(v1)
-    v2_u = unit_vector(v2)
-    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
+# Canonical implementations live in np_utils; re-exported here under their historic names.
+unit_vector = np_unit_vector
+angle_between = np_angle_between
 
 
 def get_to_space(a, b, c) -> tuple[np.ndarray, np.ndarray]:

TPTBox/spine/spinestats/poi_fun/endplates.py

@@ -8,6 +8,7 @@
 from stl.mesh import Mesh
 
 from TPTBox import NII, POI, Location, Logger_Interface, Print_Logger
+from TPTBox.core.poi_fun.pixel_based_point_finder import cdist_to_point
 from TPTBox.core.vert_constants import Vertebra_Instance
 
 _log = Print_Logger()
@@ -67,9 +68,7 @@ def _ray_cast_to_mesh(mesh: Mesh | trimesh.Trimesh, origin: np.ndarray, directio
         if len(locations) == 0:
             return None
         # closest hit
-        # d = np.linalg.norm(locations - origin, axis=1)
-        # return locations[np.argmin(d)]
-        d = np.linalg.norm(locations - origin, axis=1)
+        d = cdist_to_point(origin, locations)
         return origin + d.mean() * direction
     # numpy-stl fallback
     ts = []
@@ -133,7 +132,7 @@ def _local_curvature_grid(
     """
     world_pts = poi.local_to_global_arr(voxel_pts_full)
     casted_world = np.asarray(poi.local_to_global(tuple(casted_full)), dtype=float)
-    dists = np.linalg.norm(world_pts - casted_world, axis=1)
+    dists = cdist_to_point(casted_world, world_pts)
     neighborhood = world_pts[dists <= radius]
     if len(neighborhood) < 3:
         return 0.0