POI Strategies¶

Internal modules that implement the various strategies for computing Points of Interest from segmentation volumes.

Ray Casting¶

TPTBox.core.poi_fun.ray_casting ¶

unit_vector ¶

unit_vector(vector: ndarray) -> np.ndarray

Returns the unit vector of the vector.

Source code in TPTBox/core/poi_fun/ray_casting.py

def unit_vector(vector: np.ndarray) -> np.ndarray:
    """Returns the unit vector of the vector."""
    return vector / np.linalg.norm(vector)

trilinear_interpolate ¶

trilinear_interpolate(volume: ndarray, x: float, y: float, z: float) -> float

Perform trilinear interpolation of a 3D volume at a sub-voxel coordinate.

Parameters:

Name	Type	Description	Default
`volume`	`ndarray`	3D array to interpolate.	required
`x`	`float`	Sub-voxel coordinate along the first axis.	required
`y`	`float`	Sub-voxel coordinate along the second axis.	required
`z`	`float`	Sub-voxel coordinate along the third axis.	required

Returns:

Type	Description
`float`	Interpolated scalar value at (x, y, z), or 0.0 if the point is outside
`float`	the valid interior of the volume.

Source code in TPTBox/core/poi_fun/ray_casting.py

def trilinear_interpolate(volume: np.ndarray, x: float, y: float, z: float) -> float:
    """Perform trilinear interpolation of a 3D volume at a sub-voxel coordinate.

    Args:
        volume: 3D array to interpolate.
        x: Sub-voxel coordinate along the first axis.
        y: Sub-voxel coordinate along the second axis.
        z: Sub-voxel coordinate along the third axis.

    Returns:
        Interpolated scalar value at (x, y, z), or 0.0 if the point is outside
        the valid interior of the volume.
    """
    xi, yi, zi = np.floor(x).astype(int), np.floor(y).astype(int), np.floor(z).astype(int)
    if xi < 0 or yi < 0 or zi < 0 or xi >= volume.shape[0] - 1 or yi >= volume.shape[1] - 1 or zi >= volume.shape[2] - 1:
        return 0.0

    xd, yd, zd = x - xi, y - yi, z - zi
    c000 = volume[xi, yi, zi]
    c100 = volume[xi + 1, yi, zi]
    c010 = volume[xi, yi + 1, zi]
    c110 = volume[xi + 1, yi + 1, zi]
    c001 = volume[xi, yi, zi + 1]
    c101 = volume[xi + 1, yi, zi + 1]
    c011 = volume[xi, yi + 1, zi + 1]
    c111 = volume[xi + 1, yi + 1, zi + 1]

    c00 = c000 * (1 - xd) + c100 * xd
    c01 = c001 * (1 - xd) + c101 * xd
    c10 = c010 * (1 - xd) + c110 * xd
    c11 = c011 * (1 - xd) + c111 * xd
    c0 = c00 * (1 - yd) + c10 * yd
    c1 = c01 * (1 - yd) + c11 * yd
    return c0 * (1 - zd) + c1 * zd

max_distance_ray_cast_convex_npfast ¶

max_distance_ray_cast_convex_npfast(region_array: ndarray, start_coord: ndarray, direction_vector: ndarray, acc_delta: float = 5e-05) -> np.ndarray

Find the exit point of a ray inside a convex region using bisection (fast path).

Uses trilinear interpolation and bisection search to locate the last point inside the region along the given direction. Prints a debug line at each bisection step (development helper — see max_distance_ray_cast_convex_np for the production version without debug output).

Parameters:

Name	Type	Description	Default
`region_array`	`ndarray`	Binary 3D numpy array where nonzero voxels define the region.	required
`start_coord`	`ndarray`	Starting coordinate `(x, y, z)` of the ray (must be inside the region).	required
`direction_vector`	`ndarray`	Direction of the ray; need not be a unit vector.	required
`acc_delta`	`float`	Bisection stops when the bracket width falls below this value. Defaults to 0.00005.	`5e-05`

Returns:

Type	Description
`ndarray`	3-element array with the approximate exit coordinate along the ray.
`ndarray`	If `start_coord` is already outside the region the start coordinate is
`ndarray`	returned unchanged.

Source code in TPTBox/core/poi_fun/ray_casting.py

def max_distance_ray_cast_convex_npfast(
    region_array: np.ndarray,
    start_coord: np.ndarray,
    direction_vector: np.ndarray,
    acc_delta: float = 0.00005,
) -> np.ndarray:
    """Find the exit point of a ray inside a convex region using bisection (fast path).

    Uses trilinear interpolation and bisection search to locate the last point
    inside the region along the given direction.  Prints a debug line at each
    bisection step (development helper — see ``max_distance_ray_cast_convex_np``
    for the production version without debug output).

    Args:
        region_array: Binary 3D numpy array where nonzero voxels define the region.
        start_coord: Starting coordinate ``(x, y, z)`` of the ray (must be inside
            the region).
        direction_vector: Direction of the ray; need not be a unit vector.
        acc_delta: Bisection stops when the bracket width falls below this value.
            Defaults to 0.00005.

    Returns:
        3-element array with the approximate exit coordinate along the ray.
        If ``start_coord`` is already outside the region the start coordinate is
        returned unchanged.
    """
    # Normalize direction
    norm_vec = direction_vector / np.sqrt((direction_vector**2).sum())

    # Quick exit if start point is outside
    if trilinear_interpolate(region_array, *start_coord) <= 0.5:
        return np.array(start_coord)

    min_v = 0.0
    max_v = np.sum(region_array.shape)
    delta = max_v - min_v

    while delta > acc_delta:
        mid = 0.5 * (max_v + min_v)
        x = start_coord[0] + norm_vec[0] * mid
        y = start_coord[1] + norm_vec[1] * mid
        z = start_coord[2] + norm_vec[2] * mid
        val = trilinear_interpolate(region_array, x, y, z)
        print(f"Raycast check at distance {mid:.2f}: value={val:.4f}")
        if val > 0.5:
            min_v = mid
        else:
            max_v = mid
        delta = max_v - min_v

    dist = 0.5 * (min_v + max_v)
    return np.array(
        [
            start_coord[0] + norm_vec[0] * dist,
            start_coord[1] + norm_vec[1] * dist,
            start_coord[2] + norm_vec[2] * dist,
        ]
    )

max_distance_ray_cast_convex_np ¶

max_distance_ray_cast_convex_np(region: ndarray, start_coord: COORDINATE | ndarray, direction_vector: ndarray, acc_delta: float = 5e-05, max_v: int | None = None) -> np.ndarray | None

Find the exit point of a ray inside a convex region (numpy array input).

Uses RegularGridInterpolator and bisection to locate the boundary of a convex binary region along the given direction vector.

Parameters:

Name	Type	Description	Default
`region`	`ndarray`	Binary 3D numpy array where values > 0.5 are considered inside.	required
`start_coord`	`COORDINATE \| ndarray`	Starting coordinate `(x, y, z)` of the ray. If already outside the region the start coordinate is returned unchanged.	required
`direction_vector`	`ndarray`	Direction of the ray; will be normalised internally.	required
`acc_delta`	`float`	Bisection convergence threshold in voxel units. Defaults to 0.00005.	`5e-05`
`max_v`	`int \| None`	Upper bound for the bisection search (in voxels along the ray). Defaults to `sum(region.shape)` when `None`.	`None`

Returns:

Type	Description
`ndarray \| None`	3-element numpy array with the approximate exit coordinate, or the start
`ndarray \| None`	coordinate if it lies outside the region. Returns `None` when
`ndarray \| None`	`start_coord` is `None`.

Source code in TPTBox/core/poi_fun/ray_casting.py

def max_distance_ray_cast_convex_np(
    region: np.ndarray,
    start_coord: COORDINATE | np.ndarray,
    direction_vector: np.ndarray,
    acc_delta: float = 0.00005,
    max_v: int | None = None,
) -> np.ndarray | None:
    """Find the exit point of a ray inside a convex region (numpy array input).

    Uses ``RegularGridInterpolator`` and bisection to locate the boundary of a
    convex binary region along the given direction vector.

    Args:
        region: Binary 3D numpy array where values > 0.5 are considered inside.
        start_coord: Starting coordinate ``(x, y, z)`` of the ray.  If already
            outside the region the start coordinate is returned unchanged.
        direction_vector: Direction of the ray; will be normalised internally.
        acc_delta: Bisection convergence threshold in voxel units.
            Defaults to 0.00005.
        max_v: Upper bound for the bisection search (in voxels along the ray).
            Defaults to ``sum(region.shape)`` when ``None``.

    Returns:
        3-element numpy array with the approximate exit coordinate, or the start
        coordinate if it lies outside the region.  Returns ``None`` when
        ``start_coord`` is ``None``.
    """
    start_point_np = np.asarray(start_coord)
    if start_point_np is None:
        return None

    """Convex assumption!"""
    # Compute a normal vector, that defines the plane direction
    normal_vector = np.asarray(direction_vector)
    normal_vector = normal_vector / norm(normal_vector)
    # Create a function to interpolate within the mask array
    interpolator = RegularGridInterpolator([np.arange(region.shape[i]) for i in range(3)], region)

    def is_inside(distance):
        coords = [start_point_np[i] + normal_vector[i] * distance for i in [0, 1, 2]]
        if any(i < 0 for i in coords):
            return 0
        if any(coords[i] > region.shape[i] - 1 for i in range(len(coords))):
            return 0
        # Evaluate the mask value at the interpolated coordinates
        mask_value = interpolator(coords)
        return mask_value > 0.5

    if not is_inside(0):
        return start_point_np
    count = 0
    min_v = 0
    if max_v is None:
        max_v = sum(region.shape)
    delta = max_v * 2
    while acc_delta < delta:
        bisection = (max_v - min_v) / 2 + min_v
        if is_inside(bisection):
            min_v = bisection
        else:
            max_v = bisection
        delta = max_v - min_v
        count += 1
    return start_point_np + normal_vector * ((min_v + max_v) / 2)

max_distance_ray_cast_convex ¶

max_distance_ray_cast_convex(region: NII, start_coord: COORDINATE | ndarray, direction_vector: ndarray, acc_delta: float = 5e-05) -> np.ndarray | None

Find the exit point of a ray inside a convex NII region.

Wraps the underlying numpy logic by extracting the voxel array from an NII object and delegating to bisection-based ray casting.

Parameters:

Name	Type	Description	Default
`region`	`NII`	`NII` object whose nonzero voxels define the convex region.	required
`start_coord`	`COORDINATE \| ndarray`	Starting coordinate `(x, y, z)` of the ray in voxel space. If already outside the region, the start coordinate is returned unchanged.	required
`direction_vector`	`ndarray`	Direction of the ray; will be normalised internally.	required
`acc_delta`	`float`	Bisection convergence threshold in voxel units. Defaults to 0.00005.	`5e-05`

Returns:

Type	Description
`ndarray \| None`	3-element numpy array with the approximate exit coordinate, or the start
`ndarray \| None`	coordinate if it lies outside the region. Returns `None` when
`ndarray \| None`	`start_coord` is `None`.

Source code in TPTBox/core/poi_fun/ray_casting.py

def max_distance_ray_cast_convex(
    region: NII,
    start_coord: COORDINATE | np.ndarray,
    direction_vector: np.ndarray,
    acc_delta: float = 0.00005,
) -> np.ndarray | None:
    """Find the exit point of a ray inside a convex NII region.

    Wraps the underlying numpy logic by extracting the voxel array from an
    ``NII`` object and delegating to bisection-based ray casting.

    Args:
        region: ``NII`` object whose nonzero voxels define the convex region.
        start_coord: Starting coordinate ``(x, y, z)`` of the ray in voxel
            space.  If already outside the region, the start coordinate is
            returned unchanged.
        direction_vector: Direction of the ray; will be normalised internally.
        acc_delta: Bisection convergence threshold in voxel units.
            Defaults to 0.00005.

    Returns:
        3-element numpy array with the approximate exit coordinate, or the start
        coordinate if it lies outside the region.  Returns ``None`` when
        ``start_coord`` is ``None``.
    """
    start_point_np = np.asarray(start_coord)
    if start_point_np is None:
        return None

    """Convex assumption!"""
    # Compute a normal vector, that defines the plane direction
    normal_vector = np.asarray(direction_vector)
    normal_vector = normal_vector / norm(normal_vector)
    # Create a function to interpolate within the mask array
    interpolator = RegularGridInterpolator([np.arange(region.shape[i]) for i in range(3)], region.get_array())

    def is_inside(distance):
        coords = [start_point_np[i] + normal_vector[i] * distance for i in [0, 1, 2]]
        if any(i < 0 for i in coords):
            return 0
        if any(coords[i] > region.shape[i] - 1 for i in range(len(coords))):
            return 0
        # Evaluate the mask value at the interpolated coordinates
        mask_value = interpolator(coords)
        return mask_value > 0.5

    if not is_inside(0):
        return start_point_np
    count = 0
    min_v = 0
    max_v = sum(region.shape)
    delta = max_v * 2
    while acc_delta < delta:
        bisection = (max_v - min_v) / 2 + min_v
        if is_inside(bisection):
            min_v = bisection
        else:
            max_v = bisection
        delta = max_v - min_v
        count += 1
    return start_point_np + normal_vector * ((min_v + max_v) / 2)

ray_cast_pixel_lvl ¶

ray_cast_pixel_lvl(start_point_np: ndarray, normal_vector: ndarray, shape: ndarray | tuple[int, ...], two_sided: bool = False) -> tuple[np.ndarray | None, np.ndarray | None]

Cast a ray at pixel/voxel resolution and return all integer voxel coordinates along it.

Traces the ray starting at start_point_np in the direction of normal_vector until it exits the volume defined by shape. When two_sided is True, the ray is also cast in the opposite direction and the two halves are concatenated.

Parameters:

Name	Type	Description	Default
`start_point_np`	`ndarray`	Starting voxel coordinate `(x, y, z)` of the ray.	required
`normal_vector`	`ndarray`	Direction vector of the ray; will be normalised to a unit vector internally.	required
`shape`	`ndarray \| tuple[int, ...]`	Volume shape `(dim0, dim1, dim2)` used to clip the ray.	required
`two_sided`	`bool`	When `True`, cast the ray in both directions and concatenate the results. Defaults to `False`.	`False`

Returns:

Type	Description
`ndarray \| None`	A tuple `(plane_coords, arange)` where
`ndarray \| None`	`plane_coords` is an integer array of shape `(N, 3)` with the voxel indices visited by the ray, or `None` on failure.
`tuple[ndarray \| None, ndarray \| None]`	`arange` is a float array of shape `(N,)` with the distance values along the ray for each voxel, or `None` on failure.

Source code in TPTBox/core/poi_fun/ray_casting.py

def ray_cast_pixel_lvl(
    start_point_np: np.ndarray,
    normal_vector: np.ndarray,
    shape: np.ndarray | tuple[int, ...],
    two_sided: bool = False,
) -> tuple[np.ndarray | None, np.ndarray | None]:
    """Cast a ray at pixel/voxel resolution and return all integer voxel coordinates along it.

    Traces the ray starting at ``start_point_np`` in the direction of
    ``normal_vector`` until it exits the volume defined by ``shape``.  When
    ``two_sided`` is ``True``, the ray is also cast in the opposite direction and
    the two halves are concatenated.

    Args:
        start_point_np: Starting voxel coordinate ``(x, y, z)`` of the ray.
        normal_vector: Direction vector of the ray; will be normalised to a
            unit vector internally.
        shape: Volume shape ``(dim0, dim1, dim2)`` used to clip the ray.
        two_sided: When ``True``, cast the ray in both directions and concatenate
            the results.  Defaults to ``False``.

    Returns:
        A tuple ``(plane_coords, arange)`` where

        * ``plane_coords`` is an integer array of shape ``(N, 3)`` with the
          voxel indices visited by the ray, or ``None`` on failure.
        * ``arange`` is a float array of shape ``(N,)`` with the distance
          values along the ray for each voxel, or ``None`` on failure.
    """
    normal_vector = unit_vector(normal_vector)

    def _calc_pixels(normal_vector, start_point_np):
        # Make a plane through start_point with the norm of "normal_vector", which is shifted by "shift" along the norm
        start_point_np = start_point_np.copy()
        num_pixel = np.abs(np.floor(np.max((np.array(shape) - start_point_np) / normal_vector))).item()
        arange = np.arange(0, min(num_pixel, 1000), step=1, dtype=float)
        coords = [start_point_np[i] + normal_vector[i] * arange for i in [0, 1, 2]]

        # Clip coordinates to region bounds
        for i in [0, 1, 2]:
            cut_off = (shape[i] <= np.floor(coords[i])).sum()
            if cut_off == 0:
                cut_off = (np.floor(coords[i]) <= 0).sum()
            if cut_off != 0:
                coords = [c[:-cut_off] for c in coords]
                arange = arange[:-cut_off]
        # Convert coordinates to integers for indexing
        int_coords = [c.astype(int) for c in coords]

        return np.stack(int_coords, -1), arange

    plane_coords, arange = _calc_pixels(normal_vector, start_point_np)
    if two_sided:
        plane_coords2, arange2 = _calc_pixels(-normal_vector, start_point_np)
        arange2 = -arange2
        plane_coords = np.concatenate([plane_coords, plane_coords2])
        arange = np.concatenate([arange, arange2]) - np.min(arange2)

    return plane_coords, arange

add_ray_to_img ¶

add_ray_to_img(start_point: ndarray | COORDINATE, normal_vector: ndarray, seg: NII, add_to_img: bool = True, inplace: bool = False, value: int = 0, dilate: int = 1) -> NII | None

Paint a ray into a segmentation NII image.

Casts a ray from start_point along normal_vector at voxel resolution, fills the visited voxels with distance values (or a fixed value), optionally dilates the ray, and optionally composites it onto the source segmentation.

Parameters:

Name	Type	Description	Default
`start_point`	`ndarray \| COORDINATE`	Starting voxel coordinate of the ray.	required
`normal_vector`	`ndarray`	Direction of the ray.	required
`seg`	`NII`	Segmentation `NII` that defines the volume shape and is used as the base when `add_to_img` is `True`.	required
`add_to_img`	`bool`	If `True`, the ray is overlaid on `seg` and the composite image is returned. If `False`, only the ray image is returned. Defaults to `True`.	`True`
`inplace`	`bool`	Modify `seg` in place when `add_to_img` is `True`. Defaults to `False`.	`False`
`value`	`int`	Fixed voxel value written along the ray. When 0, distance values along the ray are used instead. Defaults to 0.	`0`
`dilate`	`int`	Morphological dilation radius applied to the ray image. Set to 0 to skip dilation. Defaults to 1.	`1`

Returns:

Type	Description
`NII \| None`	The modified `NII` image, or `None` if the ray cast produced no
`NII \| None`	valid coordinates.

Source code in TPTBox/core/poi_fun/ray_casting.py

def add_ray_to_img(
    start_point: np.ndarray | COORDINATE,
    normal_vector: np.ndarray,
    seg: NII,
    add_to_img: bool = True,
    inplace: bool = False,
    value: int = 0,
    dilate: int = 1,
) -> NII | None:
    """Paint a ray into a segmentation NII image.

    Casts a ray from ``start_point`` along ``normal_vector`` at voxel
    resolution, fills the visited voxels with distance values (or a fixed
    ``value``), optionally dilates the ray, and optionally composites it onto
    the source segmentation.

    Args:
        start_point: Starting voxel coordinate of the ray.
        normal_vector: Direction of the ray.
        seg: Segmentation ``NII`` that defines the volume shape and is used as
            the base when ``add_to_img`` is ``True``.
        add_to_img: If ``True``, the ray is overlaid on ``seg`` and the
            composite image is returned.  If ``False``, only the ray image is
            returned.  Defaults to ``True``.
        inplace: Modify ``seg`` in place when ``add_to_img`` is ``True``.
            Defaults to ``False``.
        value: Fixed voxel value written along the ray.  When 0, distance
            values along the ray are used instead.  Defaults to 0.
        dilate: Morphological dilation radius applied to the ray image.
            Set to 0 to skip dilation.  Defaults to 1.

    Returns:
        The modified ``NII`` image, or ``None`` if the ray cast produced no
        valid coordinates.
    """
    start_point = np.array(start_point)
    plane_coords, arange = ray_cast_pixel_lvl(start_point, normal_vector, shape=seg.shape)
    if plane_coords is None:
        return None
    selected_arr = np.zeros(seg.shape, dtype=seg.dtype)
    selected_arr[plane_coords[..., 0], plane_coords[..., 1], plane_coords[..., 2]] = arange if value == 0 else value
    ray = seg.set_array(selected_arr)
    if dilate != 0:
        ray.dilate_msk_(dilate)
    if add_to_img:
        if not inplace:
            seg = seg.copy()
        seg[ray != 0] = ray[ray != 0]
        return seg
    return ray

add_spline_to_img ¶

add_spline_to_img(seg: NII, poi: POI, location: int = 50, add_to_img: bool = True, override_seg: bool = True, value: int = 100, dilate: int = 2) -> NII

Fit and paint a POI spline onto a segmentation NII image.

Computes a cubic spline through the POI centroids at location, converts it to voxel coordinates, paints each spline point with value, dilates the result, and (optionally) composites it onto seg.

Parameters:

Name	Type	Description	Default
`seg`	`NII`	Base segmentation `NII` image used for shape and affine.	required
`poi`	`POI`	`POI` object from which the spline is computed.	required
`location`	`int`	Subregion ID used as the spline anchor points. Defaults to 50 (`Vertebra_Corpus`).	`50`
`add_to_img`	`bool`	If `True`, the spline is overlaid on `seg`. Defaults to `True`.	`True`
`override_seg`	`bool`	When overlaying, whether to overwrite existing nonzero voxels in `seg`. Defaults to `True`.	`True`
`value`	`int`	Voxel label written along the spline. Defaults to 100.	`100`
`dilate`	`int`	Morphological dilation radius. Defaults to 2.	`2`

Returns:

Type	Description
`NII`	The composite or standalone spline `NII` image.

Source code in TPTBox/core/poi_fun/ray_casting.py

def add_spline_to_img(
    seg: NII,
    poi: POI,
    location: int = 50,
    add_to_img: bool = True,
    override_seg: bool = True,
    value: int = 100,
    dilate: int = 2,
) -> NII:
    """Fit and paint a POI spline onto a segmentation NII image.

    Computes a cubic spline through the POI centroids at ``location``,
    converts it to voxel coordinates, paints each spline point with ``value``,
    dilates the result, and (optionally) composites it onto ``seg``.

    Args:
        seg: Base segmentation ``NII`` image used for shape and affine.
        poi: ``POI`` object from which the spline is computed.
        location: Subregion ID used as the spline anchor points.
            Defaults to 50 (``Vertebra_Corpus``).
        add_to_img: If ``True``, the spline is overlaid on ``seg``.
            Defaults to ``True``.
        override_seg: When overlaying, whether to overwrite existing nonzero
            voxels in ``seg``.  Defaults to ``True``.
        value: Voxel label written along the spline.  Defaults to 100.
        dilate: Morphological dilation radius.  Defaults to 2.

    Returns:
        The composite or standalone spline ``NII`` image.
    """
    cor, _ = poi.fit_spline(location=location, vertebra=True)
    spline = seg.copy() * 0
    # spline.rescale_()
    for x, y, z in cor:
        spline[round(x), round(y), round(z)] = value
    spline.dilate_msk_(dilate)
    # spline.resample_from_to_(a)
    if add_to_img:
        if override_seg:
            seg[spline != 0] = spline[spline != 0]
        else:
            cond = np.logical_and(spline != 0, seg == 0)
            seg[cond] = spline[cond]
        return seg
    return spline

shift_point ¶

shift_point(poi: POI, vert_id: int, bb: tuple[slice, slice, slice] | None, start_point: Location = Location.Vertebra_Corpus, direction: DIRECTIONS | None = 'R', log: Logger_Interface = _log) -> np.ndarray | None

Shift a POI start point laterally by a fraction of the vertebra width.

Estimates the vertebra width from articular-process or costal-process POIs and returns a new coordinate displaced from start_point along direction by a scaled fraction of that width. Sacrum vertebrae without arcus are skipped.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object with pre-computed anatomical landmarks.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice] \| None`	Bounding-box slices used to convert global POI coordinates to local coordinates within the cropped sub-volume.	required
`start_point`	`Location`	Location enum member or already-resolved numpy coordinate from which the shift originates. Defaults to `Location.Vertebra_Corpus`.	`Vertebra_Corpus`
`direction`	`DIRECTIONS \| None`	Anatomical direction letter (`"R"`, `"L"`, etc.) for the shift, or `None` to return the raw local start point without shifting. Defaults to `"R"`.	`'R'`
`log`	`Logger_Interface`	Logger used for warning and error messages.	`_log`

Returns:

Type	Description
`ndarray \| None`	Shifted voxel coordinate as a 3-element numpy array, or `None` when
`ndarray \| None`	required POIs are missing or the vertebra is excluded.

Source code in TPTBox/core/poi_fun/ray_casting.py

def shift_point(
    poi: POI,
    vert_id: int,
    bb: tuple[slice, slice, slice] | None,
    start_point: Location = Location.Vertebra_Corpus,
    direction: DIRECTIONS | None = "R",
    log: Logger_Interface = _log,
) -> np.ndarray | None:
    """Shift a POI start point laterally by a fraction of the vertebra width.

    Estimates the vertebra width from articular-process or costal-process POIs
    and returns a new coordinate displaced from ``start_point`` along
    ``direction`` by a scaled fraction of that width.  Sacrum vertebrae without
    arcus are skipped.

    Args:
        poi: ``POI`` object with pre-computed anatomical landmarks.
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices used to convert global POI coordinates to local
            coordinates within the cropped sub-volume.
        start_point: Location enum member or already-resolved numpy coordinate
            from which the shift originates.  Defaults to
            ``Location.Vertebra_Corpus``.
        direction: Anatomical direction letter (``"R"``, ``"L"``, etc.) for the
            shift, or ``None`` to return the raw local start point without
            shifting.  Defaults to ``"R"``.
        log: Logger used for warning and error messages.

    Returns:
        Shifted voxel coordinate as a 3-element numpy array, or ``None`` when
        required POIs are missing or the vertebra is excluded.
    """
    if vert_id in sacrum_w_o_arcus:
        return

    if direction is None:
        return to_local_np(start_point, bb, poi, vert_id, log)
    sup_articular_right = sup_articular_left = None
    if Vertebra_Instance.is_sacrum(vert_id):
        if (vert_id, Location.Costal_Process_Left) not in poi:
            return
        sup_articular_right = to_local_np(Location.Costal_Process_Right, bb, poi, vert_id, log)
        sup_articular_left = to_local_np(Location.Costal_Process_Left, bb, poi, vert_id, log)
        factor = 6.0
    if sup_articular_left is None or sup_articular_right is None:
        sup_articular_right = to_local_np(Location.Superior_Articular_Right, bb, poi, vert_id, log, verbose=False)
        sup_articular_left = to_local_np(Location.Superior_Articular_Left, bb, poi, vert_id, log, verbose=False)
        factor = 3.0
    if sup_articular_left is None or sup_articular_right is None:
        sup_articular_right = to_local_np(Location.Inferior_Articular_Right, bb, poi, vert_id, log)
        sup_articular_left = to_local_np(Location.Inferior_Articular_Left, bb, poi, vert_id, log)
        factor = 2.0
    if sup_articular_left is None or sup_articular_right is None:
        return
    if vert_id <= 11:
        factor *= (12 - vert_id) / 11 + 1
    vertebra_width = np.linalg.norm(sup_articular_right - sup_articular_left)
    shift = vertebra_width / factor
    normal_vector = get_direction(direction, poi, vert_id)  # / np.array(poi.zoom)
    normal_vector = normal_vector / norm(normal_vector)
    start_point_np = to_local_np(start_point, bb, poi, vert_id, log) if isinstance(start_point, Location) else start_point
    if start_point_np is None:
        return None
    return start_point_np + normal_vector * shift

max_distance_ray_cast_convex_poi ¶

max_distance_ray_cast_convex_poi(poi: POI, region: NII, vert_id: int, bb: tuple[slice, slice, slice] | None, normal_vector_points: tuple[Location, Location] | DIRECTIONS = 'R', start_point: Location | ndarray = Location.Vertebra_Corpus, log: Logger_Interface = _log, acc_delta: float = 5e-05) -> np.ndarray | None

Find the boundary point of a convex region along a direction derived from POI landmarks.

Resolves the ray start coordinate and direction vector from the POI object, then delegates to :func:max_distance_ray_cast_convex to locate the exit point.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object with pre-computed anatomical landmarks and direction vectors.	required
`region`	`NII`	`NII` image whose nonzero voxels define the convex region to cast the ray into.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice] \| None`	Bounding-box slices that map between global POI coordinates and the local sub-volume coordinate system.	required
`normal_vector_points`	`tuple[Location, Location] \| DIRECTIONS`	Either a `DIRECTIONS` string (e.g. `"R"`, `"P"`) that is resolved from the vertebra orientation, or a 2-tuple of `Location` members whose vector difference defines the ray direction. Defaults to `"R"`.	`'R'`
`start_point`	`Location \| ndarray`	Starting point of the ray: a `Location` enum member (resolved via the POI) or a pre-computed numpy coordinate. Defaults to `Location.Vertebra_Corpus`.	`Vertebra_Corpus`
`log`	`Logger_Interface`	Logger used for warning and error messages.	`_log`
`acc_delta`	`float`	Bisection convergence threshold in voxel units. Defaults to 0.00005.	`5e-05`

Returns:

Type	Description
`ndarray \| None`	3-element numpy array with the approximate exit coordinate in the local
`ndarray \| None`	sub-volume coordinate system, or `None` when the start point or
`ndarray \| None`	direction cannot be resolved.

Source code in TPTBox/core/poi_fun/ray_casting.py

def max_distance_ray_cast_convex_poi(
    poi: POI,
    region: NII,
    vert_id: int,
    bb: tuple[slice, slice, slice] | None,
    normal_vector_points: tuple[Location, Location] | DIRECTIONS = "R",
    start_point: Location | np.ndarray = Location.Vertebra_Corpus,
    log: Logger_Interface = _log,
    acc_delta: float = 0.00005,
) -> np.ndarray | None:
    """Find the boundary point of a convex region along a direction derived from POI landmarks.

    Resolves the ray start coordinate and direction vector from the ``POI``
    object, then delegates to :func:`max_distance_ray_cast_convex` to locate
    the exit point.

    Args:
        poi: ``POI`` object with pre-computed anatomical landmarks and direction
            vectors.
        region: ``NII`` image whose nonzero voxels define the convex region to
            cast the ray into.
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices that map between global POI coordinates and the
            local sub-volume coordinate system.
        normal_vector_points: Either a ``DIRECTIONS`` string (e.g. ``"R"``,
            ``"P"``) that is resolved from the vertebra orientation, or a
            2-tuple of ``Location`` members whose vector difference defines the
            ray direction.  Defaults to ``"R"``.
        start_point: Starting point of the ray: a ``Location`` enum member
            (resolved via the POI) or a pre-computed numpy coordinate.
            Defaults to ``Location.Vertebra_Corpus``.
        log: Logger used for warning and error messages.
        acc_delta: Bisection convergence threshold in voxel units.
            Defaults to 0.00005.

    Returns:
        3-element numpy array with the approximate exit coordinate in the local
        sub-volume coordinate system, or ``None`` when the start point or
        direction cannot be resolved.
    """
    start_point_np = to_local_np(start_point, bb, poi, vert_id, log) if isinstance(start_point, Location) else start_point
    if start_point_np is None:
        return None

    """Convex assumption!"""
    # Compute a normal vector, that defines the plane direction
    if isinstance(normal_vector_points, str):
        try:
            normal_vector = get_direction(normal_vector_points, poi, vert_id)
        except KeyError:
            if vert_id not in sacrum_w_o_arcus:
                log.on_fail(f"region={vert_id},DIRECTIONS={normal_vector_points} is missing")
            return
            # raise KeyError(f"region={label},subregion={loc.value} is missing.")
        normal_vector /= norm(normal_vector)

    else:
        try:
            b = to_local_np(normal_vector_points[1], bb, poi, vert_id, log)
            if b is None:
                return None
            a = to_local_np(normal_vector_points[0], bb, poi, vert_id, log)
            normal_vector = b - a
            normal_vector = normal_vector / norm(normal_vector)
            log.on_fail(f"ray_cast used with old normal_vector_points {normal_vector_points}")
        except TypeError as e:
            log.on_fail("TypeError", e)
            return None
    return max_distance_ray_cast_convex(region, start_point_np, normal_vector, acc_delta)

calculate_pca_normal_np ¶

calculate_pca_normal_np(segmentation: ndarray, pca_component: int, zoom: tuple[float, ...] | ndarray | None = None, verbose: bool = False) -> np.ndarray

Computes the normal vector of a segmented region using Principal Component Analysis (PCA).

Parameters:¶

segmentation : np.ndarray A binary mask where nonzero values indicate the segmented region. pca_component : int, optional The principal component index to return as the normal vector. - 0: The primary axis (direction of greatest variance, often the main elongation). - 1: The secondary axis (orthogonal to the primary, capturing the second-most variance). - 2: The third axis (typically the normal direction to the structure in 3D). zoom : tuple or array-like, optional If provided, scales the normal vector by the inverse of the voxel size to account for anisotropic resolution. verbose : bool, optional If True, prints the principal component vectors for debugging.

Returns:¶

normal_vector : np.ndarray The selected principal component as a normal vector.

Usage:¶

Use pca_component=2 when you want the normal to a surface-like structure. If analyzing an elongated structure (e.g., a vessel or bone), pca_component=0 gives the primary axis, while pca_component=1 provides the secondary direction.

Source code in TPTBox/core/poi_fun/ray_casting.py

def calculate_pca_normal_np(
    segmentation: np.ndarray,
    pca_component: int,
    zoom: tuple[float, ...] | np.ndarray | None = None,
    verbose: bool = False,
) -> np.ndarray:
    """Computes the normal vector of a segmented region using Principal Component Analysis (PCA).

    Parameters:
    ----------
    segmentation : np.ndarray
        A binary mask where nonzero values indicate the segmented region.
    pca_component : int, optional
        The principal component index to return as the normal vector.
        - `0`: The primary axis (direction of greatest variance, often the main elongation).
        - `1`: The secondary axis (orthogonal to the primary, capturing the second-most variance).
        - `2`: The third axis (typically the normal direction to the structure in 3D).
    zoom : tuple or array-like, optional
        If provided, scales the normal vector by the inverse of the voxel size to account for anisotropic resolution.
    verbose : bool, optional
        If True, prints the principal component vectors for debugging.

    Returns:
    -------
    normal_vector : np.ndarray
        The selected principal component as a normal vector.

    Usage:
    ------
    Use `pca_component=2` when you want the normal to a surface-like structure.
    If analyzing an elongated structure (e.g., a vessel or bone), `pca_component=0` gives the primary axis,
    while `pca_component=1` provides the secondary direction.
    """
    # Get indices of segmented region (assuming segmentation is a binary mask)
    points = np.argwhere(segmentation > 0)
    # Perform PCA to find the principal axes
    pca = PCA(n_components=3)
    pca.fit(points)
    # First, second, and third principal components
    if verbose:
        print(f"Main Axis (PC1): {pca.components_[0]}")
        print(f"Secondary Axis (PC2): {pca.components_[1]}")
        print(f"Third Axis (PC3): {pca.components_[2]}")
    normal_vector = pca.components_[pca_component]
    if zoom is not None:
        normal_vector = normal_vector / np.array(zoom)
    return normal_vector

set_label_above_3_point_plane ¶

set_label_above_3_point_plane(array: NII | ndarray, p1: ndarray | list[float], p2: ndarray | list[float], p3: ndarray | list[float], value: float = 0, invert: Literal[-1, 1] = 1, mask: ndarray | NII | bool = True, inplace: bool = False) -> NII | np.ndarray

Set all values in a 3D array above a plane defined by three non-collinear points to a specified value.

Parameters:¶

array : NII | np.ndarray A 3D NumPy array or an NII object representing the volume data. p1, p2, p3 : array-like Three (x, y, z) points defining the plane. value : int or float, optional (default=0) The value to set for all elements above the plane. inf : Literal[-1, 1], optional (default=1) Controls the direction of "above": - 1 means values superior to the plane (default). - -1 means values inferior to the plane. If the input is an NII object with an inferior-superior orientation, this will be adjusted accordingly.

Returns:¶

np.ndarray The modified 3D array with values set above the plane.

Notes:¶

The plane is defined by the equation ax + by + cz + d = 0, where (a, b, c) is the normal vector.
Uses meshgrid to construct a 3D grid and determine which values lie above the plane.

Example:¶

import numpy as np

data = np.random.rand(300, 300, 300)
p1 = [100, 100, 50]
p2 = [100, 50, 100]
p3 = [50, 100, 100]

result = set_above_3_point_plane(data, p1, p2, p3, value=0)

Source code in TPTBox/core/poi_fun/ray_casting.py

def set_label_above_3_point_plane(
    array: NII | np.ndarray,
    p1: np.ndarray | list[float],
    p2: np.ndarray | list[float],
    p3: np.ndarray | list[float],
    value: float = 0,
    invert: Literal[-1, 1] = 1,
    mask: np.ndarray | NII | bool = True,
    inplace: bool = False,
) -> NII | np.ndarray:
    """Set all values in a 3D array above a plane defined by three non-collinear points to a specified value.

    Parameters:
    -----------
    array : NII | np.ndarray
        A 3D NumPy array or an NII object representing the volume data.
    p1, p2, p3 : array-like
        Three (x, y, z) points defining the plane.
    value : int or float, optional (default=0)
        The value to set for all elements above the plane.
    inf : Literal[-1, 1], optional (default=1)
        Controls the direction of "above":
        - `1` means values superior to the plane (default).
        - `-1` means values inferior to the plane.
        If the input is an NII object with an inferior-superior orientation, this will be adjusted accordingly.

    Returns:
    --------
    np.ndarray
        The modified 3D array with values set above the plane.

    Notes:
    ------
    - The plane is defined by the equation `ax + by + cz + d = 0`, where `(a, b, c)` is the normal vector.
    - Uses `meshgrid` to construct a 3D grid and determine which values lie above the plane.

    Example:
    --------
    ```python
    import numpy as np

    data = np.random.rand(300, 300, 300)
    p1 = [100, 100, 50]
    p2 = [100, 50, 100]
    p3 = [50, 100, 100]

    result = set_above_3_point_plane(data, p1, p2, p3, value=0)
    ```
    """
    import numpy as np

    if not inplace:
        array = array.copy()
    if isinstance(array, NII) and array.orientation[array.get_axis("S")] == "I":
        # for NII inf = 1 means Superior
        invert *= -1

    # Define three 3D points that form a plane
    p1 = np.array(p1)
    p2 = np.array(p2)
    p3 = np.array(p3)

    # Compute the normal vector of the plane
    normal = np.cross(p2 - p1, p3 - p1)
    a, b, c = normal
    d = -np.dot(normal, p1)

    # Define a 3D grid of points
    shape = array.shape
    x, y, z = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]), np.arange(shape[2]), indexing="ij")

    # Compute the plane equation for each (x, y)
    plane_z = (-a * x - b * y - d) / c

    # Create the 3D array and set values above the plane to 0
    array[np.logical_and(mask, z * invert > plane_z)] = value
    return array

Vertebra Direction¶

TPTBox.core.poi_fun.vertebra_direction ¶

calc_orientation_of_vertebra_PIR ¶

calc_orientation_of_vertebra_PIR(poi: POI | None, vert: NII, subreg: NII, spline_subreg_point_id=Location.Vertebra_Corpus, source_subreg_point_id=Location.Vertebra_Corpus, subreg_id=Location.Spinal_Canal, do_fill_back: bool = False, spine_plot_path: None | str = None, save_normals_in_info=False, _orientation_version=0) -> tuple[POI, NII | None]

Calculate the orientation of vertebrae using PIR (Posterior, Inferior, Right) DIRECTIONS.

Parameters:

Name	Type	Description	Default
`poi`	`POI \| None`	Point of interest. If None, computed from `vert` and `subreg`.	required
`vert`	`NII`	Vertebra (full).	required
`subreg`	`NII`	Subregion (full).	required
`spline_subreg_point_id`	`Location`	Subregion point ID for spline computation. Defaults to Location.Vertebra_Corpus.	`Vertebra_Corpus`
`source_subreg_point_id`	`Location`	Source subregion point ID. Defaults to Location.Vertebra_Corpus.	`Vertebra_Corpus`
`subreg_id`	`Location`	Subregion ID. Defaults to Location.Spinal_Canal.	`Spinal_Canal`
`do_fill_back`	`bool`	Whether to fill back. Defaults to False.	`False`
`spine_plot_path`	`None \| str`	Path to spine plot. Defaults to None.	`None`
`save_normals_in_info`	`bool`	Whether to save normals in info. Defaults to False.	`False`

Returns:

Type	Description
`tuple[POI, NII \| None]`	Tuple[POI, NII \| None]: Point of interest and filled back NII.

Source code in TPTBox/core/poi_fun/vertebra_direction.py

def calc_orientation_of_vertebra_PIR(
    poi: POI | None,
    vert: NII,
    subreg: NII,
    spline_subreg_point_id=Location.Vertebra_Corpus,
    source_subreg_point_id=Location.Vertebra_Corpus,
    subreg_id=Location.Spinal_Canal,
    do_fill_back: bool = False,
    spine_plot_path: None | str = None,
    save_normals_in_info=False,
    _orientation_version=0,
) -> tuple[POI, NII | None]:
    """Calculate the orientation of vertebrae using PIR (Posterior, Inferior, Right) DIRECTIONS.

    Args:
        poi (POI | None): Point of interest. If None, computed from `vert` and `subreg`.
        vert (NII): Vertebra (full).
        subreg (NII): Subregion (full).
        spline_subreg_point_id (Location, optional): Subregion point ID for spline computation. Defaults to Location.Vertebra_Corpus.
        source_subreg_point_id (Location, optional): Source subregion point ID. Defaults to Location.Vertebra_Corpus.
        subreg_id (Location, optional): Subregion ID. Defaults to Location.Spinal_Canal.
        do_fill_back (bool, optional): Whether to fill back. Defaults to False.
        spine_plot_path (None | str, optional): Path to spine plot. Defaults to None.
        save_normals_in_info (bool, optional): Whether to save normals in info. Defaults to False.

    Returns:
        Tuple[POI, NII | None]: Point of interest and filled back NII.
    """
    assert poi is None or poi.zoom is not None
    from TPTBox import calc_centroids

    if _orientation_version != 0:
        warn("out dated _orientation_version; Set is to 0", stacklevel=1)
    # Step 1 compute the up direction
    # check if label 50 is already computed in POI
    if poi is None or spline_subreg_point_id.value not in poi.keys_subregion():
        poi = calc_poi_from_subreg_vert(vert, subreg, extend_to=poi, subreg_id=spline_subreg_point_id)
    # compute Spline in ISO space
    poi_iso = poi.rescale().reorient()
    body_spline, body_spline_der = poi_iso.fit_spline(location=spline_subreg_point_id, vertebra=True)
    # Step 2 compute the back direction by Spinosus_Process or arcus
    intersection_target = [Location.Spinosus_Process, Location.Arcus_Vertebrae]
    # We compute everything in iso space
    subreg_iso = subreg.rescale().reorient()
    if _orientation_version == 0:
        target_labels = subreg_iso.extract_label(intersection_target).get_array()
        # we want to see more of the Spinosus_Process and Arcus_Vertebrae than we cut with the plane. Should reduce randomness.
        # The ideal solution would be to make a projection onto the plane. Instead we fill values that have a vertical distanc of 10 mm up and down. This approximates the projection on to the plane.
        # Without this we have the chance to miss most of the arcus and spinosus, witch leads to instability in the direction.
        # TODO this will fail if the vertebra is not roughly aligned with S/I-direction
        for _ in range(15):
            target_labels[:, :-1] += target_labels[:, 1:]
            target_labels[:, 1:] += target_labels[:, :-1]
        target_labels = np.clip(target_labels, 0, 1)
    elif _orientation_version == 1:
        target_labels = subreg_iso.extract_label(intersection_target).get_array()
    elif _orientation_version == 2:
        target_labels = subreg_iso.extract_label(list(range(40, 49))).get_array()
    else:
        raise NotImplementedError(_orientation_version)
    out = target_labels * 0
    fill_back_nii = subreg_iso.copy() if do_fill_back else None
    fill_back = out.copy() if do_fill_back else None
    down_vector: dict[int, np.ndarray] = {}
    # Draw a plain with the up_vector an cut it with intersection_target
    for reg_label, _, cords in poi_iso.extract_subregion(source_subreg_point_id).items():
        # calculate_normal_vector
        distances = np.sqrt(np.sum((body_spline - np.array(cords)) ** 2, -1))
        normal_vector_post = body_spline_der[np.argmin(distances)]
        normal_vector_post /= np.linalg.norm(normal_vector_post)
        down_vector[reg_label] = normal_vector_post.copy()
        # create_plane_coords
        # The main axis will be treated differently
        idx = [_plane_dict[i] for i in subreg_iso.orientation]
        axis = idx.index(_plane_dict["S"])
        # assert axis == np.argmax(np.abs(normal_vector)).item()
        dims = [0, 1, 2]
        dims.remove(axis)
        dim1, dim2 = dims
        # Make a plane through start_point with the norm of "normal_vector", which is shifted by "shift" along the norm
        start_point_np = np.array(cords)
        start_point_np[axis] = start_point_np[axis]
        shift_total = -start_point_np.dot(normal_vector_post)
        xx, yy = np.meshgrid(range(subreg_iso.shape[dim1]), range(subreg_iso.shape[dim2]))  # type: ignore
        zz = (-normal_vector_post[dim1] * xx - normal_vector_post[dim2] * yy - shift_total) * 1.0 / normal_vector_post[axis]
        z_max = subreg_iso.shape[axis] - 1
        zz[zz < 0] = 0
        zz[zz > z_max] = 0
        plane_coords = np.zeros([xx.shape[0], xx.shape[1], 3])
        plane_coords[:, :, axis] = zz
        plane_coords[:, :, dim1] = xx
        plane_coords[:, :, dim2] = yy
        plane_coords = plane_coords.astype(int)
        # create_subregion
        # 1 where the selected subreg is, else 0
        select = subreg_iso.get_array() * 0
        select[plane_coords[:, :, 0], plane_coords[:, :, 1], plane_coords[:, :, 2]] = 1
        out[out == 0] += (target_labels * select * reg_label)[out == 0]

        if fill_back is not None:
            fill_back[np.logical_and(select == 1, fill_back == 0)] = reg_label
    if fill_back is not None and fill_back_nii is not None:
        subreg_sar = subreg_iso.set_array(fill_back).reorient(("S", "A", "R"))
        fill_back = subreg_sar.get_array()
        x_slice = np.ones_like(fill_back[0]) * np.max(fill_back) + 1
        for i in range(fill_back.shape[0]):
            curr_slice = fill_back[i]
            cond = np.where(curr_slice != 0)
            x_slice[cond] = np.minimum(curr_slice[cond], x_slice[cond])
            fill_back[i] = x_slice
        subreg_sar.set_array(fill_back).reorient(poi.orientation).rescale_(poi.zoom)
        arr = subreg_sar.get_array()
        fill_back_nii.set_array_(arr)

    ret = calc_centroids(subreg_iso.set_array(out), second_stage=subreg_id, extend_to=poi_iso.copy(), inplace=True)

    poi._vert_orientation_pir = {}
    if save_normals_in_info:
        poi.info["vert_orientation_PIR"] = poi._vert_orientation_pir
    # calc posterior vector and the crossproduct
    for vert_id, normal_down in down_vector.items():
        try:
            # get two points and compute the direction:
            a = np.array(ret[vert_id : subreg_id.value]) - 1
            b = np.array(ret[vert_id : source_subreg_point_id.value]) - 1
            normal_vector_post = a - b
            normal_vector_post = normal_vector_post / norm(normal_vector_post)
            poi._vert_orientation_pir[vert_id] = (
                normal_vector_post,
                normal_down,
                np.cross(normal_vector_post, normal_down),
            )

            ### MAKE DIRECTIONS POIs ###
            ret[vert_id, Location.Vertebra_Direction_Posterior] = tuple(ret[vert_id, source_subreg_point_id] + normal_vector_post * 10)
            ret[vert_id, Location.Vertebra_Direction_Inferior] = tuple(ret[vert_id, source_subreg_point_id] + normal_down * 10)
            ret[vert_id, Location.Vertebra_Direction_Right] = tuple(
                ret[vert_id:source_subreg_point_id] + np.cross(normal_vector_post, normal_down * 10)
            )
        except KeyError as e:
            if vert_id not in sacrum_w_o_direction:
                _log.on_fail(f"calc_orientation_of_vertebra_PIR {vert_id=} - KeyError=", e)

    # if make_thicker:
    ret.remove(*ret.extract_subregion(subreg_id).keys())

    if spine_plot_path is not None:
        make_spine_plot(ret, body_spline, vert, spine_plot_path)

    ret = ret.resample_from_to(poi)  # type: ignore

    return poi.join_right_(ret), fill_back_nii

get_direction ¶

get_direction(d: DIRECTIONS, poi: POI, vert_id: int, verbose: bool = True) -> np.ndarray

Return the unit vector for an anatomical direction of a vertebra.

Retrieves the pre-computed PIR orientation vectors for vertebra vert_id from poi and returns the one corresponding to direction d. C1 has no independent direction and falls back to C2.

Parameters:

Name	Type	Description	Default
`d`	`DIRECTIONS`	Anatomical direction letter (`"P"`, `"A"`, `"I"`, `"S"`, `"R"`, or `"L"`).	required
`poi`	`POI`	`POI` object with pre-computed `Vertebra_Direction_*` landmarks.	required
`vert_id`	`int`	Vertebra identifier. If 1 (C1), direction is taken from C2.	required
`verbose`	`bool`	Emit a warning when C1 falls back to C2. Defaults to `True`.	`True`

Returns:

Type	Description
`ndarray`	Unit vector (3-element numpy array) pointing in direction `d` in the
`ndarray`	image coordinate frame.

Raises:

Type	Description
`KeyError`	If orientation data for the requested vertebra is absent.

Source code in TPTBox/core/poi_fun/vertebra_direction.py

def get_direction(d: DIRECTIONS, poi: POI, vert_id: int, verbose: bool = True) -> np.ndarray:
    """Return the unit vector for an anatomical direction of a vertebra.

    Retrieves the pre-computed PIR orientation vectors for vertebra ``vert_id``
    from ``poi`` and returns the one corresponding to direction ``d``.  C1 has
    no independent direction and falls back to C2.

    Args:
        d: Anatomical direction letter (``"P"``, ``"A"``, ``"I"``, ``"S"``,
            ``"R"``, or ``"L"``).
        poi: ``POI`` object with pre-computed ``Vertebra_Direction_*`` landmarks.
        vert_id: Vertebra identifier.  If 1 (C1), direction is taken from C2.
        verbose: Emit a warning when C1 falls back to C2.  Defaults to ``True``.

    Returns:
        Unit vector (3-element numpy array) pointing in direction ``d`` in the
        image coordinate frame.

    Raises:
        KeyError: If orientation data for the requested vertebra is absent.
    """
    if vert_id == 1:
        _log.on_warning("C1 has no direction computation use C2 as direction", verbose=verbose)
        vert_id = 2
    P, I, R = get_vert_direction_PIR(poi, vert_id, to_pir=False)  # noqa: N806
    if d == "P":
        return P
    elif d == "A":
        return -P
    elif d == "I":
        return I
    elif d == "S":
        return -I
    elif d == "R":
        return R
    elif d == "L":
        return -R
    else:
        never_called(d)

get_vert_direction_PIR ¶

get_vert_direction_PIR(poi: POI, vert_id: int, do_norm: bool = True, to_pir: bool = True) -> Vertebra_Orientation

Retrieve the vertebra orientation vectors from the POI.

Must be called after :func:calc_orientation_of_vertebra_PIR has populated the Vertebra_Direction_* POI entries.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object with `Vertebra_Corpus`, `Vertebra_Direction_Posterior`, `Vertebra_Direction_Inferior`, and `Vertebra_Direction_Right` entries for the given vertebra.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`do_norm`	`bool`	Normalise each direction vector to unit length. Defaults to `True`.	`True`
`to_pir`	`bool`	Rescale and reorient the POI to isotropic PIR space before computing the directions. When `True` the result is also cached in `poi._vert_orientation_pir`. Defaults to `True`.	`True`

Returns:

Type	Description
`Vertebra_Orientation`	Tuple of three unit vectors `(P, I, R)` representing the Posterior,
`Vertebra_Orientation`	Inferior, and Right directions of the vertebra in the (optionally
`Vertebra_Orientation`	reoriented) coordinate frame.

Source code in TPTBox/core/poi_fun/vertebra_direction.py

def get_vert_direction_PIR(poi: POI, vert_id: int, do_norm: bool = True, to_pir: bool = True) -> Vertebra_Orientation:
    """Retrieve the vertebra orientation vectors from the POI.

    Must be called after :func:`calc_orientation_of_vertebra_PIR` has populated
    the ``Vertebra_Direction_*`` POI entries.

    Args:
        poi: ``POI`` object with ``Vertebra_Corpus``,
            ``Vertebra_Direction_Posterior``, ``Vertebra_Direction_Inferior``,
            and ``Vertebra_Direction_Right`` entries for the given vertebra.
        vert_id: Vertebra identifier (integer label).
        do_norm: Normalise each direction vector to unit length.
            Defaults to ``True``.
        to_pir: Rescale and reorient the POI to isotropic PIR space before
            computing the directions.  When ``True`` the result is also cached
            in ``poi._vert_orientation_pir``.  Defaults to ``True``.

    Returns:
        Tuple of three unit vectors ``(P, I, R)`` representing the Posterior,
        Inferior, and Right directions of the vertebra in the (optionally
        reoriented) coordinate frame.
    """
    if vert_id in poi._vert_orientation_pir and to_pir:
        return poi._vert_orientation_pir[vert_id]  # Elusive buffer of iso/PIR directions.
    poi = poi.extract_subregion(
        Location.Vertebra_Corpus,
        Location.Vertebra_Direction_Posterior,
        Location.Vertebra_Direction_Inferior,
        Location.Vertebra_Direction_Right,
    )
    if to_pir:
        poi = poi.rescale(verbose=False).reorient(verbose=False)

    def n(x):
        return x / norm(x) if do_norm else x

    center = np.array(poi[vert_id : Location.Vertebra_Corpus])
    post = np.array(poi[vert_id : Location.Vertebra_Direction_Posterior])
    down = np.array(poi[vert_id : Location.Vertebra_Direction_Inferior])
    right = np.array(poi[vert_id : Location.Vertebra_Direction_Right])
    out = n(post - center), n(down - center), n(right - center)
    if to_pir:
        poi._vert_orientation_pir[vert_id] = out

    return out

get_vert_direction_matrix ¶

get_vert_direction_matrix(poi: POI, vert_id: int, to_pir: bool = False) -> tuple[np.ndarray, np.ndarray]

Return the change-of-basis matrices between the image frame and the vertebra PIR frame.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object with pre-computed vertebra direction landmarks.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`to_pir`	`bool`	Whether to convert the POI to isotropic PIR space before computing. Defaults to `False`.	`False`

Returns:

Type	Description
`ndarray`	A tuple `(to_vert_orient, from_vert_orient)` where
`ndarray`	`to_vert_orient` transforms image-frame vectors into the vertebra PIR frame (shape `(3, 3)`).
`tuple[ndarray, ndarray]`	`from_vert_orient` transforms vertebra-frame vectors back into the image frame (shape `(3, 3)`).

Source code in TPTBox/core/poi_fun/vertebra_direction.py

def get_vert_direction_matrix(poi: POI, vert_id: int, to_pir: bool = False) -> tuple[np.ndarray, np.ndarray]:
    """Return the change-of-basis matrices between the image frame and the vertebra PIR frame.

    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``.

    Returns:
        A tuple ``(to_vert_orient, from_vert_orient)`` where

        * ``to_vert_orient`` transforms image-frame vectors into the vertebra
          PIR frame (shape ``(3, 3)``).
        * ``from_vert_orient`` transforms vertebra-frame vectors back into the
          image frame (shape ``(3, 3)``).
    """
    P, I, R = get_vert_direction_PIR(poi, vert_id=vert_id, to_pir=to_pir)  # noqa: N806
    from_vert_orient = np.stack([P, I, R], axis=1)
    to_vert_orient = np.linalg.inv(from_vert_orient)
    return to_vert_orient, from_vert_orient

calc_center_spinal_cord ¶

calc_center_spinal_cord(poi: POI, subreg: NII, spline_subreg_point_id: Location | list[Location] = Location.Vertebra_Corpus, source_subreg_point_id: Location = Location.Vertebra_Corpus, subreg_id=Location.Spinal_Canal, intersection_target: list[Location] | None = None, _fill_inplace: NII | None = None, add_dense=False) -> POI

Calculate the center of the spinal cord within a specified region.

Parameters: - poi (POI): Point of Interest object containing relevant data. - subreg (NII): Neuroimaging subregion data for analysis. - spline_subreg_point_id (int, optional): The height of the region to analyze (default is Location.Vertebra_Corpus). - subreg_id (int, optional): The identifier for the subregion (default is Location.Spinal_Canal). - intersection_target (List[Location], optional): A list of locations for intersection analysis (default is [Location.Spinal_Cord, Location.Spinal_Canal]).

Returns: - POI: Point of Interest object with calculated centroids.

This function calculates the center of the spinal cord within a specified region using a spline fit to the given Point of Interest (POI) data and a NII. It first extracts a subregion of interest based on the provided spline_subreg_point_id and then calculates the center using geometric operations.

Note: - The poi object should contain the relevant data for spline fitting. - The subreg object should be a NII subregion containing spine with number 60 and 61. - The subreg_id parameter is an optional identifier for the subregion. - The intersection_target parameter allows you to specify the locations to intersect (e.g., spinal cord and spinal canal).

Example usage:

poi = POI(...)
subreg = NII(...)
updated_poi = calc_center_spinal_cord(
    poi,
    subreg,
    spline_subreg_point_id=Location.Vertebra_Corpus,
    subreg_id=Location.Spinal_Canal,
    intersection_target=[Location.Spinal_Cord, Location.Spinal_Canal],
)

Source code in TPTBox/core/poi_fun/vertebra_direction.py

def calc_center_spinal_cord(
    poi: POI,
    subreg: NII,
    spline_subreg_point_id: Location | list[Location] = Location.Vertebra_Corpus,
    source_subreg_point_id: Location = Location.Vertebra_Corpus,
    subreg_id=Location.Spinal_Canal,
    intersection_target: list[Location] | None = None,
    _fill_inplace: NII | None = None,
    add_dense=False,
) -> POI:
    """Calculate the center of the spinal cord within a specified region.

    Parameters:
    - poi (POI): Point of Interest object containing relevant data.
    - subreg (NII): Neuroimaging subregion data for analysis.
    - spline_subreg_point_id (int, optional): The height of the region to analyze (default is Location.Vertebra_Corpus).
    - subreg_id (int, optional): The identifier for the subregion (default is Location.Spinal_Canal).
    - intersection_target (List[Location], optional): A list of locations for intersection analysis (default is [Location.Spinal_Cord, Location.Spinal_Canal]).

    Returns:
    - POI: Point of Interest object with calculated centroids.

    This function calculates the center of the spinal cord within a specified region using a spline fit to the
    given Point of Interest (POI) data and a NII. It first extracts a subregion of
    interest based on the provided spline_subreg_point_id and then calculates the center using geometric operations.

    Note:
    - The `poi` object should contain the relevant data for spline fitting.
    - The `subreg` object should be a NII subregion containing spine with number 60 and 61.
    - The `subreg_id` parameter is an optional identifier for the subregion.
    - The `intersection_target` parameter allows you to specify the locations to intersect (e.g., spinal cord and spinal canal).

    Example usage:
    ```python
    poi = POI(...)
    subreg = NII(...)
    updated_poi = calc_center_spinal_cord(
        poi,
        subreg,
        spline_subreg_point_id=Location.Vertebra_Corpus,
        subreg_id=Location.Spinal_Canal,
        intersection_target=[Location.Spinal_Cord, Location.Spinal_Canal],
    )
    ```
    """
    from TPTBox import calc_centroids

    if intersection_target is None:
        intersection_target = [Location.Spinal_Cord, Location.Spinal_Canal]
    assert _fill_inplace is None or subreg == _fill_inplace
    poi_iso = poi.rescale()
    if add_dense and (2, Location.Dens_axis) in poi_iso:
        sx = spline_subreg_point_id[0] if isinstance(spline_subreg_point_id, Sequence) else spline_subreg_point_id
        poi_iso[1, sx] = poi_iso[2, Location.Dens_axis]
        poi_iso[1, sx] = poi_iso[2, Location.Dens_axis]
    subreg_iso = subreg.rescale()
    body_spline, body_spline_der = poi_iso.fit_spline(location=spline_subreg_point_id, vertebra=True)
    target_labels = subreg_iso.extract_label(intersection_target).get_array()
    out = target_labels * 0
    fill_back = out.copy() if _fill_inplace is not None else None
    for reg_label, _, cords in poi_iso.extract_subregion(source_subreg_point_id).items():
        # calculate_normal_vector
        distances = np.sqrt(np.sum((body_spline - np.array(cords)) ** 2, -1))
        normal_vector = body_spline_der[np.argmin(distances)]
        normal_vector /= np.linalg.norm(normal_vector)
        # create_plane_coords
        # The main axis will be treated differently
        idx = [_plane_dict[i] for i in subreg_iso.orientation]
        axis = idx.index(_plane_dict["S"])
        # assert axis == np.argmax(np.abs(normal_vector)).item()
        dims = [0, 1, 2]
        dims.remove(axis)
        dim1, dim2 = dims
        # Make a plane through start_point with the norm of "normal_vector", which is shifted by "shift" along the norm
        start_point_np = np.array(cords)
        start_point_np[axis] = start_point_np[axis]
        shift_total = -start_point_np.dot(normal_vector)
        xx, yy = np.meshgrid(range(subreg_iso.shape[dim1]), range(subreg_iso.shape[dim2]))  # type: ignore
        zz = (-normal_vector[dim1] * xx - normal_vector[dim2] * yy - shift_total) * 1.0 / normal_vector[axis]
        z_max = subreg_iso.shape[axis] - 1
        zz[zz < 0] = 0
        zz[zz > z_max] = 0
        plane_coords = np.zeros([xx.shape[0], xx.shape[1], 3])
        plane_coords[:, :, axis] = zz
        plane_coords[:, :, dim1] = xx
        plane_coords[:, :, dim2] = yy
        plane_coords = plane_coords.astype(int)
        # create_subregion
        # 1 where the selected subreg is, else 0
        select = subreg_iso.get_array() * 0
        select[plane_coords[:, :, 0], plane_coords[:, :, 1], plane_coords[:, :, 2]] = 1
        out += target_labels * select * reg_label

        if fill_back is not None:
            fill_back[np.logical_and(select == 1, fill_back == 0)] = reg_label
    if fill_back is not None:
        assert _fill_inplace is not None
        subreg_iso2 = subreg_iso.set_array(fill_back).reorient(("S", "A", "R"))
        fill_back = subreg_iso2.get_array()
        x_slice = np.ones_like(fill_back[0]) * np.max(fill_back) + 1
        for i in range(fill_back.shape[0]):
            curr_slice = fill_back[i]
            cond = np.where(curr_slice != 0)
            x_slice[cond] = np.minimum(curr_slice[cond], x_slice[cond])
            fill_back[i] = x_slice

        arr = subreg_iso2.set_array(fill_back).reorient(poi.orientation).rescale_(poi.zoom).get_array()
        _fill_inplace.set_array_(arr)
    ret = calc_centroids(subreg_iso.set_array(out), second_stage=subreg_id, extend_to=poi_iso.extract_subregion(subreg_id), inplace=True)
    ret.rescale_(poi.zoom)
    return poi.join_left_(ret)

Vertebra Non-Centroid POIs¶

TPTBox.core.poi_fun.vertebra_pois_non_centroids ¶

Strategy_Pattern ¶

Implements the Strategy design pattern by encapsulating different strategies as callable objects.

Parameters:

Name	Type	Description	Default
`target`	`Location`	The target location for which this strategy is defined.	required
`strategy`	`Callable`	The strategy function that implements the desired behavior.	required
`prerequisite`	`set[Location] \| None`	A set of prerequisite locations that must be satisfied before applying this strategy. Defaults to None.	`None`
`**args`		Additional keyword arguments to be passed to the strategy function.	`{}`

Attributes:

Name	Type	Description
`target`	`Location`	The target location for which this strategy is defined.
`args`	`dict`	Additional keyword arguments to be passed to the strategy function.
`prerequisite`	`set[Location]`	A set of prerequisite locations that must be satisfied before applying this strategy.
`strategy`	`Callable`	The strategy function that implements the desired behavior.

Note

The strategy function should accept the following arguments: - poi (POI): The point of interest. - current_subreg (NII): The current subregion. - vert_id (int): The vertex ID. - bb: The bounding box. - log (Logger_Interface, optional): The logger interface. Defaults to _log, which should be defined globally.

Example

def strategy_function(poi, current_subreg, location, log, vert_id, bb, **kwargs): ... # Strategy implementation ... pass strategy = Strategy_Pattern(target_location, strategy_function, prerequisite={prerequisite_location}, additional_arg=value) result = strategy(poi, current_subreg, vert_id, bb)

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

class Strategy_Pattern:
    """Implements the Strategy design pattern by encapsulating different strategies as callable objects.

    Args:
        target (Location): The target location for which this strategy is defined.
        strategy (Callable): The strategy function that implements the desired behavior.
        prerequisite (set[Location] | None, optional): A set of prerequisite locations that must be satisfied before applying this strategy. Defaults to None.
        **args: Additional keyword arguments to be passed to the strategy function.

    Attributes:
        target (Location): The target location for which this strategy is defined.
        args (dict): Additional keyword arguments to be passed to the strategy function.
        prerequisite (set[Location]): A set of prerequisite locations that must be satisfied before applying this strategy.
        strategy (Callable): The strategy function that implements the desired behavior.

    Note:
        The strategy function should accept the following arguments:
        - poi (POI): The point of interest.
        - current_subreg (NII): The current subregion.
        - vert_id (int): The vertex ID.
        - bb: The bounding box.
        - log (Logger_Interface, optional): The logger interface. Defaults to _log, which should be defined globally.

    Example:
        >>> def strategy_function(poi, current_subreg, location, log, vert_id, bb, **kwargs):
        ...     # Strategy implementation
        ...     pass
        >>> strategy = Strategy_Pattern(target_location, strategy_function, prerequisite={prerequisite_location}, additional_arg=value)
        >>> result = strategy(poi, current_subreg, vert_id, bb)
    """

    def __init__(
        self,
        target: Location,
        strategy: Callable,
        prerequisite: set[Location] | None = None,
        prio=0,
        sakrum=False,
        **args,
    ) -> None:
        self.target = target
        self.args = args
        self.sacrum = sakrum
        if prerequisite is None:
            prerequisite = set()
        if "direction" in args.keys():
            prerequisite.add(Location.Vertebra_Direction_Inferior)
        for i in args.values():
            if isinstance(i, Location):
                prerequisite.add(i)
            elif isinstance(i, Sequence):
                for j in i:
                    if isinstance(j, Location):
                        prerequisite.add(j)
        self.prerequisite = prerequisite
        self.strategy = strategy
        all_poi_functions[target.value] = self
        self._prio = prio

    def __call__(
        self,
        poi: POI,
        current_subreg: NII,
        vert_id: int,
        bb,
        log: Logger_Interface = _log,
    ):
        """Execute the POI computation strategy for the given vertebra, returning the result or None on failure."""
        try:
            return self.strategy(
                poi=poi,
                current_subreg=current_subreg,
                location=self.target,
                log=log,
                vert_id=vert_id,
                bb=bb,
                **self.args,
            )
        except Exception:
            _log.print_error()
            return None

    def prority(self) -> int:
        """Return the scheduling priority for this strategy (lower value runs first)."""
        return self.target.value + self._prio

prority ¶

prority() -> int

Return the scheduling priority for this strategy (lower value runs first).

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

def prority(self) -> int:
    """Return the scheduling priority for this strategy (lower value runs first)."""
    return self.target.value + self._prio

Strategy_Pattern_Side_Effect ¶

Bases: Strategy_Pattern

Strategy marker for POIs computed as a side effect of another strategy.

Registers target as a POI that does not need its own compute call because it is produced implicitly when prerequisite is computed. The callable is a no-op; the dependency is tracked in :data:pois_computed_by_side_effect.

Parameters:

Name	Type	Description	Default
`target`	`Location`	The `Location` that is produced as a side effect.	required
`prerequisite`	`Location`	The `Location` whose computation produces `target`.	required
`**args`		Additional keyword arguments forwarded to :class:`Strategy_Pattern`.	`{}`

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

class Strategy_Pattern_Side_Effect(Strategy_Pattern):
    """Strategy marker for POIs computed as a side effect of another strategy.

    Registers ``target`` as a POI that does not need its own compute call
    because it is produced implicitly when ``prerequisite`` is computed.  The
    callable is a no-op; the dependency is tracked in
    :data:`pois_computed_by_side_effect`.

    Args:
        target: The ``Location`` that is produced as a side effect.
        prerequisite: The ``Location`` whose computation produces ``target``.
        **args: Additional keyword arguments forwarded to
            :class:`Strategy_Pattern`.
    """

    def __init__(self, target: Location, prerequisite: Location, **args) -> None:
        super().__init__(target, _strategy_side_effect, {prerequisite}, **args)
        pois_computed_by_side_effect[target.value] = prerequisite

Strategy_Computed_Before ¶

Bases: Strategy_Pattern

Strategy marker for POIs that must be computed before the pipeline runs.

Use this when a Location is computed by an earlier, separate step (e.g. spinal canal centroids) rather than by a generic strategy function. The callable is a no-op; the class simply registers the prerequisite dependencies.

Parameters:

Name	Type	Description	Default
`target`	`Location`	The `Location` that is already computed before the pipeline.	required
`*prerequisite`	`Location`	One or more `Location` members that must be present before `target` can be used.	`()`
`**args`		Additional keyword arguments forwarded to :class:`Strategy_Pattern`.	`{}`

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

class Strategy_Computed_Before(Strategy_Pattern):
    """Strategy marker for POIs that must be computed before the pipeline runs.

    Use this when a ``Location`` is computed by an earlier, separate step
    (e.g. spinal canal centroids) rather than by a generic strategy function.
    The callable is a no-op; the class simply registers the prerequisite
    dependencies.

    Args:
        target: The ``Location`` that is already computed before the pipeline.
        *prerequisite: One or more ``Location`` members that must be present
            before ``target`` can be used.
        **args: Additional keyword arguments forwarded to
            :class:`Strategy_Pattern`.
    """

    def __init__(self, target: Location, *prerequisite: Location, **args) -> None:
        super().__init__(target, _strategy_side_effect, set(prerequisite), **args)

run_poi_pipeline ¶

run_poi_pipeline(vert: NII, subreg: NII, poi_path: Path, logger: Logger_Interface = _log) -> None

Compute all POIs from vertebra and subregion segmentations and save them to disk.

Wraps :func:~TPTBox.calc_poi_from_subreg_vert for every Location member and writes the result as a JSON file.

Parameters:

Name	Type	Description	Default
`vert`	`NII`	Vertebra instance segmentation `NII` (one label per vertebra).	required
`subreg`	`NII`	Subregion segmentation `NII` (spine substructure labels).	required
`poi_path`	`Path`	Output path for the JSON POI file. Used as a buffer file so that an existing file is read and extended rather than recomputed from scratch.	required
`logger`	`Logger_Interface`	Logger instance for progress and error messages.	`_log`

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

def run_poi_pipeline(vert: NII, subreg: NII, poi_path: Path, logger: Logger_Interface = _log) -> None:
    """Compute all POIs from vertebra and subregion segmentations and save them to disk.

    Wraps :func:`~TPTBox.calc_poi_from_subreg_vert` for every ``Location``
    member and writes the result as a JSON file.

    Args:
        vert: Vertebra instance segmentation ``NII`` (one label per vertebra).
        subreg: Subregion segmentation ``NII`` (spine substructure labels).
        poi_path: Output path for the JSON POI file.  Used as a buffer file
            so that an existing file is read and extended rather than
            recomputed from scratch.
        logger: Logger instance for progress and error messages.
    """
    poi = calc_poi_from_subreg_vert(vert, subreg, buffer_file=poi_path, save_buffer_file=True, subreg_id=list(Location), verbose=logger)
    poi.save(poi_path)

add_prerequisites ¶

add_prerequisites(locs: Sequence[Location]) -> list[Location]

Expand a list of Location members with all transitive prerequisites.

Walks the :data:all_poi_functions registry and iteratively adds any Location entries that are required by the provided list (and their requirements, and so on) until no new entries are found.

Parameters:

Name	Type	Description	Default
`locs`	`Sequence[Location]`	Initial sequence of `Location` members for which POIs are to be computed.	required

Returns:

Type	Description
`list[Location]`	Sorted list of `Location` members including all original entries and
`list[Location]`	every transitive prerequisite, ordered by `Location.value`.

Raises:

Type	Description
`UserWarning`	If a dependency cycle is detected (loop limit of 1000 iterations reached).

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

def add_prerequisites(locs: Sequence[Location]) -> list[Location]:
    """Expand a list of ``Location`` members with all transitive prerequisites.

    Walks the :data:`all_poi_functions` registry and iteratively adds any
    ``Location`` entries that are required by the provided list (and their
    requirements, and so on) until no new entries are found.

    Args:
        locs: Initial sequence of ``Location`` members for which POIs are to
            be computed.

    Returns:
        Sorted list of ``Location`` members including all original entries and
        every transitive prerequisite, ordered by ``Location.value``.

    Raises:
        UserWarning: If a dependency cycle is detected (loop limit of 1000
            iterations reached).
    """
    addendum = set()
    locs2 = set(locs)
    loop_var = locs2
    i = 0
    while i != 1000:  # Prevent Deadlock
        for l in loop_var:
            if l.value in all_poi_functions:
                for prereq in all_poi_functions[l.value].prerequisite:
                    if prereq not in locs:
                        addendum.add(prereq)
        if len(addendum) == 0:
            break
        locs2 = addendum | locs2
        loop_var = addendum
        addendum = set()
        i += 1
    else:
        warnings.warn("Deadlock in add_prerequisites", stacklevel=10)
    return sorted(list(locs2), key=lambda x: x.value)  # type: ignore # noqa: C414

compute_non_centroid_pois ¶

compute_non_centroid_pois(poi: POI, locations: Sequence[Location] | Location, vert: NII, subreg: NII, _vert_ids: Sequence[int] | None = None, log: Logger_Interface = _log, verbose: bool | None = True, _orientation_version: int = 0) -> None

Compute all requested non-centroid anatomical landmarks and store them in poi in place.

Runs the full non-centroid POI pipeline:

Vertebra orientation (PIR direction vectors) — always computed first if Location.Vertebra_Direction_Inferior is requested.
Global landmarks: spinal canal / cord centres, intervertebral disc POIs, dense-axis tip.
Per-vertebra landmarks via the registered :class:Strategy_Pattern functions (extreme points, ray casts, corner finders, etc.).
Intervertebral disc (IVD) POIs.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will be modified in place. Must contain centroid information for `vert` and `subreg`.	required
`locations`	`Sequence[Location] \| Location`	One or more `Location` members identifying the POIs to compute. A single `Location` is accepted and converted to a list internally.	required
`vert`	`NII`	Vertebra instance segmentation `NII` (one label per vertebra).	required
`subreg`	`NII`	Spine subregion segmentation `NII`.	required
`_vert_ids`	`Sequence[int] \| None`	Explicit list of vertebra integer labels to process. When `None`, all labels present in `vert` are used.	`None`
`log`	`Logger_Interface`	Logger for progress and error messages.	`_log`
`verbose`	`bool \| None`	Verbosity passed to logging calls. `None` silences everything, `True` enables all messages.	`True`
`_orientation_version`	`int`	Internal version selector for the vertebra orientation algorithm; keep at 0 (default) for production use.	`0`

Source code in TPTBox/core/poi_fun/vertebra_pois_non_centroids.py

def compute_non_centroid_pois(  # noqa: C901
    poi: POI,
    locations: Sequence[Location] | Location,
    vert: NII,
    subreg: NII,
    _vert_ids: Sequence[int] | None = None,
    log: Logger_Interface = _log,
    verbose: bool | None = True,
    _orientation_version: int = 0,
) -> None:
    """Compute all requested non-centroid anatomical landmarks and store them in ``poi`` in place.

    Runs the full non-centroid POI pipeline:

    1. Vertebra orientation (PIR direction vectors) — always computed first if
       ``Location.Vertebra_Direction_Inferior`` is requested.
    2. Global landmarks: spinal canal / cord centres, intervertebral disc
       POIs, dense-axis tip.
    3. Per-vertebra landmarks via the registered :class:`Strategy_Pattern`
       functions (extreme points, ray casts, corner finders, etc.).
    4. Intervertebral disc (IVD) POIs.

    Args:
        poi: ``POI`` object that will be modified in place.  Must contain
            centroid information for ``vert`` and ``subreg``.
        locations: One or more ``Location`` members identifying the POIs to
            compute.  A single ``Location`` is accepted and converted to a
            list internally.
        vert: Vertebra instance segmentation ``NII`` (one label per vertebra).
        subreg: Spine subregion segmentation ``NII``.
        _vert_ids: Explicit list of vertebra integer labels to process.
            When ``None``, all labels present in ``vert`` are used.
        log: Logger for progress and error messages.
        verbose: Verbosity passed to logging calls.  ``None`` silences
            everything, ``True`` enables all messages.
        _orientation_version: Internal version selector for the vertebra
            orientation algorithm; keep at 0 (default) for production use.
    """
    if _vert_ids is None:
        _vert_ids = vert.unique()

    locations = list(locations) if isinstance(locations, Sequence) else [locations]
    ### STEP 1 Vert Direction###
    assert 52 not in poi.keys_region()

    if Location.Vertebra_Direction_Inferior in locations:
        log.on_text("Compute Vertebra DIRECTIONS", verbose=verbose)
        ### Calc vertebra direction; We always need them, so we just compute them. ###
        sub_regions = poi.keys_subregion()
        if any(a.value not in sub_regions for a in vert_directions):
            poi, _ = calc_orientation_of_vertebra_PIR(
                poi, vert, subreg, do_fill_back=False, save_normals_in_info=False, _orientation_version=_orientation_version
            )
            [locations.remove(i) for i in vert_directions if i in locations]

    locations = [pois_computed_by_side_effect.get(l.value, l) for l in locations]
    locations = sorted(
        set(locations),
        key=lambda x: all_poi_functions[x.value].prority() if x.value in all_poi_functions else x.value,
    )  # type: ignore
    log.on_text("Calc pois from subregion id", {l.name for l in locations}, verbose=verbose)
    ### DENSE ###
    if Location.Dens_axis in locations and 2 in _vert_ids and (2, Location.Dens_axis.value) not in poi:
        a = subreg * vert.extract_label(2)
        bb = a.compute_crop()
        a = a.apply_crop(bb)
        s = [Location.Vertebra_Corpus, Location.Vertebra_Corpus_border]
        if a.sum() != 0:
            strategy_extreme_points(poi, a, location=Location.Dens_axis, direction=["S", "P"], vert_id=2, subreg_id=s, bb=bb)
    ### STEP 2 (Other global non centroid poi; Spinal heights ###

    if Location.Spinal_Canal in locations:
        locations.remove(Location.Spinal_Canal)
        subregs_ids = subreg.unique()
        _a = Location.Spinal_Canal.value in subregs_ids or Location.Spinal_Canal.value in subregs_ids
        if _a and Location.Spinal_Canal.value not in poi.keys_subregion():
            poi = calc_center_spinal_cord(poi, subreg, add_dense=True)

    if Location.Spinal_Cord in locations:
        locations.remove(Location.Spinal_Cord)
        subregs_ids = subreg.unique()
        v = Location.Spinal_Cord.value
        if (v in subregs_ids or Location.Spinal_Cord.value in subregs_ids) and v not in poi.keys_subregion():
            poi = calc_center_spinal_cord(
                poi,
                subreg,
                source_subreg_point_id=Location.Vertebra_Disc,
                subreg_id=Location.Spinal_Cord,
                add_dense=True,
                intersection_target=[Location.Spinal_Cord],
            )
    if Location.Spinal_Canal_ivd_lvl in locations:
        locations.remove(Location.Spinal_Canal_ivd_lvl)
        subregs_ids = subreg.unique()
        v = Location.Spinal_Canal_ivd_lvl.value
        if (v in subregs_ids or Location.Spinal_Cord.value in subregs_ids) and v not in poi.keys_subregion():
            poi = calc_center_spinal_cord(
                poi, subreg, source_subreg_point_id=Location.Vertebra_Disc, subreg_id=Location.Spinal_Canal_ivd_lvl, add_dense=True
            )
    # Step 3 Compute on individual Vertebras
    ivd_location = set()

    for vert_id in _vert_ids:
        if vert_id >= 39:
            continue
        current_vert = vert.extract_label(vert_id)
        bb = current_vert.compute_crop()
        current_vert.apply_crop_(bb)
        current_subreg = subreg.apply_crop(bb) * current_vert
        for location in locations:
            # TODO IVD STUFF
            if location.value <= 50:
                continue
            if (vert_id, location.value) in poi:
                continue
            if location in [
                Location.Implant_Entry_Left,
                Location.Implant_Entry_Right,
                Location.Implant_Target_Left,
                Location.Implant_Target_Right,
                Location.Spinal_Canal,
            ]:
                continue
            if location in ivd_pois_list:
                ivd_location.add(location)
                continue
            if location.value in all_poi_functions:
                fun = all_poi_functions[location.value]

                if not fun.sacrum and vert_id in _sacrum:
                    continue

                fun(poi, current_subreg, vert_id, bb=bb, log=log)
            else:
                warnings.warn(f"NotImplementedError: {location}", stacklevel=0)
    # Step 4 IVD points
    if len(ivd_location) != 0:
        poi = calculate_IVD_POI(vert, subreg, poi, ivd_location)

Pixel-Based Point Finder¶

TPTBox.core.poi_fun.pixel_based_point_finder ¶

get_nearest_neighbor ¶

get_nearest_neighbor(p: ndarray, sr_msk: ndarray, region_label: int) -> np.ndarray

Return the voxel in a label mask that is closest to a query point.

Parameters:

Name	Type	Description	Default
`p`	`ndarray`	Query point as a 1-D or column array of shape `(3,)` or `(3, 1)`.	required
`sr_msk`	`ndarray`	3-D integer array containing segmentation labels.	required
`region_label`	`int`	Label value whose voxels are searched for the nearest neighbour.	required

Returns:

Type	Description
`ndarray`	1-D integer array `(x, y, z)` of the voxel in `sr_msk` with label
`ndarray`	`region_label` that minimises the Euclidean distance to `p`.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def get_nearest_neighbor(
    p: np.ndarray,
    sr_msk: np.ndarray,
    region_label: int,
) -> np.ndarray:
    """Return the voxel in a label mask that is closest to a query point.

    Args:
        p: Query point as a 1-D or column array of shape ``(3,)`` or ``(3, 1)``.
        sr_msk: 3-D integer array containing segmentation labels.
        region_label: Label value whose voxels are searched for the nearest
            neighbour.

    Returns:
        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()]

max_distance_ray_cast_pixel_level ¶

max_distance_ray_cast_pixel_level(poi: POI, region: NII, vert_id: int, bb: tuple[slice, slice, slice], normal_vector_points: tuple[Location, Location] | DIRECTIONS = 'R', start_point: Location | ndarray = Location.Vertebra_Corpus, two_sided=False, log: Logger_Interface = _log) -> tuple[int, ...] | None

Calculate the maximum distance ray cast in a region.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	Point of interest.	required
`region`	`NII`	Region to cast rays in.	required
`vert_id`	`int`	Label of the region in `region`.	required
`bb`	`Tuple[slice, slice, slice]`	Bounding box coordinates.	required
`normal_vector_points`	`Union[Tuple[Location, Location], DIRECTIONS]`	Points defining the normal vector or the direction. Defaults to "R".	`'R'`
`start_point`	`Location`	Starting point of the ray. Defaults to Location.Vertebra_Corpus.	`Vertebra_Corpus`
`log`	`Logger_Interface`	Logger interface. Defaults to _log.	`_log`

Returns:

Type	Description
`tuple[int, ...] \| None`	Tuple[int, int, int]: The coordinates of the maximum distance ray cast.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def max_distance_ray_cast_pixel_level(
    poi: POI,
    region: NII,
    vert_id: int,
    bb: tuple[slice, slice, slice],
    normal_vector_points: tuple[Location, Location] | DIRECTIONS = "R",
    start_point: Location | np.ndarray = Location.Vertebra_Corpus,
    two_sided=False,
    log: Logger_Interface = _log,
) -> tuple[int, ...] | None:
    """Calculate the maximum distance ray cast in a region.

    Args:
        poi (POI): Point of interest.
        region (NII): Region to cast rays in.
        vert_id (int): Label of the region in `region`.
        bb (Tuple[slice, slice, slice]): Bounding box coordinates.
        normal_vector_points (Union[Tuple[Location, Location], DIRECTIONS], optional):
            Points defining the normal vector or the direction. Defaults to "R".
        start_point (Location, optional): Starting point of the ray. Defaults to Location.Vertebra_Corpus.
        log (Logger_Interface, optional): Logger interface. Defaults to _log.

    Returns:
        Tuple[int, int, int]: The coordinates of the maximum distance ray cast.
    """
    plane_coords, arange = ray_cast_pixel_level_from_poi(
        poi,
        region,
        vert_id,
        bb,
        normal_vector_points,
        start_point,
        log=log,
        two_sided=two_sided,
    )
    if plane_coords is None:
        return None
    selected_arr = np.zeros(region.shape)
    selected_arr[plane_coords[..., 0], plane_coords[..., 1], plane_coords[..., 2]] = arange
    selected_arr = selected_arr * region.get_array()
    out = tuple(np.unravel_index(np.argmax(selected_arr, axis=None), selected_arr.shape))
    return out

ray_cast_pixel_level_from_poi ¶

ray_cast_pixel_level_from_poi(poi: POI, region: NII, vert_id: int, bb: tuple[slice, slice, slice], normal_vector_points: tuple[Location, Location] | DIRECTIONS = 'R', start_point: Location | ndarray = Location.Vertebra_Corpus, log: Logger_Interface = _log, two_sided: bool = False) -> tuple[np.ndarray | None, np.ndarray | None]

Perform pixel-level ray casting in a region using POI-derived direction.

Resolves the ray start coordinate and direction vector from the POI object, then delegates to :func:~TPTBox.core.poi_fun.ray_casting.ray_cast_pixel_lvl to enumerate all voxels along the ray at integer resolution.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object providing landmark coordinates and orientation.	required
`region`	`NII`	`NII` image defining the volume shape for boundary clipping.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice]`	Bounding-box slices mapping global POI coordinates to local sub-volume indices.	required
`normal_vector_points`	`tuple[Location, Location] \| DIRECTIONS`	Ray direction — a `DIRECTIONS` string resolved from the vertebra orientation, or a 2-tuple of `Location` members whose difference defines the direction. Defaults to `"R"`.	`'R'`
`start_point`	`Location \| ndarray`	Ray origin — a `Location` member resolved from `poi` or a pre-computed numpy coordinate. Defaults to `Location.Vertebra_Corpus`.	`Vertebra_Corpus`
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`
`two_sided`	`bool`	Cast the ray in both directions. Defaults to `False`.	`False`

Returns:

Type	Description
`ndarray \| None`	`(plane_coords, arange)` where `plane_coords` is an integer array of
`ndarray \| None`	shape `(N, 3)` with visited voxel indices and `arange` is a float
`tuple[ndarray \| None, ndarray \| None]`	array of distance values, or `(None, None)` on failure.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def ray_cast_pixel_level_from_poi(
    poi: POI,
    region: NII,
    vert_id: int,
    bb: tuple[slice, slice, slice],
    normal_vector_points: tuple[Location, Location] | DIRECTIONS = "R",
    start_point: Location | np.ndarray = Location.Vertebra_Corpus,
    log: Logger_Interface = _log,
    two_sided: bool = False,
) -> tuple[np.ndarray | None, np.ndarray | None]:
    """Perform pixel-level ray casting in a region using POI-derived direction.

    Resolves the ray start coordinate and direction vector from the ``POI``
    object, then delegates to :func:`~TPTBox.core.poi_fun.ray_casting.ray_cast_pixel_lvl`
    to enumerate all voxels along the ray at integer resolution.

    Args:
        poi: ``POI`` object providing landmark coordinates and orientation.
        region: ``NII`` image defining the volume shape for boundary clipping.
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices mapping global POI coordinates to local
            sub-volume indices.
        normal_vector_points: Ray direction — a ``DIRECTIONS`` string resolved
            from the vertebra orientation, or a 2-tuple of ``Location`` members
            whose difference defines the direction.  Defaults to ``"R"``.
        start_point: Ray origin — a ``Location`` member resolved from ``poi``
            or a pre-computed numpy coordinate.  Defaults to
            ``Location.Vertebra_Corpus``.
        log: Logger for warning and error messages.
        two_sided: Cast the ray in both directions.  Defaults to ``False``.

    Returns:
        ``(plane_coords, arange)`` where ``plane_coords`` is an integer array of
        shape ``(N, 3)`` with visited voxel indices and ``arange`` is a float
        array of distance values, or ``(None, None)`` on failure.
    """
    from TPTBox.core.poi_fun.ray_casting import ray_cast_pixel_lvl

    start_point_np = to_local_np(start_point, bb, poi, vert_id, log) if isinstance(start_point, Location) else start_point
    if start_point_np is None:
        return None, None

    # Compute a normal vector, that defines the plane direction
    if isinstance(normal_vector_points, str):
        normal_vector = get_direction(normal_vector_points, poi, vert_id)  # / np.array(poi.zoom)
        normal_vector = normal_vector / norm(normal_vector)
    else:
        try:
            b = to_local_np(normal_vector_points[1], bb, poi, vert_id, log)
            if b is None:
                return None, None
            a = to_local_np(normal_vector_points[0], bb, poi, vert_id, log)
            normal_vector = b - a
            normal_vector = normal_vector / norm(normal_vector)
            log.on_fail(f"ray_cast used with old normal_vector_points {normal_vector_points}")
        except TypeError as e:
            log.on_fail("TypeError", e)
            return None, None
    return ray_cast_pixel_lvl(start_point_np, normal_vector, region.shape, two_sided=two_sided)

get_extreme_point_by_vert_direction ¶

get_extreme_point_by_vert_direction(poi: POI, region: NII, vert_id: int, direction: Sequence[DIRECTIONS] | DIRECTIONS | tuple[DIRECTIONS, float] | Sequence[tuple[DIRECTIONS, float]] = 'I') -> np.ndarray | None

Find the voxel in a binary region that is most extreme along a vertebra-relative direction.

Projects all nonzero voxels in region into the vertebra PIR frame and returns the voxel index with the maximum weighted score along the requested direction(s).

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object providing the vertebra orientation matrix for `vert_id`.	required
`region`	`NII`	Binary `NII` image (value 1 for foreground) aligned to the same space as `poi`.	required
`vert_id`	`int`	Vertebra identifier used to look up the orientation matrix.	required
`direction`	`Sequence[DIRECTIONS] \| DIRECTIONS \| tuple[DIRECTIONS, float] \| Sequence[tuple[DIRECTIONS, float]]`	Anatomical direction(s) to optimise. May be: A single direction letter (e.g. `"I"`). A sequence of direction letters (e.g. `["P", "I"]`). A `(direction, weight)` tuple or a sequence of such tuples for weighted combinations. Defaults to `"I"`.	`'I'`

Returns:

Type	Description
`ndarray \| None`	Integer voxel index array `(x, y, z)` of the most extreme point, or
`ndarray \| None`	`None` if the orientation matrix cannot be retrieved.

Note

region must contain binary values (1 for foreground, 0 for background). Zoom scaling from poi is applied to ensure physically meaningful distances.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def get_extreme_point_by_vert_direction(
    poi: POI,
    region: NII,
    vert_id: int,
    direction: Sequence[DIRECTIONS] | DIRECTIONS | tuple[DIRECTIONS, float] | Sequence[tuple[DIRECTIONS, float]] = "I",
) -> np.ndarray | None:
    """Find the voxel in a binary region that is most extreme along a vertebra-relative direction.

    Projects all nonzero voxels in ``region`` into the vertebra PIR frame and
    returns the voxel index with the maximum weighted score along the requested
    direction(s).

    Args:
        poi: ``POI`` object providing the vertebra orientation matrix for
            ``vert_id``.
        region: Binary ``NII`` image (value 1 for foreground) aligned to the
            same space as ``poi``.
        vert_id: Vertebra identifier used to look up the orientation matrix.
        direction: Anatomical direction(s) to optimise.  May be:

            * A single direction letter (e.g. ``"I"``).
            * A sequence of direction letters (e.g. ``["P", "I"]``).
            * A ``(direction, weight)`` tuple or a sequence of such tuples for
              weighted combinations.

            Defaults to ``"I"``.

    Returns:
        Integer voxel index array ``(x, y, z)`` of the most extreme point, or
        ``None`` if the orientation matrix cannot be retrieved.

    Note:
        ``region`` must contain binary values (1 for foreground, 0 for
        background).  Zoom scaling from ``poi`` is applied to ensure
        physically meaningful distances.
    """
    direction_: Sequence[DIRECTIONS] | Sequence[tuple[DIRECTIONS, float]] = direction if isinstance(direction, Sequence) else (direction,)  # type: ignore

    direction_weight: list[tuple[DIRECTIONS, float]] = [(d[0], 1.0) if isinstance(d, str) else (d[0], d[1]) for d in direction_]  # type: ignore
    assert poi.orientation == region.orientation, (poi.orientation, region.orientation)
    try:
        to_reference_frame, from_reference_frame = get_vert_direction_matrix(poi, vert_id=vert_id, to_pir=False)
    except KeyError:
        return None
    pc = np.stack(np.where(region.get_array() == 1))
    cords = to_reference_frame @ pc  # 3,n; 3 = P,I,R of vert
    a = [_get_sub_array_by_direction(d, cords) * poi.zoom[poi.get_axis(d)] * w for d, w in direction_weight]

    idx = np.argmax(sum(a))
    return pc[:, idx]

project_pois_onto_segmentation_surface ¶

project_pois_onto_segmentation_surface(poi: POI, seg: NII, connectivity: int = 1, dilated_surface: bool = False) -> POI

Projects points of interest (POI) onto a segmentation surface.

This function computes the surface points of a segmentation volume and projects the given points of interest (POI) onto these computed surface points.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	The points of interest to be projected.	required
`seg`	`NII`	A segmentation volume object containing the target surface.	required
`connectivity`	`int`	The connectivity level for defining the surface. Default is 1, where 1 denotes face-connectivity.	`1`
`dilated_surface`	`bool`	Whether to compute a dilated version of the surface, expanding the surface area. Default is False.	`False`

Returns:

Name	Type	Description
`POI`	`POI`	The points of interest projected onto the segmentation surface.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def project_pois_onto_segmentation_surface(
    poi: POI,
    seg: NII,
    connectivity: int = 1,
    dilated_surface: bool = False,
) -> POI:
    """Projects points of interest (POI) onto a segmentation surface.

    This function computes the surface points of a segmentation volume and
    projects the given points of interest (POI) onto these computed surface points.

    Args:
        poi (POI): The points of interest to be projected.
        seg (NII): A segmentation volume object containing the target surface.
        connectivity (int, optional): The connectivity level for defining the surface.
            Default is 1, where 1 denotes face-connectivity.
        dilated_surface (bool, optional): Whether to compute a dilated version of the
            surface, expanding the surface area. Default is False.

    Returns:
        POI: The points of interest projected onto the segmentation surface.
    """
    point_set = seg.compute_surface_points(
        connectivity=connectivity,
        dilated_surface=dilated_surface,
    )
    return project_pois_onto_set_of_points(poi, point_set)

project_pois_onto_set_of_points ¶

project_pois_onto_set_of_points(poi: POI, point_set: list[COORDINATE]) -> POI

Projects points of interest (POI) onto the nearest points in a given set.

For each point in the POI, this function finds the closest point in the provided point set and updates the POI coordinates to align with the nearest points.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	The points of interest to be projected.	required
`point_set`	`list[COORDINATE]`	A list of coordinates representing the target points for projection.	required

Returns:

Name	Type	Description
`POI`	`POI`	The updated POI with coordinates projected onto the nearest points in
	`POI`	the provided point set.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def project_pois_onto_set_of_points(poi: POI, point_set: list[COORDINATE]) -> POI:
    """Projects points of interest (POI) onto the nearest points in a given set.

    For each point in the POI, this function finds the closest point in the
    provided point set and updates the POI coordinates to align with the nearest points.

    Args:
        poi (POI): The points of interest to be projected.
        point_set (list[COORDINATE]): A list of coordinates representing the target points
            for projection.

    Returns:
        POI: The updated POI with coordinates projected onto the nearest points in
        the provided point set.
    """
    poi_n = poi.copy()
    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

    return poi_n

cdist_to_point ¶

cdist_to_point(point: COORDINATE | ndarray, a: ndarray) -> np.ndarray

Compute Euclidean distances from a single point to every row in an array.

Parameters:

Name	Type	Description	Default
`point`	`COORDINATE \| ndarray`	Query point as a sequence or 1-D array of length 3.	required
`a`	`ndarray`	Array of shape `(N, 3)` containing the candidate points.	required

Returns:

Type	Description
`ndarray`	1-D float array of shape `(N,)` with the Euclidean distance from
`ndarray`	`point` to each row of `a`.

Source code in TPTBox/core/poi_fun/pixel_based_point_finder.py

def cdist_to_point(
    point: COORDINATE | np.ndarray,
    a: np.ndarray,
) -> np.ndarray:
    """Compute Euclidean distances from a single point to every row in an array.

    Args:
        point: Query point as a sequence or 1-D array of length 3.
        a: Array of shape ``(N, 3)`` containing the candidate points.

    Returns:
        1-D float array of shape ``(N,)`` with the Euclidean distance from
        ``point`` to each row of ``a``.
    """
    return cdist([point], a)[0]

Strategies¶

TPTBox.core.poi_fun.strategies ¶

strategy_extreme_points ¶

strategy_extreme_points(poi: POI, current_subreg: NII, location: Location, direction: Sequence[DIRECTIONS] | DIRECTIONS, vert_id: int, subreg_id: Location | list[Location], bb: tuple[slice, slice, slice] | None, log: Logger_Interface = _log) -> None

Strategy: place a POI at the anatomically extreme voxel of a subregion.

Extracts the label(s) subreg_id from current_subreg, finds the voxel most extreme in direction relative to the vertebra orientation, and stores the result in poi.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will receive the new landmark at `(vert_id, location.value)`.	required
`current_subreg`	`NII`	Cropped subregion `NII` aligned to the vertebra bounding box.	required
`location`	`Location`	Target `Location` enum member identifying the landmark slot.	required
`direction`	`Sequence[DIRECTIONS] \| DIRECTIONS`	Anatomical direction(s) to optimise. Accepts a single direction string, a sequence of strings, or weighted `(direction, weight)` tuples.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`subreg_id`	`Location \| list[Location]`	Label(s) to extract from `current_subreg` before searching for the extreme point.	required
`bb`	`tuple[slice, slice, slice] \| None`	Bounding-box slices mapping local to global voxel coordinates.	required
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`

Source code in TPTBox/core/poi_fun/strategies.py

def strategy_extreme_points(
    poi: POI,
    current_subreg: NII,
    location: Location,
    direction: Sequence[DIRECTIONS] | DIRECTIONS,
    vert_id: int,
    subreg_id: Location | list[Location],
    bb: tuple[slice, slice, slice] | None,
    log: Logger_Interface = _log,
) -> None:
    """Strategy: place a POI at the anatomically extreme voxel of a subregion.

    Extracts the label(s) ``subreg_id`` from ``current_subreg``, finds the
    voxel most extreme in ``direction`` relative to the vertebra orientation,
    and stores the result in ``poi``.

    Args:
        poi: ``POI`` object that will receive the new landmark at
            ``(vert_id, location.value)``.
        current_subreg: Cropped subregion ``NII`` aligned to the vertebra
            bounding box.
        location: Target ``Location`` enum member identifying the landmark slot.
        direction: Anatomical direction(s) to optimise.  Accepts a single
            direction string, a sequence of strings, or weighted
            ``(direction, weight)`` tuples.
        vert_id: Vertebra identifier (integer label).
        subreg_id: Label(s) to extract from ``current_subreg`` before
            searching for the extreme point.
        bb: Bounding-box slices mapping local to global voxel coordinates.
        log: Logger for warning and error messages.
    """
    region = current_subreg.extract_label(subreg_id)
    if region.sum() == 0:
        log.on_fail(f"reg={vert_id},subreg={subreg_id} is missing (extreme_points); {current_subreg.unique()}")
        return
    # extreme_point = get_extreme_point(poi, region, vert_id, bb, anti_point)

    extreme_point = get_extreme_point_by_vert_direction(poi, region, vert_id, direction)
    if extreme_point is None:
        return
    poi[vert_id, location.value] = tuple(a.start + b for a, b in zip_strict(bb, extreme_point))

strategy_line_cast ¶

strategy_line_cast(poi: POI, vert_id: int, current_subreg: NII, location: Location, start_point: Location | ndarray, regions_loc: list[Location] | Location, normal_vector_points: tuple[Location, Location] | DIRECTIONS, bb: tuple[slice, slice, slice] | None, log: Logger_Interface = _log) -> None

Strategy: place a POI at the ray-cast exit point of a subregion surface.

Extracts the combined label mask for regions_loc, then casts a ray from start_point along normal_vector_points using :func:~TPTBox.core.poi_fun.ray_casting.max_distance_ray_cast_convex_poi and stores the resulting boundary coordinate in poi.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will receive the new landmark at `(vert_id, location.value)`.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`current_subreg`	`NII`	Cropped subregion `NII` aligned to the vertebra bounding box.	required
`location`	`Location`	Target `Location` enum member identifying the landmark slot.	required
`start_point`	`Location \| ndarray`	Ray origin — a `Location` member resolved from `poi` or a pre-computed numpy coordinate.	required
`regions_loc`	`list[Location] \| Location`	Label(s) to extract from `current_subreg` as the region to cast the ray into.	required
`normal_vector_points`	`tuple[Location, Location] \| DIRECTIONS`	Ray direction, given either as a `DIRECTIONS` string or a pair of `Location` members whose difference defines the direction.	required
`bb`	`tuple[slice, slice, slice] \| None`	Bounding-box slices mapping local to global voxel coordinates.	required
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`

Source code in TPTBox/core/poi_fun/strategies.py

def strategy_line_cast(
    poi: POI,
    vert_id: int,
    current_subreg: NII,
    location: Location,
    start_point: Location | np.ndarray,
    regions_loc: list[Location] | Location,
    normal_vector_points: tuple[Location, Location] | DIRECTIONS,
    bb: tuple[slice, slice, slice] | None,
    log: Logger_Interface = _log,
) -> None:
    """Strategy: place a POI at the ray-cast exit point of a subregion surface.

    Extracts the combined label mask for ``regions_loc``, then casts a ray
    from ``start_point`` along ``normal_vector_points`` using
    :func:`~TPTBox.core.poi_fun.ray_casting.max_distance_ray_cast_convex_poi`
    and stores the resulting boundary coordinate in ``poi``.

    Args:
        poi: ``POI`` object that will receive the new landmark at
            ``(vert_id, location.value)``.
        vert_id: Vertebra identifier (integer label).
        current_subreg: Cropped subregion ``NII`` aligned to the vertebra
            bounding box.
        location: Target ``Location`` enum member identifying the landmark slot.
        start_point: Ray origin — a ``Location`` member resolved from ``poi`` or
            a pre-computed numpy coordinate.
        regions_loc: Label(s) to extract from ``current_subreg`` as the region
            to cast the ray into.
        normal_vector_points: Ray direction, given either as a ``DIRECTIONS``
            string or a pair of ``Location`` members whose difference defines
            the direction.
        bb: Bounding-box slices mapping local to global voxel coordinates.
        log: Logger for warning and error messages.
    """
    region = current_subreg.extract_label(regions_loc)
    # if legacy_code:
    #    horizontal_plane_landmarks_old(poi, region, label_id, bb, log)
    # else:
    extreme_point = max_distance_ray_cast_convex_poi(poi, region, vert_id, bb, normal_vector_points, start_point, log=log)
    if extreme_point is None:
        return
    poi[vert_id, location.value] = tuple(a.start + b for a, b in zip_strict(bb, extreme_point))

strategy_find_corner ¶

strategy_find_corner(poi: POI, current_subreg: NII, vert_id: int, bb: tuple[slice, slice, slice], location: Location, vec1: Location, vec2: Location, start_point: Location | ndarray = Location.Vertebra_Corpus, log: Logger_Interface = _log, shift_direction: DIRECTIONS | None = None) -> None

Strategy: place a POI at the corner of a vertebral body defined by two direction vectors.

Shifts the start point laterally if shift_direction is given, then uses a 2-D bisection search within the vertebral body (Vertebra_Corpus and Vertebra_Corpus_border) to find the furthest reachable corner in the directions of vec1 and vec2.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will receive the new landmark at `(vert_id, location.value)`.	required
`current_subreg`	`NII`	Cropped subregion `NII` aligned to the vertebra bounding box.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice]`	Bounding-box slices mapping local to global voxel coordinates.	required
`location`	`Location`	Target `Location` enum member identifying the landmark slot.	required
`vec1`	`Location`	First direction `Location` (e.g. anterior or superior edge point).	required
`vec2`	`Location`	Second direction `Location` (e.g. superior or inferior edge point).	required
`start_point`	`Location \| ndarray`	Ray origin — a `Location` member resolved from `poi` or a pre-computed numpy coordinate. Defaults to `Location.Vertebra_Corpus`.	`Vertebra_Corpus`
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`
`shift_direction`	`DIRECTIONS \| None`	Optional lateral shift direction applied to `start_point` before the corner search (e.g. `"R"` for right). `None` skips the shift.	`None`

Source code in TPTBox/core/poi_fun/strategies.py

def strategy_find_corner(
    poi: POI,
    current_subreg: NII,
    vert_id: int,
    bb: tuple[slice, slice, slice],
    location: Location,
    vec1: Location,
    vec2: Location,
    start_point: Location | np.ndarray = Location.Vertebra_Corpus,
    log: Logger_Interface = _log,
    shift_direction: DIRECTIONS | None = None,
) -> None:
    """Strategy: place a POI at the corner of a vertebral body defined by two direction vectors.

    Shifts the start point laterally if ``shift_direction`` is given, then
    uses a 2-D bisection search within the vertebral body (``Vertebra_Corpus``
    and ``Vertebra_Corpus_border``) to find the furthest reachable corner in
    the directions of ``vec1`` and ``vec2``.

    Args:
        poi: ``POI`` object that will receive the new landmark at
            ``(vert_id, location.value)``.
        current_subreg: Cropped subregion ``NII`` aligned to the vertebra
            bounding box.
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices mapping local to global voxel coordinates.
        location: Target ``Location`` enum member identifying the landmark slot.
        vec1: First direction ``Location`` (e.g. anterior or superior edge
            point).
        vec2: Second direction ``Location`` (e.g. superior or inferior edge
            point).
        start_point: Ray origin — a ``Location`` member resolved from ``poi``
            or a pre-computed numpy coordinate.  Defaults to
            ``Location.Vertebra_Corpus``.
        log: Logger for warning and error messages.
        shift_direction: Optional lateral shift direction applied to
            ``start_point`` before the corner search (e.g. ``"R"`` for right).
            ``None`` skips the shift.
    """
    if vert_id in sacrum_w_o_arcus:
        return

    start_point = shift_point(poi, vert_id, bb, start_point, direction=shift_direction, log=log)
    location2 = [Location.Vertebra_Corpus, Location.Vertebra_Corpus_border]
    current_subreg = current_subreg.extract_label(location2, keep_label=False)
    if current_subreg.sum() == 0:
        return
    corner_point = _find_corner_point(poi, current_subreg, vert_id, bb, start_point, vec1=vec1, vec2=vec2, log=log, location=location)

    if corner_point is None:
        return
    poi[vert_id, location.value] = tuple(a.start + b for a, b in zip_strict(bb, corner_point))

strategy_ligament_attachment_point_flava ¶

strategy_ligament_attachment_point_flava(poi: POI, current_subreg: NII, vert_id: int, bb: tuple[slice, slice, slice], location: Location, goal: Location | ndarray, log: Logger_Interface = _log, delta: float = 1e-07, shift_direction: Literal['S', 'I'] = 'S', dir2: Literal['A'] = 'A') -> None

Strategy: place a ligamentum-flavum attachment POI on the arcus / spinosus surface.

Starting from the spinosus process (or arcus if spinosus is absent), performs a 2-D bisection search moving in shift_direction and dir2 while staying inside the arcus/spinosus mask and approaching goal. The search stops when the bracket width falls below delta.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will receive the new landmark at `(vert_id, location)`.	required
`current_subreg`	`NII`	Cropped subregion `NII` aligned to the vertebra bounding box (must contain arcus and spinosus labels).	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice]`	Bounding-box slices mapping local to global voxel coordinates.	required
`location`	`Location`	Target `Location` enum member identifying the landmark slot.	required
`goal`	`Location \| ndarray`	Target `Location` or pre-computed coordinate toward which the search advances (e.g. an anterior longitudinal ligament point).	required
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`
`delta`	`float`	Convergence threshold for the bisection loop. Defaults to 0.0000001.	`1e-07`
`shift_direction`	`Literal['S', 'I']`	Primary direction of search along the arcus (`"S"` = superior, `"I"` = inferior). Defaults to `"S"`.	`'S'`
`dir2`	`Literal['A']`	Secondary search direction (always anterior, `"A"`). Defaults to `"A"`.	`'A'`

Source code in TPTBox/core/poi_fun/strategies.py

def strategy_ligament_attachment_point_flava(
    poi: POI,
    current_subreg: NII,
    vert_id: int,
    bb: tuple[slice, slice, slice],
    location: Location,
    goal: Location | np.ndarray,
    log: Logger_Interface = _log,
    delta: float = 0.0000001,
    shift_direction: Literal["S", "I"] = "S",
    dir2: Literal["A"] = "A",
) -> None:
    """Strategy: place a ligamentum-flavum attachment POI on the arcus / spinosus surface.

    Starting from the spinosus process (or arcus if spinosus is absent),
    performs a 2-D bisection search moving in ``shift_direction`` and ``dir2``
    while staying inside the arcus/spinosus mask and approaching ``goal``.
    The search stops when the bracket width falls below ``delta``.

    Args:
        poi: ``POI`` object that will receive the new landmark at
            ``(vert_id, location)``.
        current_subreg: Cropped subregion ``NII`` aligned to the vertebra
            bounding box (must contain arcus and spinosus labels).
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices mapping local to global voxel coordinates.
        location: Target ``Location`` enum member identifying the landmark slot.
        goal: Target ``Location`` or pre-computed coordinate toward which the
            search advances (e.g. an anterior longitudinal ligament point).
        log: Logger for warning and error messages.
        delta: Convergence threshold for the bisection loop.
            Defaults to 0.0000001.
        shift_direction: Primary direction of search along the arcus
            (``"S"`` = superior, ``"I"`` = inferior).  Defaults to ``"S"``.
        dir2: Secondary search direction (always anterior, ``"A"``).
            Defaults to ``"A"``.
    """
    if vert_id in sacrum_w_o_arcus:
        return
    try:
        normal_vector1 = get_direction(shift_direction, poi, vert_id)  # / np.array(poi.zoom)
        v1 = normal_vector1 / norm(normal_vector1)

        normal_vector2 = get_direction(dir2, poi, vert_id)  # / np.array(poi.zoom)
        v2 = normal_vector2 / norm(normal_vector2)
    except KeyError:
        return
    if isinstance(goal, Location):
        start_point_np = to_local_np(Location.Spinosus_Process, bb, poi, vert_id, log, verbose=False)
        if start_point_np is None:
            start_point_np = to_local_np(Location.Arcus_Vertebrae, bb, poi, vert_id, log, verbose=True)
    else:
        start_point_np = goal

    goal_np = to_local_np(goal, bb, poi, vert_id, log) if isinstance(goal, Location) else goal
    if goal_np is None or start_point_np is None:
        return

    region = current_subreg.extract_label([Location.Arcus_Vertebrae, Location.Spinosus_Process]).get_array()
    interpolator = RegularGridInterpolator([np.arange(region.shape[i]) for i in range(3)], region)
    dist_curr = norm(start_point_np - goal_np)

    def is_inside_and_closer(f1=0.0, f2=0.0):
        coords = start_point_np + f1 * v1 + f2 * v2
        if any(i < 0 for i in coords) or any(coords[i] > region.shape[i] - 1 for i in range(len(coords))):
            return False
        if interpolator(coords) > 0.5:
            dist_new = norm(coords - goal_np)
            return dist_curr > dist_new
        else:
            return False

    v1_n: float = 0.5
    v2_n: float = 0.5
    f1 = f2 = 0
    # Refine factors using delta
    while delta < (v1_n + v2_n) / 2:
        changed = False
        if is_inside_and_closer(v1_n + f1, f2):
            f1 += v1_n
            changed = True
            coords = start_point_np + f1 * v1 + f2 * v2
            dist_curr = norm(coords - goal_np)
        if is_inside_and_closer(f1, f2 + v2_n):
            f2 += v2_n
            changed = True
            coords = start_point_np + f1 * v1 + f2 * v2
            dist_curr = norm(coords - goal_np)

        if not changed:
            v1_n /= 2
            v2_n /= 2

    coords = start_point_np + f1 * v1 + f2 * v2
    poi[vert_id, location] = tuple(x + y.start for x, y in zip(coords, bb))

strategy_shifted_line_cast ¶

strategy_shifted_line_cast(poi: POI, current_subreg: NII, location: Location, vert_id: int, bb: tuple[slice, slice, slice] | None, regions_loc: list[Location] | Location, normal_vector_points: tuple[Location, Location] | DIRECTIONS, start_point: Location = Location.Vertebra_Corpus, log: Logger_Interface = _log, direction: DIRECTIONS = 'R', do_shift: bool = True) -> None

Strategy: lateral-shift a start point then ray-cast to the subregion surface.

Optionally displaces start_point laterally by a vertebra-width fraction via :func:~TPTBox.core.poi_fun.ray_casting.shift_point, then delegates to :func:strategy_line_cast to find the boundary along normal_vector_points.

Parameters:

Name	Type	Description	Default
`poi`	`POI`	`POI` object that will receive the new landmark at `(vert_id, location.value)`.	required
`current_subreg`	`NII`	Cropped subregion `NII` aligned to the vertebra bounding box.	required
`location`	`Location`	Target `Location` enum member identifying the landmark slot.	required
`vert_id`	`int`	Vertebra identifier (integer label).	required
`bb`	`tuple[slice, slice, slice] \| None`	Bounding-box slices mapping local to global voxel coordinates.	required
`regions_loc`	`list[Location] \| Location`	Label(s) to extract from `current_subreg` as the target region for the ray cast.	required
`normal_vector_points`	`tuple[Location, Location] \| DIRECTIONS`	Ray direction — a `DIRECTIONS` string or a pair of `Location` members.	required
`start_point`	`Location`	Ray origin before shifting. Defaults to `Location.Vertebra_Corpus`.	`Vertebra_Corpus`
`log`	`Logger_Interface`	Logger for warning and error messages.	`_log`
`direction`	`DIRECTIONS`	Lateral direction for the shift (`"R"` or `"L"`). Defaults to `"R"`.	`'R'`
`do_shift`	`bool`	When `False`, skip the lateral shift and use `start_point` directly. Defaults to `True`.	`True`

Source code in TPTBox/core/poi_fun/strategies.py

def strategy_shifted_line_cast(
    poi: POI,
    current_subreg: NII,
    location: Location,
    vert_id: int,
    bb: tuple[slice, slice, slice] | None,
    regions_loc: list[Location] | Location,
    normal_vector_points: tuple[Location, Location] | DIRECTIONS,
    start_point: Location = Location.Vertebra_Corpus,
    log: Logger_Interface = _log,
    direction: DIRECTIONS = "R",
    do_shift: bool = True,
) -> None:
    """Strategy: lateral-shift a start point then ray-cast to the subregion surface.

    Optionally displaces ``start_point`` laterally by a vertebra-width fraction
    via :func:`~TPTBox.core.poi_fun.ray_casting.shift_point`, then delegates to
    :func:`strategy_line_cast` to find the boundary along ``normal_vector_points``.

    Args:
        poi: ``POI`` object that will receive the new landmark at
            ``(vert_id, location.value)``.
        current_subreg: Cropped subregion ``NII`` aligned to the vertebra
            bounding box.
        location: Target ``Location`` enum member identifying the landmark slot.
        vert_id: Vertebra identifier (integer label).
        bb: Bounding-box slices mapping local to global voxel coordinates.
        regions_loc: Label(s) to extract from ``current_subreg`` as the target
            region for the ray cast.
        normal_vector_points: Ray direction — a ``DIRECTIONS`` string or a pair
            of ``Location`` members.
        start_point: Ray origin before shifting.  Defaults to
            ``Location.Vertebra_Corpus``.
        log: Logger for warning and error messages.
        direction: Lateral direction for the shift (``"R"`` or ``"L"``).
            Defaults to ``"R"``.
        do_shift: When ``False``, skip the lateral shift and use ``start_point``
            directly.  Defaults to ``True``.
    """
    if vert_id in sacrum_w_o_arcus:
        return
    try:
        cords = shift_point(
            poi,
            vert_id,
            bb,
            start_point,
            direction=direction if do_shift else None,
            log=log,
        )
    except KeyError:
        log.on_fail(f"region={vert_id},subregion={location} is missing. (shifted_line_cast)")
        return
    if cords is None:
        return
    strategy_line_cast(
        start_point=cords,
        regions_loc=regions_loc,
        normal_vector_points=normal_vector_points,
        poi=poi,
        vert_id=vert_id,
        current_subreg=current_subreg,
        location=location,
        bb=bb,
        log=_log,
    )

Save / Load¶

TPTBox.core.poi_fun.save_load ¶

save_poi ¶

save_poi(poi: POI | POI_Global, out_path: Path | str, make_parents=False, additional_info: dict | None = None, resample_reference: Has_Grid | None = None, verbose: logging = True, save_hint=2) -> None

Saves the POIs to a JSON file.

Parameters:

Name	Type	Description	Default
`out_path`	`Path \| str`	The path where the JSON file will be saved.	required
`make_parents`	`bool`	If True, create any necessary parent directories for the output file. Defaults to False.	`False`
`verbose`	`bool`	If True, print status messages to the console. Defaults to True.	`True`
`save_hint`		0 Default, 1 Gruber, 2 POI (readable), 10 ISO-POI (outdated)	`2`

Returns:

Type	Description
`None`	None

Raises:

Type	Description
`TypeError`	If any of the POIs have an invalid type.

Example

POIs = Centroids(...) POIs.save("output/POIs.json")

Source code in TPTBox/core/poi_fun/save_load.py

def save_poi(
    poi: POI | POI_Global,
    out_path: Path | str,
    make_parents=False,
    additional_info: dict | None = None,
    resample_reference: Has_Grid | None = None,
    verbose: logging = True,
    save_hint=2,
) -> None:
    """Saves the POIs to a JSON file.

    Args:
        out_path (Path | str): The path where the JSON file will be saved.
        make_parents (bool, optional): If True, create any necessary parent directories for the output file.
            Defaults to False.
        verbose (bool, optional): If True, print status messages to the console. Defaults to True.
        save_hint: 0 Default, 1 Gruber, 2 POI (readable), 10 ISO-POI (outdated)

    Returns:
        None

    Raises:
        TypeError: If any of the POIs have an invalid type.

    Example:
        >>> POIs = Centroids(...)
        >>> POIs.save("output/POIs.json")
    """
    _file_types = ["json"]
    file_ending = Path(out_path).name.split(".")[-1]
    if file_ending not in _file_types:
        raise ValueError(f"Not supported file ending for POI: {file_ending} not in {_file_types}")
    if make_parents:
        Path(out_path).parent.mkdir(exist_ok=True, parents=True)

    poi.sort()
    out_path = str(out_path)
    if len(poi.centroids) == 0:
        log.print("POIs empty, not saved:", out_path, ltype=Log_Type.FAIL, verbose=verbose)
        return
    json_object, print_add = _poi_to_dict_list(poi, additional_info, save_hint, resample_reference, verbose)

    # Problem with python 3 and int64 serialization.
    def convert(o):
        if isinstance(o, np.integer):
            return int(o)
        if isinstance(o, np.floating):
            return float(o)
        if isinstance(o, np.ndarray):
            return o.tolist()
        if isinstance(o, Path):
            return str(o.absolute())
        raise TypeError(type(o))

    try:
        with open(out_path, "w") as f:
            json.dump(json_object, f, default=convert, indent=4)
    except TypeError:
        Path(out_path).unlink(missing_ok=True)
        raise
    log.print("POIs saved:", out_path, print_add, ltype=Log_Type.SAVE, verbose=verbose)

load_poi ¶

load_poi(ctd_path: POI_Reference, verbose=True) -> POI | POI_Global

Load POIs from a file or a BIDS file object.

Parameters:

Name	Type	Description	Default
`ctd_path`	`Centroid_Reference`	Path to a file or BIDS file object from which to load POIs. Alternatively, it can be a tuple containing the following items: - vert: str, the name of the vertebra. - subreg: str, the name of the subregion. - ids: list[int \| Location], a list of integers and/or Location objects used to filter the POIs.	required

Returns:

Type	Description
`POI \| POI_Global`	A Centroids object containing the loaded POIs.

Raises:

Type	Description
`AssertionError`	If `ctd_path` is not a recognized type.

Source code in TPTBox/core/poi_fun/save_load.py

def load_poi(ctd_path: POI_Reference, verbose=True) -> POI | POI_Global:  # noqa: ARG001
    """Load POIs from a file or a BIDS file object.

    Args:
        ctd_path (Centroid_Reference): Path to a file or BIDS file object from which to load POIs.
            Alternatively, it can be a tuple containing the following items:
            - vert: str, the name of the vertebra.
            - subreg: str, the name of the subregion.
            - ids: list[int | Location], a list of integers and/or Location objects used to filter the POIs.

    Returns:
        A Centroids object containing the loaded POIs.

    Raises:
        AssertionError: If `ctd_path` is not a recognized type.

    """
    from TPTBox import POI, calc_poi_from_subreg_vert

    # already loaded
    if isinstance(ctd_path, POI):
        return ctd_path
    if isinstance(ctd_path, tuple):
        vert = ctd_path[0]
        subreg = ctd_path[1]
        ids: list[int | Location] = ctd_path[2]  # type: ignore
        return calc_poi_from_subreg_vert(vert, subreg, subreg_id=ids)

    ### Check Datatype ###
    dict_list = _open_file(ctd_path)
    ### New Spine_r has a dict instead of a dict list. ###
    if isinstance(dict_list, dict):
        if "markups" in dict_list:
            return _load_mkr_POI(dict_list)
        else:
            return _load_form_POI_spine_r2(dict_list)
    ### format_POI_old has no META header ###
    if "direction" not in dict_list[0] and "vert_label" in dict_list[0]:
        return _load_format_POI_old(dict_list)  # This file if used in the old POI-pipeline and is deprecated
    ### Ours Global Space ##
    format_ = dict_list[0].get("format", None)
    format_ = format_key2value[format_] if format_ is not None else None
    level_one_info = _register_lvl[dict_list[0].get("level_one_info", Vertebra_Instance.__name__)]
    level_two_info = _register_lvl[dict_list[0].get("level_two_info", Location.__name__)]
    info = {k: v for k, v in dict_list[0].items() if k not in ctd_info_blacklist}
    if format_ in (FORMAT_GLOBAL, FORMAT_PLST):
        from TPTBox import POI_Global

        centroids = POI_Descriptor()
        itk_coords = global_spacing_name_key2value[dict_list[0].get("coordinate_system", "nib")]
        _load_POI_centroids(dict_list, centroids, level_one_info, level_two_info)

        return POI_Global(centroids, itk_coords=itk_coords, level_one_info=level_one_info, level_two_info=level_two_info, info=info)

    ### Ours ###
    assert "direction" in dict_list[0] or format_ in [FORMAT_PLST], (
        f'File format error: first index must be a "Direction" but got {dict_list[0]}'
    )
    axcode: AX_CODES = tuple(dict_list[0].get("direction", ("U", "U", "U")))  # type: ignore
    zoom: ZOOMS = dict_list[0].get("zoom", None)  # type: ignore
    shape = dict_list[0].get("shape", None)  # type: ignore
    shape = tuple(shape) if shape is not None else None
    origin = dict_list[0].get("origin", None)
    origin = tuple(origin) if origin is not None else None
    rotation: ROTATION = dict_list[0].get("rotation", None)

    centroids = POI_Descriptor()
    if format_ in (FORMAT_DOCKER, FORMAT_GRUBER) or format_ is None:
        _load_docker_centroids(dict_list, centroids, format_)
    elif format_ in (FORMAT_POI,):
        _load_POI_centroids(dict_list, centroids, level_one_info, level_two_info)
    else:
        raise NotImplementedError(format_)
    return POI(
        centroids=centroids,
        orientation=axcode,
        zoom=zoom,
        shape=shape,  # type: ignore
        format=format_,
        info=info,
        origin=origin,  # type: ignore
        rotation=rotation,  # type: ignore
        level_one_info=level_one_info,
        level_two_info=level_two_info,
    )  # type: ignore