Spine

Spine-specific utilities: modular 2D snapshot generation and statistical measurements (distances, angles, intervertebral disc POIs).

2D Snapshots

TPTBox.spine.snapshot2D.snapshot_modular

Visualization_Type

Bases: Enum

Enum selecting how voxels along a projection axis are reduced to a 2D pixel value.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

class Visualization_Type(Enum):
    """Enum selecting how voxels along a projection axis are reduced to a 2D pixel value."""

    Slice = auto()
    Maximum_Intensity = auto()
    Maximum_Intensity_Colored_Depth = auto()
    Mean_Intensity = auto()

Snapshot_Frame dataclass

Dataclass bundling an image, optional segmentation/centroids, and per-frame rendering settings.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

@dataclass(init=True)
class Snapshot_Frame:
    """Dataclass bundling an image, optional segmentation/centroids, and per-frame rendering settings."""

    # Content
    image: Image_Reference
    segmentation: Image_Reference | None = None
    centroids: POI_Reference | None = None
    # Views
    sagittal: bool = True
    coronal: bool = False
    axial: bool = False

    axial_heights: list[float | int] | None = None
    # Title: str = for all views same, list entry for each view, None = no title
    title: str | list[str] | None = None
    # Image mode, cmap
    mode: Image_Modes = "MINMAX"
    cmap: ListedColormap | str = field(default_factory=lambda: cm_itk)
    alpha: float = 0.3
    # Pre-procesing
    crop_msk: bool = False
    crop_img: bool = False
    ignore_cdt_for_centering: bool = False
    rescale_to_iso: bool = True
    ignore_seg_for_centering: bool = False
    # Type, post-processing
    visualization_type: Visualization_Type = Visualization_Type.Slice
    only_mask_area: bool = False
    image_threshold: float | None = None  # everything below this threshold is set to min value of the img
    denoise_threshold: float | None = None  # threshold like above, but set for a smoothed img version
    gauss_filter: bool = False  # applies a gauss filter to the img
    cor_savgol_filter: bool = False  # applies a savgol_filter on the curve projection spline interpolation in coronal plane

    hide_segmentation: bool = False
    hide_centroids: bool = False
    hide_centroid_labels: bool = False
    poi_labelmap: dict[int, str] = field(default_factory=lambda: {k: v for k, v in v_idx2name.items() if k < 35})
    force_show_cdt: bool = False  # Shows the centroid computed by a segmentation, if no centroids are provided
    curve_location: Location | None = None  # Location.Vertebra_Corpus
    show_these_subreg_poi: list[int | Location] | None = None

sag_cor_curve_projection

sag_cor_curve_projection(ctd_list: POI, img_data: ndarray, ctd_fallback: POI, cor_savgol_filter: bool = False, curve_location: Location = Location.Vertebra_Corpus) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Makes a curve projection (spline interpolation) over the spline through the given centroids.

Parameters:

Name	Type	Description	Default
ctd_list	POI	given Centroids	required
img_data	ndarray	given img_data	required
cor_savgol_filter	bool	If true, will perform the savgol filter also in coronal view	False
curve_location	Location	Location of the curve's centroids to be used.	Vertebra_Corpus

Returns:

Name	Type	Description
x_ctd	tuple[ndarray, ndarray, ndarray]	ctd x values sorted (nparray), y_cord: interpolated y coords (nparray), z_cord: interpolated z coords (nparray)

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def sag_cor_curve_projection(
    ctd_list: POI,
    img_data: np.ndarray,
    ctd_fallback: POI,
    cor_savgol_filter: bool = False,
    curve_location: Location = Location.Vertebra_Corpus,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Makes a curve projection (spline interpolation) over the spline through the given centroids.

    Args:
        ctd_list: given Centroids
        img_data: given img_data
        cor_savgol_filter: If true, will perform the savgol filter also in coronal view
        curve_location: Location of the curve's centroids to be used.

    Returns:
        x_ctd: ctd x values sorted (nparray), y_cord: interpolated y coords (nparray), z_cord: interpolated z coords (nparray)
    """
    assert ctd_list is not None

    if 26 in ctd_list.centroids.keys_region() and cor_savgol_filter:
        warnings.warn(
            "Sacrum centroid present with cor_savgol_filter might overshadow the sacrum in coronal view",
            UserWarning,
            stacklevel=4,
        )
    # Sagittal and coronal projections of a curved plane defined by centroids
    # Note: Will assume IPL orientation!
    # if x-direction (=S/I) is not fully incremental, a straight, not an interpolated plane will be returned
    order = v_idx_order
    order += [i for i in range(256) if i not in v_idx_order]
    # ctd_list.sorting_list = v_idx_order
    ctd_list.round_(3)

    ctd_list.sort()
    if curve_location == Location.Vertebra_Corpus:
        s = ctd_list.keys_subregion()
        ids = 50 if 50 in s else (40 if 40 in s else 0)
        l = [v for k1, k2, v in ctd_list.items() if k2 == ids]
    else:
        l = [v for k1, k2, v in ctd_list.items() if k2 == curve_location.value]

    # TODO make the selected subregion flexible.
    if len(l) <= 3:
        l = list(ctd_list.values())
    if len(l) <= 3:
        l = list(ctd_fallback.values())
    # throw out all centroids that are not in correct up-down order
    x_cur = l[0][0] - 1
    throw_idx = []
    for idx, c in enumerate(l.copy()):
        if c[0] > x_cur:
            x_cur = c[0]
        else:
            throw_idx.append(idx)
    l = [i for idx, i in enumerate(l) if idx not in throw_idx] if len(l) - len(throw_idx) >= 3 else sorted(l, key=lambda x: x[0])

    ctd_arr = np.transpose(np.asarray(l))
    shp = img_data.shape
    x_ctd = np.rint(ctd_arr[0]).astype(int)
    y_ctd = np.rint(ctd_arr[1]).astype(int)
    z_ctd = np.rint(ctd_arr[2]).astype(int)
    # axl_plane = np.zeros((shp[1], shp[2]))
    try:
        f_sag = interp1d(x_ctd, y_ctd, kind="quadratic")
        f_cor = interp1d(x_ctd, z_ctd, kind="quadratic")
    except Exception:
        f_sag = interp1d(x_ctd, y_ctd, kind="linear")
        f_cor = interp1d(x_ctd, z_ctd, kind="linear")
        # print("quadratic", l, x_ctd, len(l), len(throw_idx))
        # exit()
    window_size = int((max(x_ctd) - min(x_ctd)) / 2)
    poly_order = 3
    if window_size % 2 == 0:
        window_size += 1
    window_size = max(window_size, poly_order + 1)
    y_cord = np.array([np.rint(f_sag(x)).astype(int) for x in range(min(x_ctd), max(x_ctd))])
    y_cord = np.rint(savgol_filter(y_cord, window_size, poly_order)).astype(int) if cor_savgol_filter else y_cord

    z_cord = np.array([np.rint(f_cor(x)).astype(int) for x in range(min(x_ctd), max(x_ctd))])
    z_cord = np.rint(savgol_filter(z_cord, window_size, poly_order)).astype(int)

    y_cord[y_cord < 0] = 0  # handle out-of-volume interpolations
    y_cord[y_cord >= shp[1]] = shp[1] - 1
    z_cord[z_cord < 0] = 0
    z_cord[z_cord >= shp[2]] = shp[2] - 1
    return x_ctd, y_cord, z_cord

curve_projected_slice

curve_projected_slice(x_ctd: ndarray, img_data: ndarray, y_cord: ndarray, z_cord: ndarray, axial_heights: list[float] | None) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Extract sagittal, coronal, and axial slice planes along a curved spinal path.

For each axial position the sagittal slice is taken at the interpolated coronal coordinate and vice versa. Outside the range of the centroids the boundary values are extrapolated.

Parameters:

Name	Type	Description	Default
x_ctd	ndarray	Sorted x-indices of the vertebral centroids (integer array).	required
img_data	ndarray	3-D image volume in IPL orientation with shape (X, Y, Z).	required
y_cord	ndarray	Interpolated y (anterior-posterior) coordinates for each x position between the first and last centroid.	required
z_cord	ndarray	Interpolated z (left-right) coordinates for each x position between the first and last centroid.	required
axial_heights	list[float] | None	Heights for the axial fallback planes passed through to curve_projection_axial_fallback.	required

Returns:

Type	Description
ndarray	A tuple (sag_plane, cor_plane, axl_plane) where each element is a
ndarray	2-D numpy array representing the sagittal, coronal, and axial
ndarray	projections respectively.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def curve_projected_slice(
    x_ctd: np.ndarray,
    img_data: np.ndarray,
    y_cord: np.ndarray,
    z_cord: np.ndarray,
    axial_heights: list[float] | None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Extract sagittal, coronal, and axial slice planes along a curved spinal path.

    For each axial position the sagittal slice is taken at the interpolated
    coronal coordinate and vice versa.  Outside the range of the centroids the
    boundary values are extrapolated.

    Args:
        x_ctd: Sorted x-indices of the vertebral centroids (integer array).
        img_data: 3-D image volume in IPL orientation with shape (X, Y, Z).
        y_cord: Interpolated y (anterior-posterior) coordinates for each x
            position between the first and last centroid.
        z_cord: Interpolated z (left-right) coordinates for each x position
            between the first and last centroid.
        axial_heights: Heights for the axial fallback planes passed through to
            ``curve_projection_axial_fallback``.

    Returns:
        A tuple ``(sag_plane, cor_plane, axl_plane)`` where each element is a
        2-D numpy array representing the sagittal, coronal, and axial
        projections respectively.
    """
    shp = img_data.shape
    cor_plane = np.zeros((shp[0], shp[2]))
    sag_plane = np.zeros((shp[0], shp[1]))
    for x in range(shp[0] - 1):
        if x < min(x_ctd):
            cor_plane[x, :] = img_data[x, y_cord[0], :]
            sag_plane[x, :] = img_data[x, :, z_cord[0]]
        elif x >= max(x_ctd):
            cor_plane[x, :] = img_data[x, y_cord[-1], :]
            sag_plane[x, :] = img_data[x, :, z_cord[-1]]
        else:
            cor_plane[x, :] = img_data[x, y_cord[x - min(x_ctd)], :]
            sag_plane[x, :] = img_data[x, :, z_cord[x - min(x_ctd)]]
    return (
        sag_plane,
        cor_plane,
        curve_projection_axial_fallback(img_data, x_ctd, heights=axial_heights),
    )

curve_projected_mean

curve_projected_mean(img_data: ndarray, zms: tuple[float, float, float], x_ctd: ndarray, y_cord: ndarray, ctd_list: POI, thick_t: tuple[int, int] = (100, 300), axial_heights: list[float] | None = None) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Compute mean-intensity curved-planar projections along the spinal curve.

For each axial position the sagittal mean is computed over all non-zero voxels in that row, while the coronal mean is computed over a slab of width thick_t centred on the interpolated y coordinate.

Parameters:

Name	Type	Description	Default
img_data	ndarray	3-D image volume in IPL orientation with shape (X, Y, Z).	required
zms	tuple[float, float, float]	Voxel spacing (zoom) in mm for each axis (x, y, z).	required
x_ctd	ndarray	Sorted x-indices of the vertebral centroids (integer array).	required
y_cord	ndarray	Interpolated y (anterior-posterior) coordinates for each x position between the first and last centroid.	required
ctd_list	POI	POI object whose region keys are used to adjust slab thickness near the sacrum.	required
thick_t	tuple[int, int]	Anterior and posterior slab half-widths in mm used for the coronal projection (anterior_mm, posterior_mm).	(100, 300)
axial_heights	list[float] | None	Heights for the axial fallback planes.	None

Returns:

Type	Description
tuple[ndarray, ndarray, ndarray]	A tuple (sag_plane, cor_plane, axl_plane) of 2-D float arrays.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def curve_projected_mean(
    img_data: np.ndarray,
    zms: tuple[float, float, float],
    x_ctd: np.ndarray,
    y_cord: np.ndarray,
    ctd_list: POI,
    thick_t: tuple[int, int] = (100, 300),
    axial_heights: list[float] | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Compute mean-intensity curved-planar projections along the spinal curve.

    For each axial position the sagittal mean is computed over all non-zero
    voxels in that row, while the coronal mean is computed over a slab of
    width ``thick_t`` centred on the interpolated y coordinate.

    Args:
        img_data: 3-D image volume in IPL orientation with shape (X, Y, Z).
        zms: Voxel spacing (zoom) in mm for each axis ``(x, y, z)``.
        x_ctd: Sorted x-indices of the vertebral centroids (integer array).
        y_cord: Interpolated y (anterior-posterior) coordinates for each x
            position between the first and last centroid.
        ctd_list: POI object whose region keys are used to adjust slab
            thickness near the sacrum.
        thick_t: Anterior and posterior slab half-widths in mm used for the
            coronal projection ``(anterior_mm, posterior_mm)``.
        axial_heights: Heights for the axial fallback planes.

    Returns:
        A tuple ``(sag_plane, cor_plane, axl_plane)`` of 2-D float arrays.
    """
    shp = img_data.shape
    cor_plane = np.zeros((shp[0], shp[2]))
    sag_plane = np.zeros((shp[0], shp[1]))
    y_zoom = zms[1]  # 0.9 = 1px = 0.9 mm # 10cm = 112px
    thick = (*thick_t,)

    for x in range(shp[0] - 1):
        if x < min(x_ctd):  # higher
            y_ref = y_cord[0]
        elif x >= max(x_ctd):  # lower than sacrum
            y_ref = y_cord[-1]
        else:
            y_ref = y_cord[x - min(x_ctd)]

        if 23 in ctd_list and x > int(ctd_list[23][1]):
            thick = (100, 50)

        thick = [int(i // y_zoom) + int(i % y_zoom > 0) for i in thick]
        y_post_rel_to_border = y_ref + int(0.4 * (shp[1] - 1 - y_ref))  # one-third distance to border
        y_range_low = int(max(0, y_ref - thick[1]))  # sagittal left
        y_range_high = int(min(y_ref + thick[0], y_post_rel_to_border))  # sagittal right
        cor_cut = img_data[x, y_range_low:y_range_high, :]

        plane_bool = np.zeros_like(cor_cut).astype(bool)
        plane_bool[cor_cut > 0] = True
        sag = np.nansum(img_data[x, :, :], 1, where=img_data[x, :, :] > 0)
        sag_plane[x, :] = div0(sag, np.count_nonzero(img_data[x, :, :], 1), fill=0)
        cor = np.nansum(cor_cut, 0, where=plane_bool)
        cor_plane[x, :] = div0(cor, np.count_nonzero(plane_bool, 0), fill=0)
    return (
        sag_plane,
        cor_plane,
        curve_projection_axial_fallback(img_data, x_ctd, heights=axial_heights),
    )

curve_projected_mip

curve_projected_mip(img_data: ndarray, zms: tuple[float, float, float], x_ctd: ndarray, y_cord: ndarray, ctd_list: POI, thick_t: tuple[int, int] = (100, 300), make_colored_depth: bool = False, axial_heights: list[float] | None = None) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Compute maximum-intensity curved-planar projections along the spinal curve.

A slab of configurable thickness is extracted around the interpolated spinal curve for each axial position, and the maximum intensity is projected onto the output plane. Optionally the depth of the maximum is colour-encoded using the inferno colourmap.

Parameters:

Name	Type	Description	Default
img_data	ndarray	3-D image volume in IPL orientation with shape (X, Y, Z).	required
zms	tuple[float, float, float]	Voxel spacing (zoom) in mm for each axis (x, y, z).	required
x_ctd	ndarray	Sorted x-indices of the vertebral centroids (integer array).	required
y_cord	ndarray	Interpolated y (anterior-posterior) coordinates for each x position between the first and last centroid.	required
ctd_list	POI	POI object (currently unused, reserved for future thickness adaptation near the sacrum).	required
thick_t	tuple[int, int]	Anterior and posterior slab half-widths in mm (anterior_mm, posterior_mm).	(100, 300)
make_colored_depth	bool	If True the returned planes are RGB arrays colour-coded by depth; otherwise they are scalar intensity arrays.	False
axial_heights	list[float] | None	Heights for the axial fallback planes.	None

Returns:

Type	Description
ndarray	A tuple (sag_plane, cor_plane, axl_plane). When
ndarray	make_colored_depth is True the sagittal and coronal arrays
ndarray	have shape (X, Y, 3) and (X, Z, 3) respectively.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def curve_projected_mip(
    img_data: np.ndarray,
    zms: tuple[float, float, float],
    x_ctd: np.ndarray,
    y_cord: np.ndarray,
    ctd_list: POI,  # noqa: ARG001
    thick_t: tuple[int, int] = (100, 300),
    make_colored_depth: bool = False,
    axial_heights: list[float] | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Compute maximum-intensity curved-planar projections along the spinal curve.

    A slab of configurable thickness is extracted around the interpolated
    spinal curve for each axial position, and the maximum intensity is
    projected onto the output plane.  Optionally the depth of the maximum is
    colour-encoded using the ``inferno`` colourmap.

    Args:
        img_data: 3-D image volume in IPL orientation with shape (X, Y, Z).
        zms: Voxel spacing (zoom) in mm for each axis ``(x, y, z)``.
        x_ctd: Sorted x-indices of the vertebral centroids (integer array).
        y_cord: Interpolated y (anterior-posterior) coordinates for each x
            position between the first and last centroid.
        ctd_list: POI object (currently unused, reserved for future thickness
            adaptation near the sacrum).
        thick_t: Anterior and posterior slab half-widths in mm
            ``(anterior_mm, posterior_mm)``.
        make_colored_depth: If ``True`` the returned planes are RGB arrays
            colour-coded by depth; otherwise they are scalar intensity arrays.
        axial_heights: Heights for the axial fallback planes.

    Returns:
        A tuple ``(sag_plane, cor_plane, axl_plane)``.  When
        ``make_colored_depth`` is ``True`` the sagittal and coronal arrays
        have shape ``(X, Y, 3)`` and ``(X, Z, 3)`` respectively.
    """
    shp = img_data.shape
    cor_plane = np.zeros((shp[0], shp[2]))
    cor_depth_plane = np.zeros((shp[0], shp[2]))
    sag_plane = np.zeros((shp[0], shp[1]))
    sag_depth_plane = np.zeros((shp[0], shp[1]))
    y_zoom = zms[1]  # 0.9 = 1px = 0.9 mm # 10cm = 112px
    thick = (*thick_t,)

    for x in range(shp[0] - 1):
        if x < min(x_ctd):  # higher
            y_ref = y_cord[0]
        elif x >= max(x_ctd):  # lower than sacrum
            y_ref = y_cord[-1]
        else:
            y_ref = y_cord[x - min(x_ctd)]

        # if 23 in ctd_list and x > int(ctd_list[23][1]) and not make_colored_depth:
        #    thick = (100, 50)

        # TODO set y_zoom for broken sample, see if it works
        try:
            thicke = [int(i // y_zoom) + int(i % y_zoom > 0) for i in thick]
        except Exception:
            print("thick infinity bug", y_zoom, thick_t, thick)
            thicke = (*thick_t,)
        thick = thicke
        y_post_rel_to_border = y_ref + int(0.4 * (shp[1] - 1 - y_ref))  # one-third distance to border
        y_range_low = int(max(0, y_ref - thick[1]))  # sagittal left
        y_range_high = int(min(y_ref + thick[0], y_post_rel_to_border))  # sagittal right
        # print("range", y_range_low, y_range_high)
        cor_cut = img_data[x, y_range_low:y_range_high, :]

        cor_plane[x, :] = np.max(cor_cut, axis=0)  # arr[x, mip_i, :]
        cor_depth_plane[x, :] = np.argmax(cor_cut, axis=0)
        sag_plane[x, :] = np.max(img_data[x, :, :], axis=1)  # img_data[x, :, z_ref]
        sag_depth_plane[x, :] = np.argmax(img_data[x, :, :], axis=1)

    if make_colored_depth:
        cor_depth_plane = normalize_image(cor_depth_plane)
        sag_depth_plane = normalize_image(sag_depth_plane)

        cor_plane = normalize_image(cor_plane)
        sag_plane = normalize_image(sag_plane)

        cor_m_plane = np.sqrt(cor_plane) * cor_depth_plane  # sqrt
        sag_m_plane = np.sqrt(sag_plane) * sag_depth_plane
        cor_m_plane = normalize_image(cor_m_plane)
        sag_m_plane = normalize_image(sag_m_plane)

        # convert to color image
        cmap = plt.get_cmap("inferno")
        cor_plane = cmap(cor_m_plane)[..., :3]
        sag_plane = cmap(sag_m_plane)[..., :3]

    return (
        sag_plane,
        cor_plane,
        curve_projection_axial_fallback(img_data, x_ctd, heights=axial_heights),
    )

normalize_image

normalize_image(img: ndarray, v_range: tuple[float, float] | None = None) -> np.ndarray

Linearly rescale an array to the range [0, 1].

Parameters:

Name	Type	Description	Default
img	ndarray	Input array to normalise.	required
v_range	tuple[float, float] | None	Optional (min, max) clip range. When None the actual minimum and maximum of img are used.	None

Returns:

Type	Description
ndarray	A float array with values in [0, 1].

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def normalize_image(img: np.ndarray, v_range: tuple[float, float] | None = None) -> np.ndarray:
    """Linearly rescale an array to the range [0, 1].

    Args:
        img: Input array to normalise.
        v_range: Optional ``(min, max)`` clip range.  When ``None`` the
            actual minimum and maximum of ``img`` are used.

    Returns:
        A float array with values in ``[0, 1]``.
    """
    if v_range is None:
        min_v = np.min(img)
        max_v = np.max(img)
    else:
        min_v = v_range[0]
        max_v = v_range[1]
    return (img - min_v) / (max_v - min_v)

curve_projection_axial_fallback

curve_projection_axial_fallback(img_data: ndarray, x_ctd: ndarray, heights: list[float] | None) -> np.ndarray

Build an axial composite plane for display as a fallback or from explicit heights.

When heights is given each value is interpreted as an absolute slice index (if |h| >= 1) or as a fraction of the volume depth (if |h| < 1). The slices are concatenated vertically. When heights is None three slices centred on the median centroid are concatenated.

Parameters:

Name	Type	Description	Default
img_data	ndarray	3-D image volume in IPL orientation with shape (X, Y, Z).	required
x_ctd	ndarray	Sorted x-indices of the vertebral centroids (integer array).	required
heights	list[float] | None	Explicit axial heights to extract, or None to use the three slices surrounding the central centroid.	required

Returns:

Type	Description
ndarray	A 2-D numpy array representing the composed axial view.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def curve_projection_axial_fallback(
    img_data: np.ndarray,
    x_ctd: np.ndarray,
    heights: list[float] | None,
) -> np.ndarray:
    """Build an axial composite plane for display as a fallback or from explicit heights.

    When ``heights`` is given each value is interpreted as an absolute slice
    index (if ``|h| >= 1``) or as a fraction of the volume depth (if
    ``|h| < 1``).  The slices are concatenated vertically.  When ``heights``
    is ``None`` three slices centred on the median centroid are concatenated.

    Args:
        img_data: 3-D image volume in IPL orientation with shape (X, Y, Z).
        x_ctd: Sorted x-indices of the vertebral centroids (integer array).
        heights: Explicit axial heights to extract, or ``None`` to use the
            three slices surrounding the central centroid.

    Returns:
        A 2-D numpy array representing the composed axial view.
    """
    if heights is not None:
        heights = [int(abs(h if abs(h) >= 1 else img_data.shape[0] * h)) for h in (heights) if abs(h) <= img_data.shape[0]]
        axl_plane = np.concatenate([img_data[h, :, :] for h in heights], axis=0)
        return axl_plane
    else:
        # Axial
        center = x_ctd[len(x_ctd) // 2]
        center_up = x_ctd[max(0, len(x_ctd) // 2 - 1)]
        center_down = x_ctd[min(len(x_ctd) - 1, len(x_ctd) // 2 + 1)]
        try:
            axl_plane = np.concatenate(
                [
                    img_data[(center + center_up) // 2, :, :],
                    img_data[center, :, :],
                    img_data[(center + center_down) // 2, :, :],
                ],
                axis=0,
            )
        except Exception as e:
            print(e)
            axl_plane = np.zeros((1, 1))
        return axl_plane

make_isotropic2d

make_isotropic2d(arr2d: ndarray, zms2d: tuple[float, float], msk: bool = False) -> np.ndarray

Resample a 2-D array to isotropic pixel spacing.

Parameters:

Name	Type	Description	Default
arr2d	ndarray	Input 2-D array.	required
zms2d	tuple[float, float]	Voxel spacing (dy, dz) of the input array in mm. The output will have 1 mm × 1 mm pixels.	required
msk	bool	If True nearest-neighbour interpolation is used (suitable for label maps); otherwise linear interpolation is applied.	False

Returns:

Type	Description
ndarray	Resampled 2-D array with isotropic pixel spacing.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def make_isotropic2d(arr2d: np.ndarray, zms2d: tuple[float, float], msk: bool = False) -> np.ndarray:
    """Resample a 2-D array to isotropic pixel spacing.

    Args:
        arr2d: Input 2-D array.
        zms2d: Voxel spacing ``(dy, dz)`` of the input array in mm.  The
            output will have 1 mm × 1 mm pixels.
        msk: If ``True`` nearest-neighbour interpolation is used (suitable
            for label maps); otherwise linear interpolation is applied.

    Returns:
        Resampled 2-D array with isotropic pixel spacing.
    """
    if np.issubdtype(arr2d.dtype, np.floating):
        arr2d = arr2d.astype(int)
    xs = list(range(arr2d.shape[0]))
    ys = list(range(arr2d.shape[1]))
    if msk:
        interpolator = RegularGridInterpolator((xs, ys), arr2d, method="nearest")
    else:
        interpolator = RegularGridInterpolator((xs, ys), arr2d, method="linear")
    new_shp = tuple(np.rint(np.multiply(arr2d.shape, zms2d)).astype(int))
    x_mm = np.linspace(0, arr2d.shape[0] - 1, num=new_shp[0])
    y_mm = np.linspace(0, arr2d.shape[1] - 1, num=new_shp[1])
    xx, yy = np.meshgrid(x_mm, y_mm)
    pts = np.vstack([xx.ravel(), yy.ravel()]).T
    try:
        lt = interpolator(pts)
    except Exception as e:
        raise ValueError(f"Needs to be casted into a other type: arr2d {arr2d.dtype}") from e  # noqa: B904
    img = np.reshape(lt, new_shp, order="F")
    return img

make_isotropic2dpluscolor

make_isotropic2dpluscolor(arr3d: ndarray, zms2d: tuple[float, float], msk: bool = False) -> np.ndarray

Resample a 2-D or 2-D+RGB array to isotropic pixel spacing.

Delegates to :func:make_isotropic2d for scalar images. For colour images (3rd dimension of size 3) each channel is resampled independently and the results are stacked.

Parameters:

Name	Type	Description	Default
arr3d	ndarray	Input array with shape (H, W) or (H, W, 3).	required
zms2d	tuple[float, float]	Voxel spacing (dy, dz) of the input array in mm.	required
msk	bool	Passed through to :func:make_isotropic2d; use True for label maps to enable nearest-neighbour interpolation.	False

Returns:

Type	Description
ndarray	Resampled array with the same number of channels as the input.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def make_isotropic2dpluscolor(
    arr3d: np.ndarray,
    zms2d: tuple[float, float],
    msk: bool = False,
) -> np.ndarray:
    """Resample a 2-D or 2-D+RGB array to isotropic pixel spacing.

    Delegates to :func:`make_isotropic2d` for scalar images.  For colour
    images (3rd dimension of size 3) each channel is resampled independently
    and the results are stacked.

    Args:
        arr3d: Input array with shape ``(H, W)`` or ``(H, W, 3)``.
        zms2d: Voxel spacing ``(dy, dz)`` of the input array in mm.
        msk: Passed through to :func:`make_isotropic2d`; use ``True`` for
            label maps to enable nearest-neighbour interpolation.

    Returns:
        Resampled array with the same number of channels as the input.
    """
    # print(arr3d.shape)
    if arr3d.ndim == 2:
        return make_isotropic2d(arr3d, zms2d, msk=msk)

    r_img = make_isotropic2d(arr3d[:, :, 0], zms2d, msk=msk)
    # print(r_img.shape)
    g_img = make_isotropic2d(arr3d[:, :, 1], zms2d, msk=msk)
    b_img = make_isotropic2d(arr3d[:, :, 2], zms2d, msk=msk)
    img = np.stack([r_img, g_img, b_img], axis=-1)
    # print(img.shape)
    return img

get_contrasting_stroke_color

get_contrasting_stroke_color(rgb: int | list[float] | tuple[float, ...]) -> str

Return "gray" or "black" depending on the perceived luminance of rgb.

A dark foreground colour (luminance < 0.3) gets a "gray" stroke so that text remains readable on dark backgrounds; all other colours receive a "black" stroke.

Parameters:

Name	Type	Description	Default
rgb	int | list[float] | tuple[float, ...]	Either an integer label index (looked up in the colour map) or an RGB / RGBA tuple / list of floats in [0, 1].	required

Returns:

Type	Description
str	"gray" for dark colours, "black" for bright colours.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def get_contrasting_stroke_color(rgb: int | list[float] | tuple[float, ...]) -> str:
    """Return ``"gray"`` or ``"black"`` depending on the perceived luminance of ``rgb``.

    A dark foreground colour (luminance < 0.3) gets a ``"gray"`` stroke so
    that text remains readable on dark backgrounds; all other colours receive
    a ``"black"`` stroke.

    Args:
        rgb: Either an integer label index (looked up in the colour map) or an
            RGB / RGBA tuple / list of floats in ``[0, 1]``.

    Returns:
        ``"gray"`` for dark colours, ``"black"`` for bright colours.
    """
    # Convert RGBA to RGB if necessary
    if isinstance(rgb, int):
        rgb = list(get_color_by_label(rgb).rgb / 255.0)
    if len(rgb) == 4:
        rgb = rgb[:3]
    luminance = 0.299 * rgb[0] + 0.587 * rgb[1] + 0.114 * rgb[2]
    return "gray" if luminance < 0.3 else "black"

create_figure

create_figure(dpi: int, planes: list[ndarray], has_title: bool = True) -> tuple

Create a matplotlib figure sized to contain all projection planes side by side.

Parameters:

Name	Type	Description	Default
dpi	int	Output resolution in dots per inch.	required
planes	list[ndarray]	List of 2-D (or 3-D colour) arrays whose widths determine the column layout of the figure.	required
has_title	bool	When True extra vertical space is reserved for a title row above the images.	True

Returns:

Type	Description
tuple	A (fig, axs) tuple where fig is the
tuple	class:~matplotlib.figure.Figure and axs is a list of
tuple	class:~matplotlib.axes.Axes, one per plane.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def create_figure(
    dpi: int,
    planes: list[np.ndarray],
    has_title: bool = True,
) -> tuple:
    """Create a matplotlib figure sized to contain all projection planes side by side.

    Args:
        dpi: Output resolution in dots per inch.
        planes: List of 2-D (or 3-D colour) arrays whose widths determine the
            column layout of the figure.
        has_title: When ``True`` extra vertical space is reserved for a title
            row above the images.

    Returns:
        A ``(fig, axs)`` tuple where ``fig`` is the
        :class:`~matplotlib.figure.Figure` and ``axs`` is a list of
        :class:`~matplotlib.axes.Axes`, one per plane.
    """
    fig_h = round(2 * planes[0].shape[0] / dpi, 2) + (0.5 if has_title else 0)
    plane_w = [p.shape[1] for p in planes]
    w = sum(plane_w)
    fig_w = round(2 * w / dpi, 2)
    x_pos = [0]
    for x in plane_w[:-1]:
        x_pos.append(x_pos[-1] + x)
    fig, axs = plt.subplots(1, len(planes), figsize=(fig_w, fig_h))

    if not isinstance(axs, np.ndarray):
        axs: list[Axes] = [axs]  # type: ignore
    for a in axs:
        a.axis("off")
        idx = list(axs).index(a)
        a.set_position([(x_pos[idx] / w), (-0.03 if has_title else 0), plane_w[idx] / w, 1])
    return fig, axs

plot_sag_centroids

plot_sag_centroids(axs: Axes, ctd: POI, zms: tuple[float, float, float], poi_labelmap: dict[int, str], hide_centroid_labels: bool, cmap: ListedColormap = cm_itk, curve_location: Location = Location.Vertebra_Corpus, show_these_subreg_poi: list[int | Location] | None = None) -> None

Overlay centroid circles and vertebra labels on a sagittal axis.

Points are drawn as filled circles coloured by vertebra label. Vertebra names from poi_labelmap are added as bold text next to the circle when hide_centroid_labels is False. Additional arrows and text annotations stored in ctd.info["line_segments_sag"] and ctd.info["text_sag"] are also rendered.

Parameters:

Name	Type	Description	Default
axs	Axes	Matplotlib axes on which to draw.	required
ctd	POI	POI object containing centroid coordinates and optional rendering hints in its info dictionary.	required
zms	tuple[float, float, float]	Voxel spacing (dx, dy, dz) in mm used to convert voxel coordinates to pixel coordinates.	required
poi_labelmap	dict[int, str]	Mapping from vertebra region ID to display name.	required
hide_centroid_labels	bool	When True vertebra name labels are omitted.	required
cmap	ListedColormap	Colour map used to assign per-vertebra colours.	cm_itk
curve_location	Location	Sub-region location whose centroid receives the text label.	Vertebra_Corpus
show_these_subreg_poi	list[int | Location] | None	If provided, only POIs with these sub-region IDs are plotted.	None

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def plot_sag_centroids(
    axs: Axes,
    ctd: POI,
    zms: tuple[float, float, float],
    poi_labelmap: dict[int, str],
    hide_centroid_labels: bool,
    cmap: ListedColormap = cm_itk,
    curve_location: Location = Location.Vertebra_Corpus,
    show_these_subreg_poi: list[int | Location] | None = None,
) -> None:
    """Overlay centroid circles and vertebra labels on a sagittal axis.

    Points are drawn as filled circles coloured by vertebra label.  Vertebra
    names from ``poi_labelmap`` are added as bold text next to the circle when
    ``hide_centroid_labels`` is ``False``.  Additional arrows and text
    annotations stored in ``ctd.info["line_segments_sag"]`` and
    ``ctd.info["text_sag"]`` are also rendered.

    Args:
        axs: Matplotlib axes on which to draw.
        ctd: POI object containing centroid coordinates and optional rendering
            hints in its ``info`` dictionary.
        zms: Voxel spacing ``(dx, dy, dz)`` in mm used to convert voxel
            coordinates to pixel coordinates.
        poi_labelmap: Mapping from vertebra region ID to display name.
        hide_centroid_labels: When ``True`` vertebra name labels are omitted.
        cmap: Colour map used to assign per-vertebra colours.
        curve_location: Sub-region location whose centroid receives the text
            label.
        show_these_subreg_poi: If provided, only POIs with these sub-region
            IDs are plotted.
    """
    # requires v_dict = dictionary of mask labels
    ctd2 = ctd
    if show_these_subreg_poi is not None:
        ctd2 = ctd.extract_subregion(*show_these_subreg_poi)
    for k1, k2, v in ctd2.items():
        color = cmap((k1 - 1) % LABEL_MAX % cmap.N)
        backgroundcolor = get_contrasting_stroke_color(color)
        # print(k, v, (v[1] * zms[1], v[0] * zms[0]), zms)
        try:
            circle = Circle(
                (v[1] * zms[1], v[0] * zms[0]),
                3,
                edgecolor=backgroundcolor,
                facecolor=color,
                linewidth=0.5,
            )
            axs.add_patch(circle)
            if not hide_centroid_labels and k2 == curve_location.value and k1 in poi_labelmap:
                axs.text(
                    4,
                    v[0] * zms[0],
                    poi_labelmap[k1],
                    fontdict={
                        "color": color,
                        "weight": "bold",
                        "fontsize": 18,
                    },
                    # bbox=dict(boxstyle="square,pad=0.2", facecolor="gray", edgecolor="none"),
                    path_effects=[withStroke(linewidth=3.0, foreground=backgroundcolor)],
                    zorder=2,
                )
        except Exception as e:
            print(e)
    if "line_segments_sag" in ctd.info:
        for color, x, (c, d) in ctd.info["line_segments_sag"]:
            if len(x) == 2:
                v = ctd[x]
            elif len(x) == 4:
                v = (np.array(ctd[x[0], x[1]]) + np.array(ctd[x[2], x[3]])) / 2

            axs.add_patch(
                FancyArrow(
                    v[1] * zms[1],
                    v[0] * zms[0],
                    c,
                    d,
                    color=cmap(color - 1 % LABEL_MAX % cmap.N),
                )
            )
    if "text_sag" in ctd.info:
        for color, x in ctd.info["text_sag"]:
            backgroundcolor = get_contrasting_stroke_color(color)
            if not isinstance(color, int) and len(color) == 2:
                color, curve_location = color  # noqa: PLW2901
            if isinstance(x, str) or len(x) == 1:
                (text) = x
                a = zms[1] * ctd[color, curve_location][1]
                b = zms[0] * ctd[color, curve_location][0]
            elif len(x) == 3:
                (text, a, b) = x
            elif len(x) == 2:
                (text, a) = x
                b = zms[0] * ctd[color, curve_location][0]
            axs.text(
                a,
                b,
                text,
                fontdict={
                    "color": color,
                    "weight": "bold",
                    "fontsize": 18,
                },
                # bbox=dict(boxstyle="square,pad=0.2", facecolor="gray", edgecolor="none"),
                path_effects=[withStroke(linewidth=3.0, foreground=backgroundcolor)],
                zorder=2,
            )

plot_cor_centroids

plot_cor_centroids(axs: Axes, ctd: POI, zms: tuple[float, float, float], poi_labelmap: dict[int, str], hide_centroid_labels: bool, cmap: ListedColormap = cm_itk, curve_location: Location = Location.Vertebra_Corpus, show_these_subreg_poi: list[int | Location] | None = None) -> None

Overlay centroid circles and vertebra labels on a coronal axis.

Analogous to :func:plot_sag_centroids but uses the z-coordinate for the horizontal axis instead of the y-coordinate. Additional arrows and text stored in ctd.info["line_segments_cor"] and ctd.info["text_cor"] are also rendered.

Parameters:

Name	Type	Description	Default
axs	Axes	Matplotlib axes on which to draw.	required
ctd	POI	POI object containing centroid coordinates and optional rendering hints in its info dictionary.	required
zms	tuple[float, float, float]	Voxel spacing (dx, dy, dz) in mm.	required
poi_labelmap	dict[int, str]	Mapping from vertebra region ID to display name.	required
hide_centroid_labels	bool	When True vertebra name labels are omitted.	required
cmap	ListedColormap	Colour map used to assign per-vertebra colours.	cm_itk
curve_location	Location	Sub-region location whose centroid receives the text label.	Vertebra_Corpus
show_these_subreg_poi	list[int | Location] | None	If provided, only POIs with these sub-region IDs are plotted.	None

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def plot_cor_centroids(
    axs: Axes,
    ctd: POI,
    zms: tuple[float, float, float],
    poi_labelmap: dict[int, str],
    hide_centroid_labels: bool,
    cmap: ListedColormap = cm_itk,
    curve_location: Location = Location.Vertebra_Corpus,
    show_these_subreg_poi: list[int | Location] | None = None,
) -> None:
    """Overlay centroid circles and vertebra labels on a coronal axis.

    Analogous to :func:`plot_sag_centroids` but uses the z-coordinate for the
    horizontal axis instead of the y-coordinate.  Additional arrows and text
    stored in ``ctd.info["line_segments_cor"]`` and ``ctd.info["text_cor"]``
    are also rendered.

    Args:
        axs: Matplotlib axes on which to draw.
        ctd: POI object containing centroid coordinates and optional rendering
            hints in its ``info`` dictionary.
        zms: Voxel spacing ``(dx, dy, dz)`` in mm.
        poi_labelmap: Mapping from vertebra region ID to display name.
        hide_centroid_labels: When ``True`` vertebra name labels are omitted.
        cmap: Colour map used to assign per-vertebra colours.
        curve_location: Sub-region location whose centroid receives the text
            label.
        show_these_subreg_poi: If provided, only POIs with these sub-region
            IDs are plotted.
    """
    ctd2 = ctd
    if show_these_subreg_poi is not None:
        ctd2 = ctd.extract_subregion(*show_these_subreg_poi)

    # requires v_dict = dictionary of mask labels
    for k1, k2, v in ctd2.items():
        color = cmap((k1 - 1) % LABEL_MAX % cmap.N)
        backgroundcolor = get_contrasting_stroke_color(color)
        try:
            circle = Circle(
                (v[2] * zms[2], v[0] * zms[0]),
                3,
                edgecolor=backgroundcolor,
                facecolor=color,
                linewidth=0.5,
            )
            axs.add_patch(circle)
            if not hide_centroid_labels and k2 == curve_location.value and k1 in poi_labelmap:
                axs.text(
                    4,
                    v[0] * zms[0],
                    poi_labelmap[k1],
                    fontdict={
                        "color": color,
                        "weight": "bold",
                        "fontsize": 18,
                    },
                    # bbox=dict(boxstyle="square,pad=0.2", facecolor="gray", edgecolor="none"),
                    path_effects=[withStroke(linewidth=3.0, foreground=backgroundcolor)],
                    zorder=2,
                )
        except Exception:
            pass

    if "line_segments_cor" in ctd.info:
        for color, x, (c, d) in ctd.info["line_segments_cor"]:
            # if isinstance(color, int):
            #    color = list(get_color_by_label(color).rgb / 255.0)

            if len(x) == 2:
                v = ctd[x]
            elif len(x) == 4:
                v = (np.array(ctd[x[0], x[1]]) + np.array(ctd[x[2], x[3]])) / 2
            axs.add_patch(
                FancyArrow(
                    v[2] * zms[2],
                    v[0] * zms[0],
                    c,
                    d,
                    color=cmap(color - 1 % LABEL_MAX % cmap.N),
                )
            )
    if "text_cor" in ctd.info:
        for color, x in ctd.info["text_cor"]:
            if isinstance(color, int):
                color = list(get_color_by_label(color).rgb / 255.0)  # noqa: PLW2901
            backgroundcolor = get_contrasting_stroke_color(color)
            if isinstance(color, Sequence) and len(color) == 2:
                color, curve_location = color  # noqa: PLW2901
            if isinstance(x, str) or len(x) == 1:
                (text) = x
                a = zms[2] * ctd[color, curve_location][2]
                b = zms[0] * ctd[color, curve_location][0]
            elif len(x) == 3:
                (text, a, b) = x
            elif len(x) == 2:
                (text, a) = x
                b = zms[0] * ctd[color, curve_location][0]
            axs.text(
                a,
                b,
                text,
                fontdict={
                    "color": color,
                    "weight": "bold",
                    "fontsize": 18,
                },
                # bbox=dict(boxstyle="square,pad=0.2", facecolor="gray", edgecolor="none"),
                path_effects=[withStroke(linewidth=3.0, foreground=backgroundcolor)],
                zorder=2,
            )

make_2d_slice

make_2d_slice(img: Image_Reference, ctd: POI, zms: tuple[float, float, float], msk: bool, visualization_type: Visualization_Type, ctd_fallback: POI, cor_savgol_filter: bool = False, to_ax: tuple[str, str, str] = ('I', 'P', 'L'), curve_location: Location = Location.Vertebra_Corpus, rescale_to_iso: bool = True, axial_heights: list[float] | None = None) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Generate 2-D sagittal, coronal, and axial projection images from a volume.

The image is reoriented to to_ax and a spline through the vertebral centroids is computed. Depending on visualization_type slices, maximum-intensity projections, depth-coloured MIPs, or mean-intensity projections are extracted. The result is optionally resampled to isotropic pixel spacing.

Parameters:

Name	Type	Description	Default
img	Image_Reference	Source image volume.	required
ctd	POI	POI object used to compute the spinal curve.	required
zms	tuple[float, float, float]	Voxel spacing (dx, dy, dz) in mm of the reoriented image.	required
msk	bool	When True the volume is treated as a label map: zero voxels are set to NaN and nearest-neighbour resampling is applied.	required
visualization_type	Visualization_Type	One of the :class:Visualization_Type enum values controlling which projection method is used.	required
ctd_fallback	POI	Fallback POI used when ctd has too few points for spline interpolation.	required
cor_savgol_filter	bool	Whether to apply a Savitzky-Golay filter to the coronal spline in :func:sag_cor_curve_projection.	False
to_ax	tuple[str, str, str]	Target orientation axes as a 3-tuple of direction codes.	('I', 'P', 'L')
curve_location	Location	Anatomical sub-region used for spline fitting.	Vertebra_Corpus
rescale_to_iso	bool	When True all returned planes are resampled to isotropic pixel spacing using the zoom factors from zms.	True
axial_heights	list[float] | None	Optional explicit axial slice heights; passed to :func:curve_projection_axial_fallback.	None

Returns:

Type	Description
ndarray	A tuple (sag, cor, axl) of 2-D float arrays. For
ndarray	attr:Visualization_Type.Maximum_Intensity_Colored_Depth without
ndarray	msk the sagittal and coronal arrays are (H, W, 3) colour
tuple[ndarray, ndarray, ndarray]	images.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def make_2d_slice(
    img: Image_Reference,
    ctd: POI,
    zms: tuple[float, float, float],
    msk: bool,
    visualization_type: Visualization_Type,
    ctd_fallback: POI,
    cor_savgol_filter: bool = False,
    to_ax: tuple[str, str, str] = ("I", "P", "L"),
    curve_location: Location = Location.Vertebra_Corpus,
    rescale_to_iso: bool = True,
    axial_heights: list[float] | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate 2-D sagittal, coronal, and axial projection images from a volume.

    The image is reoriented to ``to_ax`` and a spline through the vertebral
    centroids is computed.  Depending on ``visualization_type`` slices,
    maximum-intensity projections, depth-coloured MIPs, or mean-intensity
    projections are extracted.  The result is optionally resampled to
    isotropic pixel spacing.

    Args:
        img: Source image volume.
        ctd: POI object used to compute the spinal curve.
        zms: Voxel spacing ``(dx, dy, dz)`` in mm of the reoriented image.
        msk: When ``True`` the volume is treated as a label map: zero voxels
            are set to ``NaN`` and nearest-neighbour resampling is applied.
        visualization_type: One of the :class:`Visualization_Type` enum values
            controlling which projection method is used.
        ctd_fallback: Fallback POI used when ``ctd`` has too few points for
            spline interpolation.
        cor_savgol_filter: Whether to apply a Savitzky-Golay filter to the
            coronal spline in :func:`sag_cor_curve_projection`.
        to_ax: Target orientation axes as a 3-tuple of direction codes.
        curve_location: Anatomical sub-region used for spline fitting.
        rescale_to_iso: When ``True`` all returned planes are resampled to
            isotropic pixel spacing using the zoom factors from ``zms``.
        axial_heights: Optional explicit axial slice heights; passed to
            :func:`curve_projection_axial_fallback`.

    Returns:
        A tuple ``(sag, cor, axl)`` of 2-D float arrays.  For
        :attr:`Visualization_Type.Maximum_Intensity_Colored_Depth` without
        ``msk`` the sagittal and coronal arrays are ``(H, W, 3)`` colour
        images.
    """
    img_nii = to_nii(img)
    img_reo = img_nii.reorient_(to_ax)
    ctd_reo = ctd.reorient(img_reo.orientation)
    img_data = img_reo.get_array()

    if visualization_type in [
        visualization_type.Slice,
        visualization_type.Maximum_Intensity,
        visualization_type.Maximum_Intensity_Colored_Depth,
        visualization_type.Mean_Intensity,
    ]:
        # Make interpolated curve
        x_ctd, y_cord, z_cord = sag_cor_curve_projection(
            ctd_reo,
            img_data=img_data,
            ctd_fallback=ctd_fallback,
            cor_savgol_filter=cor_savgol_filter,
            curve_location=curve_location,
        )
        # Calculate snapshot data values depending on visualization type
        if visualization_type == Visualization_Type.Slice:
            sag, cor, axl = curve_projected_slice(
                x_ctd=x_ctd,
                img_data=img_data,
                y_cord=y_cord,
                z_cord=z_cord,
                axial_heights=axial_heights,
            )
        elif visualization_type == Visualization_Type.Maximum_Intensity:
            sag, cor, axl = curve_projected_mip(
                img_data=img_data,
                zms=zms,
                x_ctd=x_ctd,
                y_cord=y_cord,
                ctd_list=ctd_reo,
                axial_heights=axial_heights,
            )
        elif visualization_type == Visualization_Type.Maximum_Intensity_Colored_Depth:
            sag, cor, axl = curve_projected_mip(
                img_data=img_data,
                zms=zms,
                x_ctd=x_ctd,
                y_cord=y_cord,
                ctd_list=ctd_reo,
                make_colored_depth=not msk,
                axial_heights=axial_heights,
            )
        # make isotropic
        elif visualization_type == Visualization_Type.Mean_Intensity:
            sag, cor, axl = curve_projected_mean(
                img_data=img_data,
                zms=zms,
                x_ctd=x_ctd,
                y_cord=y_cord,
                ctd_list=ctd_reo,
                axial_heights=axial_heights,
            )
        else:
            raise NotImplementedError(visualization_type)

    if rescale_to_iso:
        if sag.ndim == 2:
            sag = make_isotropic2d(sag, (zms[0], zms[1]), msk=msk)
            cor = make_isotropic2d(cor, (zms[0], zms[2]), msk=msk)
            axl = make_isotropic2d(axl, (zms[1], zms[2]), msk=msk)
        elif sag.ndim == 3:  # color also encoded
            sag = make_isotropic2dpluscolor(sag, (zms[0], zms[1]), msk=msk)
            cor = make_isotropic2dpluscolor(cor, (zms[0], zms[2]), msk=msk)
            axl = make_isotropic2dpluscolor(axl, (zms[1], zms[1]), msk=msk)
        else:
            raise ValueError(f"make_2d_slice: expected sag to be ndim 2 or 3, but got shape {sag.shape}")

    sag = sag.astype(float)
    cor = cor.astype(float)
    axl = axl.astype(float)

    if msk:
        sag[sag == 0] = np.nan
        cor[cor == 0] = np.nan
        axl[axl == 0] = np.nan
    return sag, cor, axl

div0

div0(a, b, fill=0) -> np.ndarray | float

Divide a by b, replacing divisions by zero with fill.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def div0(a, b, fill=0) -> np.ndarray | float:
    """Divide ``a`` by ``b``, replacing divisions by zero with ``fill``."""
    with np.errstate(divide="ignore", invalid="ignore"):
        c = np.true_divide(a, b)
    if np.isscalar(c):
        return c if np.isfinite(c) else fill
    else:
        c[~np.isfinite(c)] = fill
        return c

to_cdt

to_cdt(ctd_bids: POI_Reference | None) -> POI | None

Load a POI from a reference and return it, or None if it is empty.

Parameters:

Name	Type	Description	Default
ctd_bids	POI_Reference | None	A path, BIDS file, or POI object to load. None is returned as-is.	required

Returns:

Type	Description
POI | None	The loaded :class:~TPTBox.POI if it contains at least one centroid,
POI | None	otherwise None.

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def to_cdt(ctd_bids: POI_Reference | None) -> POI | None:
    """Load a POI from a reference and return it, or ``None`` if it is empty.

    Args:
        ctd_bids: A path, BIDS file, or POI object to load.  ``None`` is
            returned as-is.

    Returns:
        The loaded :class:`~TPTBox.POI` if it contains at least one centroid,
        otherwise ``None``.
    """
    if ctd_bids is None:
        return None
    ctd = POI.load(ctd_bids, allow_global=True)
    if len(ctd) > 0:  # handle case if empty centroid file given
        return ctd
    print("[!][snp] To few centroids", ctd)
    return None

create_snapshot

create_snapshot(snp_path: str | Path | list[str | Path] | list[str] | list[Path], frames: list[Snapshot_Frame], crop=False, check=False, to_ax=('I', 'P', 'L'), dpi=96, verbose: bool = False) -> None

Create virtual dx, sagittal, and coronal curved-planar CT snapshots with mask overlay.

Parameters:

Name	Type	Description	Default
snp_path	str	Path to the new jpg	required
frames	List[Snapshot_Frame]	List of Images	required
crop	bool	crop output to vertebral masks (seg-vert). Defaults to False.	False
check	bool	if true, check if snap is present and do not re-create. Defaults to False.	False
to_ax	Orientation	Sets the Orientation. Can be used for flipping the image or fixing false rotations of the original inputs.	('I', 'P', 'L')
dpi	int	Set the resolution.	96

Source code in TPTBox/spine/snapshot2D/snapshot_modular.py

def create_snapshot(  # noqa: C901
    snp_path: str | Path | list[str | Path] | list[str] | list[Path],
    frames: list[Snapshot_Frame],
    crop=False,
    check=False,
    to_ax=("I", "P", "L"),
    dpi=96,
    verbose: bool = False,
) -> None:
    """Create virtual dx, sagittal, and coronal curved-planar CT snapshots with mask overlay.

    Args:
        snp_path (str): Path to the new jpg
        frames (List[Snapshot_Frame]): List of Images
        crop (bool): crop output to vertebral masks (seg-vert). Defaults to False.
        check (bool): if true, check if snap is present and do not re-create. Defaults to False.
        to_ax (Orientation): Sets the Orientation. Can be used for flipping the image or fixing false rotations of the original inputs.
        dpi (int): Set the resolution.
    """
    # Checks if snaps already exists, does nothing if true and check is true
    exist = all(Path(i).is_file() for i in snp_path) if isinstance(snp_path, list) else Path(snp_path).is_file()
    if check and exist:
        return None

    img_list = []
    frame_list = []
    frames = [f for f in frames if f is not None]
    for frame in frames:
        # PRE-PROCESSING
        img = to_nii(frame.image)
        assert img is not None
        seg = to_nii_optional(frame.segmentation, seg=True)  # can be None
        ctd = copy.deepcopy(to_cdt(frame.centroids))
        if (crop or frame.crop_msk) and seg is not None:  # crop to segmentation
            try:
                ex_slice = seg.compute_crop()
            except ValueError:
                ex_slice = None
            img = img.copy().apply_crop_(ex_slice)
            seg = seg.copy().apply_crop_(ex_slice)
            ctd = ctd.apply_crop(ex_slice).filter_points_inside_shape() if ctd is not None else None
        if frame.crop_img:  # crops image
            ex_slice = img.compute_crop(dist=0)
            img = img.apply_crop_(ex_slice)
            seg = seg.apply_crop_(ex_slice) if seg is not None else None
            ctd = ctd.apply_crop(ex_slice).filter_points_inside_shape() if ctd is not None else None
        img = img.reorient_(to_ax)
        if seg is not None:
            seg = seg.reorient_(to_ax)
            if not seg.assert_affine(img, raise_error=False):
                seg.resample_from_to_(img)
        assert not isinstance(seg, tuple), "seg is a tuple"
        if ctd is not None:
            if ctd.is_global:
                ctd = ctd.resample_from_to(img)
            if ctd.shape is None:
                POI.load(ctd, img)
            ctd = ctd.reorient(img.orientation) if ctd.shape is not None else ctd.reorient_centroids_to(img)
            if ctd.zoom not in (img.zoom, None):
                ctd.rescale_(img.zoom)

        # Preprocessing img data
        img_data = img.get_array()
        if frame.only_mask_area:
            assert seg is not None, f"only_mask_area is set but segmentation is None, got {frame}"
            seg_mask = seg.get_seg_array()
            seg_mask[seg_mask != 0] = 1
            img_data = img_data * seg_mask

        if len(img_data.shape) == 4:
            img_data = img_data[:, :, :, 0]
        if frame.gauss_filter:
            img_data = ndimage.median_filter(img_data, size=3)
        if frame.image_threshold is not None:
            img_data[img_data < frame.image_threshold] = 0  # type: ignore
        if frame.denoise_threshold is not None:
            import torch
            from torch.nn.functional import avg_pool3d

            try:
                t_arr = torch.from_numpy(img_data.copy()).unsqueeze(0).to(torch.float32)
            except Exception:
                # can't handel uint16
                t_arr = torch.from_numpy(img_data.astype(np.int32).copy()).unsqueeze(0).to(torch.float32)

            img_data_smoothed = avg_pool3d(t_arr, kernel_size=(3, 3, 3), padding=1, stride=1).squeeze(0).numpy().astype(np.int32)
            img_data[img_data_smoothed <= frame.denoise_threshold] = 0
        img = img.set_array(img_data, verbose=False)
        # PRE-PROCESSING Done

        zms = img.zoom
        try:
            ctd_fallback = POI(
                centroids={
                    (1, 50): (0, 0, img.shape[-1] // 2),
                    (2, 50): [img.shape[0] - 1, img.shape[1] - 1, img.shape[2] // 2],
                },
                orientation=img.orientation,
            )
            if ctd is None and seg is None:
                ctd_tmp = ctd_fallback
            elif frame.ignore_seg_for_centering:
                ctd_tmp = ctd_fallback
                if frame.force_show_cdt:
                    ctd = calc_centroids(seg)
            elif frame.ignore_cdt_for_centering:
                assert seg is not None, f"ignore_cdt_for_centering requires segmentation, but got None, {frame}"
                ctd_tmp = calc_centroids(seg)  # TODO BUFFER
            elif ctd is None:
                ctd_tmp = calc_centroids(seg)  # TODO BUFFER
                if frame.force_show_cdt:
                    ctd = ctd_tmp
            else:
                ctd_tmp = ctd
        except Exception:
            print("did not manage to calc ctd_tmp\n", frame)
            raise
        if frame.curve_location is None:
            if frame.show_these_subreg_poi is not None:
                l = frame.show_these_subreg_poi[0]
                frame.curve_location = Location(l) if isinstance(l, int) else l
            else:
                frame.curve_location = Location.Vertebra_Corpus
        try:
            sag_img, cor_img, axl_img = make_2d_slice(
                img,
                ctd_tmp,
                zms,
                msk=False,
                visualization_type=frame.visualization_type,
                ctd_fallback=ctd_fallback,
                to_ax=to_ax,
                cor_savgol_filter=frame.cor_savgol_filter,
                curve_location=frame.curve_location,
                rescale_to_iso=frame.rescale_to_iso,
                axial_heights=frame.axial_heights,
            )
            if seg is not None:
                sag_seg, cor_seg, axl_seg = make_2d_slice(
                    seg,
                    ctd_tmp,
                    zms,
                    msk=True,
                    visualization_type=frame.visualization_type,
                    ctd_fallback=ctd_fallback,
                    to_ax=to_ax,
                    cor_savgol_filter=frame.cor_savgol_filter,
                    curve_location=frame.curve_location,
                    rescale_to_iso=frame.rescale_to_iso,
                    axial_heights=frame.axial_heights,
                )
            else:
                sag_seg, cor_seg, axl_seg = (None, None, None)
        except Exception:
            print(frame)
            raise
        # Conversion to 2D image done, now normalization
        try:
            max_sag = np.percentile(sag_img[sag_img != 0], 99)
        except Exception:
            max_sag = 1
        try:
            max_cor = np.percentile(cor_img[cor_img != 0], 99)
        except Exception:
            max_cor = 1
        print("max sag/cor", max_sag, max_cor) if verbose else None
        ##MRT##
        if frame.mode == "MRI":
            max_intens = max(max_sag, max_cor)  # type: ignore
            wdw = Normalize(vmin=0, vmax=max_intens, clip=True)
        ##CT##
        elif frame.mode == "CT":
            wdw = wdw_hbone
        elif frame.mode == "CTs":
            wdw = wdw_sbone
        elif frame.mode == "MINMAX":
            max_intens = max(max_sag, max_cor)  # type: ignore
            min_intens = min(cor_img.min(), sag_img.min())
            wdw = Normalize(vmin=min_intens, vmax=max_intens, clip=True)
        elif frame.mode == "None":
            wdw = None
        else:
            raise ValueError(frame.mode + "is not implemented as a Normalize mode")
        alpha = frame.alpha
        # Colormap
        cmap = frame.cmap
        if isinstance(cmap, str):
            try:
                cmap = plt.get_cmap(str(cmap))
            except Exception:
                cmap = plt.get_cmap("viridis")
        # set segmentation to none if hide_segmentation
        if frame.hide_segmentation:
            sag_seg = None
            cor_seg = None
            axl_seg = None
        # set centroid to none if hide_centroids
        if frame.hide_centroids:
            ctd = None

        if not isinstance(frame.title, list):
            frame.title = [frame.title, frame.title, frame.title]

        if frame.sagittal:
            img_list.append(sag_img)
            frame_list.append(
                (
                    sag_img,
                    sag_seg,
                    ctd,
                    wdw,
                    True,
                    alpha,
                    cmap,
                    zms,
                    frame.curve_location,
                    frame.poi_labelmap,
                    frame.hide_centroid_labels,
                    frame.title[0],
                    frame,
                )
            )
        if frame.coronal:
            img_list.append(cor_img)
            frame_list.append(
                (
                    cor_img,
                    cor_seg,
                    ctd,
                    wdw,
                    False,
                    alpha,
                    cmap,
                    zms,
                    frame.curve_location,
                    frame.poi_labelmap,
                    frame.hide_centroid_labels,
                    frame.title[1],
                    frame,
                )
            )
        if frame.axial:
            img_list.append(axl_img)
            frame_list.append(
                (
                    axl_img,
                    axl_seg,
                    None,
                    wdw,
                    True,
                    alpha,
                    cmap,
                    zms,
                    frame.curve_location,
                    frame.poi_labelmap,
                    frame.hide_centroid_labels,
                    frame.title[2],
                    frame,
                )
            )

    fig, axs = create_figure(dpi, img_list, has_title=frame.title is None)
    for ax, (
        img,
        msk,
        ctd,
        wdw,
        is_sag,
        alpha,
        cmap,
        zms,
        curve_location,
        poi_labelmap,
        hide_centroid_labels,
        title,
        frame,
    ) in zip(axs, frame_list):
        if title is not None:
            ax.set_title(title, fontdict={"fontsize": 18, "color": "black"}, loc="center")
        if img.ndim == 3:
            ax.imshow(img, norm=wdw)  # type: ignore
        else:
            ax.imshow(img, cmap=plt.cm.gray, norm=wdw)  # type: ignore
        if msk is not None:
            ax.imshow(msk, cmap=cmap, alpha=alpha, vmin=1, vmax=cmap.N)
        frame: Snapshot_Frame
        if ctd is not None:
            if is_sag:
                plot_sag_centroids(
                    ax,
                    ctd,
                    zms,
                    poi_labelmap=poi_labelmap,
                    hide_centroid_labels=hide_centroid_labels,
                    cmap=cmap,
                    curve_location=curve_location,
                    show_these_subreg_poi=frame.show_these_subreg_poi,
                )
            else:
                plot_cor_centroids(
                    ax,
                    ctd,
                    zms,
                    poi_labelmap=poi_labelmap,
                    hide_centroid_labels=hide_centroid_labels,
                    cmap=cmap,
                    curve_location=curve_location,
                    show_these_subreg_poi=frame.show_these_subreg_poi,
                )

    if not isinstance(snp_path, list):
        snp_path = [str(snp_path)]
    for path in snp_path:
        fig.savefig(str(path))
        print("[*] Snapshot saved:", path) if verbose else None
    plt.close()
    return snp_path

Snapshot Templates

TPTBox.spine.snapshot2D.snapshot_templates

mip_shot

mip_shot(ct_ref: BIDS_FILE, vert_msk: Image_Reference, subreg_ctd: POI_Reference, subreg_msk: Image_Reference | None = None, out_path: str | Path | None = None) -> None

Save a multi-panel CT snapshot combining slice, mean-intensity, and coloured-depth MIP views.

Three frames are always included:

1. Standard CT slice view with vertebra segmentation and centroids (sagittal + coronal).
2. Mean-intensity projection of the masked vertebra region only (sagittal + coronal, no segmentation overlay).
3. Coloured-depth maximum-intensity projection of the bone above a 100 HU threshold (coronal only).

An optional fourth frame with the subregion mask is appended when subreg_msk is provided.

Parameters:

Name	Type	Description	Default
ct_ref	BIDS_FILE	BIDS file pointing to the CT image.	required
vert_msk	Image_Reference	Vertebra segmentation mask.	required
subreg_ctd	POI_Reference	Subregion centroids (POI reference).	required
subreg_msk	Image_Reference | None	Optional subregion segmentation mask for an extra frame.	None
out_path	str | Path | None	Output file path. When None the path is derived automatically from ct_ref with desc=mip.	None

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def mip_shot(
    ct_ref: BIDS_FILE,
    vert_msk: Image_Reference,
    subreg_ctd: POI_Reference,
    subreg_msk: Image_Reference | None = None,
    out_path: str | Path | None = None,
) -> None:
    """Save a multi-panel CT snapshot combining slice, mean-intensity, and coloured-depth MIP views.

    Three frames are always included:

    1. Standard CT slice view with vertebra segmentation and centroids
       (sagittal + coronal).
    2. Mean-intensity projection of the masked vertebra region only
       (sagittal + coronal, no segmentation overlay).
    3. Coloured-depth maximum-intensity projection of the bone above a 100 HU
       threshold (coronal only).

    An optional fourth frame with the subregion mask is appended when
    ``subreg_msk`` is provided.

    Args:
        ct_ref: BIDS file pointing to the CT image.
        vert_msk: Vertebra segmentation mask.
        subreg_ctd: Subregion centroids (POI reference).
        subreg_msk: Optional subregion segmentation mask for an extra frame.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``ct_ref`` with ``desc=mip``.
    """
    frames = [
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode="CT",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=True,
            cor_savgol_filter=False,
            # cmap="viridis",#ListedColormap(cmap),
        ),
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode="MINMAX",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=True,
            hide_segmentation=True,
            visualization_type=Visualization_Type.Mean_Intensity,
            only_mask_area=True,
        ),
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode="None",
            sagittal=False,
            coronal=True,
            axial=False,
            crop_msk=True,
            hide_segmentation=True,
            visualization_type=Visualization_Type.Maximum_Intensity_Colored_Depth,
            image_threshold=100,
            denoise_threshold=150,
            only_mask_area=False,
        ),
    ]

    if subreg_msk is not None:
        frames.append(
            Snapshot_Frame(
                image=ct_ref,
                segmentation=subreg_msk,
                centroids=subreg_ctd,
                mode="CT",
                sagittal=True,
                coronal=True,
                axial=False,
                crop_msk=True,
                cor_savgol_filter=False,
                # cmap="viridis",#ListedColormap(cmap),
            )
        )

    if out_path is None:
        out_path = ct_ref.get_changed_path(file_type="png", bids_format="snp", info={"desc": "mip"})
    create_snapshot(snp_path=[out_path], frames=frames)

sacrum_shot

sacrum_shot(ct_ref: BIDS_FILE, vert_msk: Image_Reference, subreg_ctd: POI_Reference, add_ctd: list[POI_Reference] | None = None, mip_c: bool = False, out_path: str | Path | None = None) -> None

Save a CT snapshot focused on the sacrum region.

A primary sagittal + coronal frame is always created. Additional coronal frames are appended for each extra centroid file in add_ctd. When mip_c is True a coloured-depth MIP coronal frame is added.

Parameters:

Name	Type	Description	Default
ct_ref	BIDS_FILE	BIDS file pointing to the CT image. If an instance of :class:~TPTBox.BIDS_FILE, the image is pre-loaded and reoriented.	required
vert_msk	Image_Reference	Vertebra segmentation mask.	required
subreg_ctd	POI_Reference	Primary subregion centroids (POI reference).	required
add_ctd	list[POI_Reference] | None	Optional list of additional centroid references, each rendered as an extra coronal frame.	None
mip_c	bool	When True a coloured-depth maximum-intensity projection frame is appended.	False
out_path	str | Path | None	Output file path. When None the path is derived automatically from ct_ref with desc=sacrum.	None

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def sacrum_shot(
    ct_ref: BIDS_FILE,
    vert_msk: Image_Reference,
    subreg_ctd: POI_Reference,
    add_ctd: list[POI_Reference] | None = None,
    mip_c: bool = False,
    out_path: str | Path | None = None,
) -> None:
    """Save a CT snapshot focused on the sacrum region.

    A primary sagittal + coronal frame is always created.  Additional coronal
    frames are appended for each extra centroid file in ``add_ctd``.  When
    ``mip_c`` is ``True`` a coloured-depth MIP coronal frame is added.

    Args:
        ct_ref: BIDS file pointing to the CT image.  If an instance of
            :class:`~TPTBox.BIDS_FILE`, the image is pre-loaded and
            reoriented.
        vert_msk: Vertebra segmentation mask.
        subreg_ctd: Primary subregion centroids (POI reference).
        add_ctd: Optional list of additional centroid references, each
            rendered as an extra coronal frame.
        mip_c: When ``True`` a coloured-depth maximum-intensity projection
            frame is appended.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``ct_ref`` with ``desc=sacrum``.
    """
    if isinstance(ct_ref, BIDS_FILE):
        ct_ref = ct_ref.open_nii().reorient_()
    frames = [
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode="CT",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=True,
            cor_savgol_filter=False,
            # cmap="viridis",#ListedColormap(cmap),
        )
    ]
    if add_ctd is not None:
        for c in add_ctd:
            frames.append(  # noqa: PERF401
                Snapshot_Frame(
                    image=ct_ref,
                    segmentation=vert_msk,
                    centroids=c,
                    mode="CT",
                    sagittal=False,
                    coronal=True,
                    axial=False,
                    crop_msk=True,
                    cor_savgol_filter=False,
                    # cmap="viridis",#ListedColormap(cmap),
                )
            )

    if mip_c:
        frames.append(
            Snapshot_Frame(
                image=ct_ref,
                segmentation=vert_msk,
                centroids=subreg_ctd,
                mode="None",
                sagittal=False,
                coronal=True,
                axial=False,
                crop_msk=True,
                hide_segmentation=True,
                visualization_type=Visualization_Type.Maximum_Intensity_Colored_Depth,
                image_threshold=100,
                denoise_threshold=150,
                only_mask_area=False,
            )
        )
    if out_path is None:
        out_path = ct_ref.get_changed_path(file_type="png", bids_format="snp", info={"desc": "sacrum"})
    create_snapshot(snp_path=[out_path], frames=frames)

spline_shot

spline_shot(ct_ref: BIDS_FILE, vert_msk: Image_Reference, subreg_ctd: POI_Reference, spline_nii: NII, add_info: str = '', out_path: str | Path | None = None) -> None

Save a two-frame CT snapshot showing the spinal spline interpolation result.

Frame 1 shows the original vertebra mask overlaid on the CT (sagittal only). Frame 2 shows the spline NIfTI segmentation overlaid on the CT (sagittal + coronal) so that the spline quality can be assessed visually.

Parameters:

Name	Type	Description	Default
ct_ref	BIDS_FILE	BIDS file pointing to the CT image.	required
vert_msk	Image_Reference	Vertebra segmentation mask.	required
subreg_ctd	POI_Reference	Subregion centroids (POI reference).	required
spline_nii	NII	NIfTI image containing the fitted spline segmentation.	required
add_info	str	Optional string appended to the desc field of the auto-generated output path.	''
out_path	str | Path | None	Output file path. When None the path is derived automatically from ct_ref.	None

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def spline_shot(
    ct_ref: BIDS_FILE,
    vert_msk: Image_Reference,
    subreg_ctd: POI_Reference,
    spline_nii: NII,
    add_info: str = "",
    out_path: str | Path | None = None,
) -> None:
    """Save a two-frame CT snapshot showing the spinal spline interpolation result.

    Frame 1 shows the original vertebra mask overlaid on the CT (sagittal
    only).  Frame 2 shows the spline NIfTI segmentation overlaid on the CT
    (sagittal + coronal) so that the spline quality can be assessed visually.

    Args:
        ct_ref: BIDS file pointing to the CT image.
        vert_msk: Vertebra segmentation mask.
        subreg_ctd: Subregion centroids (POI reference).
        spline_nii: NIfTI image containing the fitted spline segmentation.
        add_info: Optional string appended to the ``desc`` field of the
            auto-generated output path.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``ct_ref``.
    """
    frames = [
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode="CT",
            sagittal=True,
            axial=False,
            crop_msk=False,
        ),
        Snapshot_Frame(
            image=ct_ref,
            segmentation=spline_nii,
            centroids=subreg_ctd,
            mode="CT",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=False,
        ),
    ]
    if out_path is None:
        out_path = ct_ref.get_changed_path(
            file_type="png",
            bids_format="snp",
            info={"desc": f"spline-interpolation-{add_info}"},
        )
    create_snapshot(snp_path=[out_path], frames=frames)

snapshot

snapshot(ref: Image_Reference, vert_msk: Image_Reference, subreg_ctd: POI_Reference, subreg_msk: Image_Reference | None = None, out_path: str | Path | list[str | Path] | list[Path] | None = None, mode: str = 'MINMAX', crop: bool = False) -> list[str | Path]

Save a standard vertebra snapshot, auto-detecting CT vs. MRI mode.

Delegates to :func:mri_snapshot after setting mode to "CT" for BIDS CT files or "MRI" for all other BIDS files.

Parameters:

Name	Type	Description	Default
ref	Image_Reference	Source image reference. A :class:~TPTBox.BIDS_FILE whose bids_format is "ct" is treated as CT; everything else is treated as MRI.	required
vert_msk	Image_Reference	Vertebra segmentation mask.	required
subreg_ctd	POI_Reference	Subregion centroids (POI reference).	required
subreg_msk	Image_Reference | None	Optional subregion segmentation mask for an extra frame.	None
out_path	str | Path | list[str | Path] | list[Path] | None	Output path(s). None triggers automatic path derivation.	None
mode	str	Intensity window mode passed to :class:~.Snapshot_Frame when ref is not a :class:~TPTBox.BIDS_FILE.	'MINMAX'
crop	bool	Whether to crop the output to the segmentation bounding box.	False

Returns:

Type	Description
list[str | Path]	List of saved output paths.

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def snapshot(
    ref: Image_Reference,
    vert_msk: Image_Reference,
    subreg_ctd: POI_Reference,
    subreg_msk: Image_Reference | None = None,
    out_path: str | Path | list[str | Path] | list[Path] | None = None,
    mode: str = "MINMAX",
    crop: bool = False,
) -> list[str | Path]:
    """Save a standard vertebra snapshot, auto-detecting CT vs. MRI mode.

    Delegates to :func:`mri_snapshot` after setting ``mode`` to ``"CT"`` for
    BIDS CT files or ``"MRI"`` for all other BIDS files.

    Args:
        ref: Source image reference.  A :class:`~TPTBox.BIDS_FILE` whose
            ``bids_format`` is ``"ct"`` is treated as CT; everything else
            is treated as MRI.
        vert_msk: Vertebra segmentation mask.
        subreg_ctd: Subregion centroids (POI reference).
        subreg_msk: Optional subregion segmentation mask for an extra frame.
        out_path: Output path(s).  ``None`` triggers automatic path derivation.
        mode: Intensity window mode passed to :class:`~.Snapshot_Frame` when
            ``ref`` is not a :class:`~TPTBox.BIDS_FILE`.
        crop: Whether to crop the output to the segmentation bounding box.

    Returns:
        List of saved output paths.
    """
    if isinstance(ref, BIDS_FILE):
        mode = "CT" if ref.bids_format == "ct" else "MRI"

    return mri_snapshot(ref, vert_msk, subreg_ctd, subreg_msk, out_path, mode, crop=crop)  # type: ignore

mri_snapshot

mri_snapshot(mrt_ref: BIDS_FILE, vert_msk: Image_Reference, subreg_ctd: POI_Reference, subreg_msk: Image_Reference | None = None, out_path: str | Path | list[str | Path] | list[Path] | None = None, mode: str = 'MRI', crop: bool = False) -> list[str | Path]

Save a two-frame (or three-frame) MRI vertebra snapshot.

Frame 1: image with centroids but without segmentation overlay (to check centroid placement without clutter). Frame 2: image with both centroids and segmentation overlay. Frame 3 (optional): image with the subregion mask overlay, added when subreg_msk is provided.

Parameters:

Name	Type	Description	Default
mrt_ref	BIDS_FILE	BIDS file pointing to the MRI image.	required
vert_msk	Image_Reference	Vertebra segmentation mask.	required
subreg_ctd	POI_Reference	Subregion centroids (POI reference).	required
subreg_msk	Image_Reference | None	Optional subregion segmentation mask for an extra frame.	None
out_path	str | Path | list[str | Path] | list[Path] | None	Output path(s). When None the path is derived automatically from mrt_ref with desc=vert.	None
mode	str	Intensity window mode; defaults to "MRI".	'MRI'
crop	bool	Whether to crop output frames to the segmentation bounding box.	False

Returns:

Type	Description
list[str | Path]	List of saved output paths.

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def mri_snapshot(
    mrt_ref: BIDS_FILE,
    vert_msk: Image_Reference,
    subreg_ctd: POI_Reference,
    subreg_msk: Image_Reference | None = None,
    out_path: str | Path | list[str | Path] | list[Path] | None = None,
    mode: str = "MRI",
    crop: bool = False,
) -> list[str | Path]:
    """Save a two-frame (or three-frame) MRI vertebra snapshot.

    Frame 1: image with centroids but without segmentation overlay (to
    check centroid placement without clutter).
    Frame 2: image with both centroids and segmentation overlay.
    Frame 3 (optional): image with the subregion mask overlay, added when
    ``subreg_msk`` is provided.

    Args:
        mrt_ref: BIDS file pointing to the MRI image.
        vert_msk: Vertebra segmentation mask.
        subreg_ctd: Subregion centroids (POI reference).
        subreg_msk: Optional subregion segmentation mask for an extra frame.
        out_path: Output path(s).  When ``None`` the path is derived
            automatically from ``mrt_ref`` with ``desc=vert``.
        mode: Intensity window mode; defaults to ``"MRI"``.
        crop: Whether to crop output frames to the segmentation bounding box.

    Returns:
        List of saved output paths.
    """
    frames = [
        Snapshot_Frame(
            image=mrt_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode=mode,
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=crop,
            hide_segmentation=True,
        ),
        Snapshot_Frame(
            image=mrt_ref,
            segmentation=vert_msk,
            centroids=subreg_ctd,
            mode=mode,
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=crop,
        ),
    ]
    if subreg_msk is not None:
        frames.append(
            Snapshot_Frame(
                image=mrt_ref,
                segmentation=subreg_msk,
                centroids=subreg_ctd,
                mode=mode,
                sagittal=True,
                coronal=True,
                axial=False,
                crop_msk=crop,
            )
        )
    if out_path is None:
        out_path = mrt_ref.get_changed_path(file_type="png", bids_format="snp", info={"desc": "vert"})
    if not isinstance(out_path, list):
        out_path = [out_path]
    create_snapshot(snp_path=out_path, frames=frames)
    return out_path

vibe_snapshot

vibe_snapshot(mrt_ref: tuple[BIDS_FILE, BIDS_FILE, BIDS_FILE, BIDS_FILE], vert_msk: Image_Reference | None, subreg_ctd: POI_Reference, subreg_msk: Image_Reference, hide_centroids: bool = False, out_path: str | Path | None = None, verbose: bool = False) -> str | Path

Save a VIBE multi-contrast MRI snapshot with four Dixon phases.

Creates one sagittal frame for each of the four VIBE contrast phases in mrt_ref, all with the subregion mask overlay. An optional fifth frame shows the first phase in sagittal + coronal with the vertebra mask.

Parameters:

Name	Type	Description	Default
mrt_ref	tuple[BIDS_FILE, BIDS_FILE, BIDS_FILE, BIDS_FILE]	Four-tuple of BIDS files corresponding to the VIBE Dixon phases (e.g. in-phase, opposed-phase, water, fat).	required
vert_msk	Image_Reference | None	Optional vertebra segmentation mask used for an extra sagittal + coronal frame. Skipped when None.	required
subreg_ctd	POI_Reference	Subregion centroids (POI reference).	required
subreg_msk	Image_Reference	Subregion segmentation mask overlaid on all four phases.	required
hide_centroids	bool	When True centroid markers are hidden.	False
out_path	str | Path | None	Output file path. When None the path is derived automatically from mrt_ref[0] with desc=vibe.	None
verbose	bool	When True the saved path is printed to stdout.	False

Returns:

Type	Description
str | Path	Path to the saved snapshot image.

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def vibe_snapshot(
    mrt_ref: tuple[BIDS_FILE, BIDS_FILE, BIDS_FILE, BIDS_FILE],
    vert_msk: Image_Reference | None,
    subreg_ctd: POI_Reference,
    subreg_msk: Image_Reference,
    hide_centroids: bool = False,
    out_path: str | Path | None = None,
    verbose: bool = False,
) -> str | Path:
    """Save a VIBE multi-contrast MRI snapshot with four Dixon phases.

    Creates one sagittal frame for each of the four VIBE contrast phases in
    ``mrt_ref``, all with the subregion mask overlay.  An optional fifth
    frame shows the first phase in sagittal + coronal with the vertebra mask.

    Args:
        mrt_ref: Four-tuple of BIDS files corresponding to the VIBE Dixon
            phases (e.g. in-phase, opposed-phase, water, fat).
        vert_msk: Optional vertebra segmentation mask used for an extra
            sagittal + coronal frame.  Skipped when ``None``.
        subreg_ctd: Subregion centroids (POI reference).
        subreg_msk: Subregion segmentation mask overlaid on all four phases.
        hide_centroids: When ``True`` centroid markers are hidden.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``mrt_ref[0]`` with ``desc=vibe``.
        verbose: When ``True`` the saved path is printed to stdout.

    Returns:
        Path to the saved snapshot image.
    """
    # print(type(mrt_ref[0]))
    frames = [
        Snapshot_Frame(
            image=mrt_ref[0],
            segmentation=subreg_msk,
            centroids=subreg_ctd,
            mode="MRI",
            sagittal=True,
            coronal=False,
            axial=False,
            crop_msk=False,
            hide_centroids=hide_centroids,
        ),
        Snapshot_Frame(
            image=mrt_ref[1],
            segmentation=subreg_msk,
            centroids=subreg_ctd,
            mode="MRI",
            sagittal=True,
            coronal=False,
            axial=False,
            crop_msk=False,
            hide_centroids=hide_centroids,
        ),
        Snapshot_Frame(
            image=mrt_ref[2],
            segmentation=subreg_msk,
            centroids=subreg_ctd,
            mode="MRI",
            sagittal=True,
            coronal=False,
            axial=False,
            crop_msk=False,
            hide_centroids=hide_centroids,
        ),
        Snapshot_Frame(
            image=mrt_ref[3],
            segmentation=subreg_msk,
            centroids=subreg_ctd,
            mode="MRI",
            sagittal=True,
            coronal=False,
            axial=False,
            crop_msk=False,
            hide_centroids=hide_centroids,
        ),
    ]
    if vert_msk is not None:
        frames.append(
            Snapshot_Frame(
                image=mrt_ref[0],
                segmentation=vert_msk,
                centroids=subreg_ctd,
                mode="MRI",
                sagittal=True,
                coronal=True,
                axial=False,
                crop_msk=False,
                hide_centroids=hide_centroids,
            )
        )
    if out_path is None:
        out_path = mrt_ref[0].get_changed_path(file_type="png", bids_format="snp", info={"desc": "vibe"})
    create_snapshot(snp_path=[out_path], frames=frames, verbose=verbose)
    return out_path

ct_mri_snapshot

ct_mri_snapshot(mrt_ref: Image_Reference, ct_ref: Image_Reference, vert_msk_mrt: Image_Reference | None = None, subreg_ctd_mrt: POI_Reference | None = None, vert_msk_ct: Image_Reference | None = None, subreg_ctd_ct: POI_Reference | None = None, out_path: str | Path | None = None, return_frames: bool = False) -> list | str | Path

Save a two-frame snapshot comparing an MRI and a CT side by side.

Frame 1: MRI in "MRI" mode with optional vertebra mask and centroids. Frame 2: CT in "CT" mode with optional vertebra mask and centroids.

Parameters:

Name	Type	Description	Default
mrt_ref	Image_Reference	MRI image reference.	required
ct_ref	Image_Reference	CT image reference.	required
vert_msk_mrt	Image_Reference | None	Optional vertebra segmentation for the MRI frame.	None
subreg_ctd_mrt	POI_Reference | None	Optional subregion centroids for the MRI frame.	None
vert_msk_ct	Image_Reference | None	Optional vertebra segmentation for the CT frame.	None
subreg_ctd_ct	POI_Reference | None	Optional subregion centroids for the CT frame.	None
out_path	str | Path | None	Output file path. When None the path is derived automatically from mrt_ref (must be a :class:~TPTBox.BIDS_FILE).	None
return_frames	bool	When True the list of :class:~.Snapshot_Frame objects is returned without saving.	False

Returns:

Type	Description
list | str | Path	The list of :class:~.Snapshot_Frame objects when return_frames
list | str | Path	is True, otherwise the path to the saved snapshot image.

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def ct_mri_snapshot(
    mrt_ref: Image_Reference,
    ct_ref: Image_Reference,
    vert_msk_mrt: Image_Reference | None = None,
    subreg_ctd_mrt: POI_Reference | None = None,
    vert_msk_ct: Image_Reference | None = None,
    subreg_ctd_ct: POI_Reference | None = None,
    out_path: str | Path | None = None,
    return_frames: bool = False,
) -> list | str | Path:
    """Save a two-frame snapshot comparing an MRI and a CT side by side.

    Frame 1: MRI in ``"MRI"`` mode with optional vertebra mask and centroids.
    Frame 2: CT in ``"CT"`` mode with optional vertebra mask and centroids.

    Args:
        mrt_ref: MRI image reference.
        ct_ref: CT image reference.
        vert_msk_mrt: Optional vertebra segmentation for the MRI frame.
        subreg_ctd_mrt: Optional subregion centroids for the MRI frame.
        vert_msk_ct: Optional vertebra segmentation for the CT frame.
        subreg_ctd_ct: Optional subregion centroids for the CT frame.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``mrt_ref`` (must be a :class:`~TPTBox.BIDS_FILE`).
        return_frames: When ``True`` the list of
            :class:`~.Snapshot_Frame` objects is returned without saving.

    Returns:
        The list of :class:`~.Snapshot_Frame` objects when ``return_frames``
        is ``True``, otherwise the path to the saved snapshot image.
    """
    frames = [
        Snapshot_Frame(
            image=mrt_ref,
            segmentation=vert_msk_mrt,
            centroids=subreg_ctd_mrt,
            mode="MRI",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=False,
        ),
        Snapshot_Frame(
            image=ct_ref,
            segmentation=vert_msk_ct,
            centroids=subreg_ctd_ct,
            mode="CT",
            sagittal=True,
            coronal=True,
            axial=False,
            crop_msk=False,
        ),
    ]
    if return_frames:
        return frames
    if out_path is None:
        assert isinstance(mrt_ref, BIDS_FILE)
        out_path = mrt_ref.get_changed_path(file_type="png", bids_format="snp", info={"desc": "vert-ct-mri"})
    create_snapshot(snp_path=[out_path], frames=frames)
    return out_path

poi_snapshot

poi_snapshot(ct_nii: BIDS_FILE, vert_msk: Image_Reference | None, subreg_ctd: POI_Reference, out_path: str | Path | None = None) -> None

Save a multi-frame CT snapshot visualising ligament attachment-point POIs.

The POIs are partitioned into lateral groups (dorsal _D, sinister _S, median) and ligament groups (ALL, PLL, FL). Seven frames are produced covering different sagittal and coronal views of these groups.

Parameters:

Name	Type	Description	Default
ct_nii	BIDS_FILE	BIDS file pointing to the CT image.	required
vert_msk	Image_Reference | None	Optional vertebra segmentation mask.	required
subreg_ctd	POI_Reference	POI reference containing the ligament attachment points.	required
out_path	str | Path | None	Output file path. When None the path is derived automatically from ct_nii with desc=poi.	None

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def poi_snapshot(
    ct_nii: BIDS_FILE,
    vert_msk: Image_Reference | None,
    subreg_ctd: POI_Reference,
    out_path: str | Path | None = None,
) -> None:
    """Save a multi-frame CT snapshot visualising ligament attachment-point POIs.

    The POIs are partitioned into lateral groups (dorsal ``_D``, sinister
    ``_S``, median) and ligament groups (ALL, PLL, FL).  Seven frames are
    produced covering different sagittal and coronal views of these groups.

    Args:
        ct_nii: BIDS file pointing to the CT image.
        vert_msk: Optional vertebra segmentation mask.
        subreg_ctd: POI reference containing the ligament attachment points.
        out_path: Output file path.  When ``None`` the path is derived
            automatically from ``ct_nii`` with ``desc=poi``.
    """
    # conversion_poi = {
    #    "SSL": 81,  # this POI is not included in our POI list
    #    "ALL_CR_S": 109,
    #    "ALL_CR": 101,
    #    "ALL_CR_D": 117,
    #    "ALL_CA_S": 111,
    #    "ALL_CA": 103,
    #    "ALL_CA_D": 119,
    #    "PLL_CR_S": 110,
    #    "PLL_CR": 102,
    #    "PLL_CR_D": 118,
    #    "PLL_CA_S": 112,
    #    "PLL_CA": 104,
    #    "PLL_CA_D": 120,
    #    "FL_CR_S": 149,
    #    "FL_CR": 125,
    #    "FL_CR_D": 141,
    #    "FL_CA_S": 151,
    #    "FL_CA": 127,
    #    "FL_CA_D": 143,
    #    "ISL_CR": 134,
    #    "ISL_CA": 136,
    #    "ITL_S": 142,
    #    "ITL_D": 144,
    # }
    poi_all = POI.load(subreg_ctd)
    sinister = {}
    median = {}
    dorsal = {}
    all_ = {}
    pll = {}
    fl = {}
    # isl = {}
    # itl = {}
    from TPTBox.core.vert_constants import Location, conversion_poi2text

    for k, k2, v in poi_all.items():
        s = k2
        if conversion_poi2text[s].endswith("_D"):
            dorsal[k] = v
        elif conversion_poi2text[s].endswith("_S"):
            sinister[k] = v
        else:
            median[k] = v
        if "ALL" in conversion_poi2text[s]:
            all_[k] = v
        elif "PLL" in conversion_poi2text[s]:
            pll[k] = v
        elif "FL" in conversion_poi2text[s]:
            fl[k] = v

    frames = [
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(dorsal),
            mode="CT",
            curve_location=Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Right,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(median),
            mode="CT",
            curve_location=Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Median,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(sinister),
            mode="CT",
            curve_location=Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Left,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(all_),
            mode="CT",
            coronal=True,
            curve_location=Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Median,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(pll),
            mode="CT",
            coronal=True,
            curve_location=Location.Ligament_Attachment_Point_Posterior_Longitudinal_Superior_Median,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(fl),
            mode="MINMAX",
            sagittal=True,
            coronal=False,
            hide_segmentation=True,
            only_mask_area=True,
            visualization_type=Visualization_Type.Mean_Intensity,
            curve_location=Location.Ligament_Attachment_Point_Flava_Superior_Median,
        ),
        Snapshot_Frame(
            image=ct_nii,
            segmentation=vert_msk,
            centroids=poi_all.copy(fl),
            mode="CT",
            sagittal=False,
            coronal=True,
            curve_location=Location.Ligament_Attachment_Point_Flava_Superior_Median,
        ),
    ]
    if out_path is None:
        out_path = ct_nii.get_changed_path(file_type="png", bids_format="snp", info={"desc": "poi"})
    create_snapshot(snp_path=[out_path], frames=frames)

make_snap

make_snap(file) -> None

Generate a POI snapshot for the given file if it matches the expected subject index.

Source code in TPTBox/spine/snapshot2D/snapshot_templates.py

def make_snap(file) -> None:
    """Generate a POI snapshot for the given file if it matches the expected subject index."""
    idx = file.stem
    if idx != "63":
        return
    # print(id)
    try:
        base_path = f"/media/data/robert/datasets/dataset-poi/derivatives/WS_{idx}"

        ct_file = BIDS_FILE(
            next(Path(base_path).glob("ses-*/sub-WS_*_ses-*_seq-*_space-aligASL_ct.nii.gz")),
            "/media/data/robert/datasets/dataset-poi/",
            verbose=False,
        )
        f = next(Path(base_path).glob("ses-*/sub-WS_*_ses-*_seq-*_seg-subreg_space-aligASL_msk.nii.gz"))
        # g = next(Path(base_path).glob("ses-*/sub-WS_*_ses-*_seq-*_seg-subreg_space-aligASL_msk.nii.gz"))
        poi_snapshot(
            ct_file,
            f,
            file,
            out_path=f"/media/data/robert/datasets/dataset-poi/snapshot/{idx}.png",
        )
    except StopIteration:
        print(idx, "stopIteration")

Distances

TPTBox.spine.spinestats.distances

compute_all_distances

compute_all_distances(poi: POI, vert: Image_Reference | None = None, subreg: Image_Reference | None = None, all_pois_computed: bool = False, recompute: bool = False) -> POI

Compute all registered anatomical distances and store them in poi.info.

Iterates over :data:distances_funs and calls :func:_compute_distance for each entry. Results are accumulated in the returned POI object under the corresponding key in poi.info.

Parameters:

Name	Type	Description	Default
poi	POI	POI object to compute distances for and to store results in.	required
vert	Image_Reference | None	Vertebra segmentation image required when POIs have not been pre-computed. May be None when all_pois_computed=True.	None
subreg	Image_Reference | None	Subregion segmentation image required when POIs have not been pre-computed. May be None when all_pois_computed=True.	None
all_pois_computed	bool	When True the function skips POI computation and uses the existing centroids in poi directly.	False
recompute	bool	When True existing entries in poi.info are overwritten; otherwise already-computed keys are skipped.	False

Returns:

Type	Description
POI	The updated :class:~TPTBox.POI object with distance dictionaries
POI	stored under the keys defined in :data:distances_funs.

Source code in TPTBox/spine/spinestats/distances.py

def compute_all_distances(
    poi: POI,
    vert: Image_Reference | None = None,
    subreg: Image_Reference | None = None,
    all_pois_computed: bool = False,
    recompute: bool = False,
) -> POI:
    """Compute all registered anatomical distances and store them in ``poi.info``.

    Iterates over :data:`distances_funs` and calls :func:`_compute_distance`
    for each entry.  Results are accumulated in the returned POI object under
    the corresponding key in ``poi.info``.

    Args:
        poi: POI object to compute distances for and to store results in.
        vert: Vertebra segmentation image required when POIs have not been
            pre-computed.  May be ``None`` when ``all_pois_computed=True``.
        subreg: Subregion segmentation image required when POIs have not been
            pre-computed.  May be ``None`` when ``all_pois_computed=True``.
        all_pois_computed: When ``True`` the function skips POI computation
            and uses the existing centroids in ``poi`` directly.
        recompute: When ``True`` existing entries in ``poi.info`` are
            overwritten; otherwise already-computed keys are skipped.

    Returns:
        The updated :class:`~TPTBox.POI` object with distance dictionaries
        stored under the keys defined in :data:`distances_funs`.
    """
    for key, (l1, l2) in distances_funs.items():
        poi = _compute_distance(poi, l1, l2, key, vert, subreg, all_pois_computed, recompute)
    return poi

Angles

TPTBox.spine.spinestats.angles

MoveTo

Bases: Enum

Enum selecting which endplate or centre point of a vertebra is used as the measurement anchor.

Source code in TPTBox/spine/spinestats/angles.py

class MoveTo(Enum):
    """Enum selecting which endplate or centre point of a vertebra is used as the measurement anchor."""

    TOP = auto()
    BOTTOM = auto()
    CENTER = auto()

    def has_point(self, v: Vertebra_Instance | int, poi: POI) -> bool:
        """Check if the given vertebra has a specific POI at the current MoveTo position.

        Args:
            v (Vertebra_Instance | int): The vertebra instance or its ID.
            poi (POI): The point of interest data structure.

        Returns:
            bool: True if the POI exists for the specified MoveTo position, otherwise False.
        """
        if isinstance(v, int):
            v = Vertebra_Instance(v)
        try:
            self.get_point(v, poi)
        except KeyError:
            return False
        return True

    def get_location(self, v: Vertebra_Instance | int, poi: POI) -> tuple:
        """Determine the anatomical POI coordinates for the current MoveTo position in a vertebra.

        Args:
            v (Vertebra_Instance | int): The vertebra instance or its ID.
            poi (POI): The point of interest data structure.

        Returns:
            tuple: The location tuple that defines the position within the vertebra.
        """
        if isinstance(v, int):
            v = Vertebra_Instance(v)
        if self == self.CENTER:
            return (v, 50)
        elif self == self.BOTTOM:
            # Test IVD
            subreg = Location.Vertebra_Disc
            if (v, subreg) in poi:
                return (v, subreg)
            # Test if it has next
            subreg = Location.Additional_Vertebral_Body_Middle_Inferior_Median
            if (v, subreg) in poi:
                return (v, subreg)
            # Test if it has next POINT
            next_vert = v.get_next_poi(poi)
            if next_vert is not None and (next_vert, 50) in poi:
                return (v, subreg, next_vert, subreg)
        elif self == self.TOP:
            prev_vert = v.get_previous_poi(poi)
            # Test IVD
            subreg = Location.Vertebra_Disc
            if prev_vert is not None and (prev_vert, subreg) in poi:
                return (prev_vert, subreg)
            # Test if it has next
            subreg = Location.Dens_axis
            if (v, subreg) in poi:
                return (v, subreg)
            # Test if it has next
            subreg = Location.Additional_Vertebral_Body_Middle_Superior_Median
            if (v, subreg) in poi:
                return (v, subreg)
            # Test if it has next POINT
            if prev_vert is not None and (prev_vert, 50) in poi:
                return (v, subreg, prev_vert, subreg)
        return (v, 50)

    def get_point(self, v: Vertebra_Instance | int, poi: POI) -> np.ndarray:
        """Retrieve the 3D coordinates of a POI in a vertebra at the current MoveTo position.

        Args:
            v (Vertebra_Instance | int): The vertebra instance or its ID.
            poi (POI): The point of interest data structure.

        Returns:
            np.ndarray: The 3D coordinates of the specified POI.

        Raises:
            NotImplementedError: If the POI cannot be determined.
        """
        a = self.get_location(v, poi)
        if len(a) == 2:
            return np.array(poi[a])
        elif len(a) == 4:
            return (np.array(poi[a[0], a[1]]) + np.array(poi[a[2], a[3]])) / 2
        raise NotImplementedError(v, poi)

has_point

has_point(v: Vertebra_Instance | int, poi: POI) -> bool

Check if the given vertebra has a specific POI at the current MoveTo position.

Parameters:

Name	Type	Description	Default
v	Vertebra_Instance | int	The vertebra instance or its ID.	required
poi	POI	The point of interest data structure.	required

Returns:

Name	Type	Description
bool	bool	True if the POI exists for the specified MoveTo position, otherwise False.

Source code in TPTBox/spine/spinestats/angles.py

def has_point(self, v: Vertebra_Instance | int, poi: POI) -> bool:
    """Check if the given vertebra has a specific POI at the current MoveTo position.

    Args:
        v (Vertebra_Instance | int): The vertebra instance or its ID.
        poi (POI): The point of interest data structure.

    Returns:
        bool: True if the POI exists for the specified MoveTo position, otherwise False.
    """
    if isinstance(v, int):
        v = Vertebra_Instance(v)
    try:
        self.get_point(v, poi)
    except KeyError:
        return False
    return True

get_location

get_location(v: Vertebra_Instance | int, poi: POI) -> tuple

Determine the anatomical POI coordinates for the current MoveTo position in a vertebra.

Parameters:

Name	Type	Description	Default
v	Vertebra_Instance | int	The vertebra instance or its ID.	required
poi	POI	The point of interest data structure.	required

Returns:

Name	Type	Description
tuple	tuple	The location tuple that defines the position within the vertebra.

Source code in TPTBox/spine/spinestats/angles.py

def get_location(self, v: Vertebra_Instance | int, poi: POI) -> tuple:
    """Determine the anatomical POI coordinates for the current MoveTo position in a vertebra.

    Args:
        v (Vertebra_Instance | int): The vertebra instance or its ID.
        poi (POI): The point of interest data structure.

    Returns:
        tuple: The location tuple that defines the position within the vertebra.
    """
    if isinstance(v, int):
        v = Vertebra_Instance(v)
    if self == self.CENTER:
        return (v, 50)
    elif self == self.BOTTOM:
        # Test IVD
        subreg = Location.Vertebra_Disc
        if (v, subreg) in poi:
            return (v, subreg)
        # Test if it has next
        subreg = Location.Additional_Vertebral_Body_Middle_Inferior_Median
        if (v, subreg) in poi:
            return (v, subreg)
        # Test if it has next POINT
        next_vert = v.get_next_poi(poi)
        if next_vert is not None and (next_vert, 50) in poi:
            return (v, subreg, next_vert, subreg)
    elif self == self.TOP:
        prev_vert = v.get_previous_poi(poi)
        # Test IVD
        subreg = Location.Vertebra_Disc
        if prev_vert is not None and (prev_vert, subreg) in poi:
            return (prev_vert, subreg)
        # Test if it has next
        subreg = Location.Dens_axis
        if (v, subreg) in poi:
            return (v, subreg)
        # Test if it has next
        subreg = Location.Additional_Vertebral_Body_Middle_Superior_Median
        if (v, subreg) in poi:
            return (v, subreg)
        # Test if it has next POINT
        if prev_vert is not None and (prev_vert, 50) in poi:
            return (v, subreg, prev_vert, subreg)
    return (v, 50)

get_point

get_point(v: Vertebra_Instance | int, poi: POI) -> np.ndarray

Retrieve the 3D coordinates of a POI in a vertebra at the current MoveTo position.

Parameters:

Name	Type	Description	Default
v	Vertebra_Instance | int	The vertebra instance or its ID.	required
poi	POI	The point of interest data structure.	required

Returns:

Type	Description
ndarray	np.ndarray: The 3D coordinates of the specified POI.

Raises:

Type	Description
NotImplementedError	If the POI cannot be determined.

Source code in TPTBox/spine/spinestats/angles.py

def get_point(self, v: Vertebra_Instance | int, poi: POI) -> np.ndarray:
    """Retrieve the 3D coordinates of a POI in a vertebra at the current MoveTo position.

    Args:
        v (Vertebra_Instance | int): The vertebra instance or its ID.
        poi (POI): The point of interest data structure.

    Returns:
        np.ndarray: The 3D coordinates of the specified POI.

    Raises:
        NotImplementedError: If the POI cannot be determined.
    """
    a = self.get_location(v, poi)
    if len(a) == 2:
        return np.array(poi[a])
    elif len(a) == 4:
        return (np.array(poi[a[0], a[1]]) + np.array(poi[a[2], a[3]])) / 2
    raise NotImplementedError(v, poi)

unit_vector

unit_vector(vector: ndarray) -> np.ndarray

Return the unit vector of the input vector.

Parameters:

Name	Type	Description	Default
vector	ndarray	Any non-zero numeric array.	required

Returns:

Type	Description
ndarray	Array with the same direction as vector but unit length.

Source code in TPTBox/spine/spinestats/angles.py

def unit_vector(vector: np.ndarray) -> np.ndarray:
    """Return the unit vector of the input vector.

    Args:
        vector: Any non-zero numeric array.

    Returns:
        Array with the same direction as ``vector`` but unit length.
    """
    return vector / np.linalg.norm(vector)

angle_between

angle_between(v1, v2) -> float

Calculates the angle in radians between two vectors.

Parameters:

Name	Type	Description	Default
v1	tuple	The first vector.	required
v2	tuple	The second vector.	required

Returns:

Name	Type	Description
float	float	The angle in radians between vectors 'v1' and 'v2'.

Examples:

>>> angle_between((1, 0, 0), (0, 1, 0))
1.5707963267948966
>>> angle_between((1, 0, 0), (1, 0, 0))
0.0
>>> angle_between((1, 0, 0), (-1, 0, 0))
3.141592653589793

Source code in TPTBox/spine/spinestats/angles.py

def angle_between(v1, v2) -> float:
    """Calculates the angle in radians between two vectors.

    Args:
        v1 (tuple): The first vector.
        v2 (tuple): The second vector.

    Returns:
        float: The angle in radians between vectors 'v1' and 'v2'.

    Examples:
        >>> angle_between((1, 0, 0), (0, 1, 0))
        1.5707963267948966
        >>> angle_between((1, 0, 0), (1, 0, 0))
        0.0
        >>> angle_between((1, 0, 0), (-1, 0, 0))
        3.141592653589793
    """
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

get_to_space

get_to_space(a, b, c) -> tuple[np.ndarray, np.ndarray]

Compute forward and inverse transformation matrices for the space defined by three orthogonal vectors.

Parameters:

Name	Type	Description	Default
a	ndarray	First orthogonal vector.	required
b	ndarray	Second orthogonal vector.	required
c	ndarray	Third orthogonal vector.	required

Returns:

Name	Type	Description
tuple	tuple[ndarray, ndarray]	A tuple containing two matrices (to_space, from_space): - to_space (np.ndarray): Transformation matrix to the canonical space. - from_space (np.ndarray): Transformation matrix from the canonical space.

Source code in TPTBox/spine/spinestats/angles.py

def get_to_space(a, b, c) -> tuple[np.ndarray, np.ndarray]:
    """Compute forward and inverse transformation matrices for the space defined by three orthogonal vectors.

    Args:
        a (np.ndarray): First orthogonal vector.
        b (np.ndarray): Second orthogonal vector.
        c (np.ndarray): Third orthogonal vector.

    Returns:
        tuple: A tuple containing two matrices (to_space, from_space):
            - to_space (np.ndarray): Transformation matrix to the canonical space.
            - from_space (np.ndarray): Transformation matrix from the canonical space.
    """
    from_space = np.stack([a, b, c], axis=1)
    to_space = np.linalg.inv(from_space)
    return to_space, from_space

cosine_distance

cosine_distance(a, b) -> float

Computes the cosine distance between two vectors.

Parameters:

Name	Type	Description	Default
a	ndarray	The first vector.	required
b	ndarray	The second vector.	required

Returns:

Name	Type	Description
float	float	The cosine distance between vectors 'a' and 'b'.

Source code in TPTBox/spine/spinestats/angles.py

def cosine_distance(a, b) -> float:
    """Computes the cosine distance between two vectors.

    Args:
        a (np.ndarray): The first vector.
        b (np.ndarray): The second vector.

    Returns:
        float: The cosine distance between vectors 'a' and 'b'.
    """
    return np.dot(a, b) / (np.linalg.norm(b) * np.linalg.norm(a))

compute_angel_between_two_points_

compute_angel_between_two_points_(poi: POI, vert_id1: Vertebra_Instance | int | None, vert_id2: Vertebra_Instance | int | None, direction: DIRECTIONS, vert_id1_mv: MoveTo = MoveTo.CENTER, vert_id2_mv: MoveTo = MoveTo.CENTER, project_2D=False, use_ivd_direction=False) -> float | None

Compute the 2D or 3D angle between two anatomical landmarks.

Useful for calculating coplanar angles, lordosis, and kyphosis depending on the direction specified.

Parameters:

Name	Type	Description	Default
poi	POI	An object representing a point of interest that supports indexing with point IDs and returns 3D coordinates. Must have methods reorient_() and rescale_() to prepare the data.	required
vert_id1	Vertebra_Instance | int | None	The identifier for the first point of interest.	required
vert_id2	Vertebra_Instance | int | None	The identifier for the second point of interest.	required
direction	DIRECTIONS	The direction in which to compute the angle. Possible values are: - "P" for Posterior, used for calculating lordosis and kyphosis. - "A" for Anterior. - "R" for Right, used for calculating coplanar angles. - "L" for Left. - "S" for Superior. - "I" for Inferior.	required
vert_id1_mv	MoveTo	MoveTo instance indicating the position to consider for the first vertebra. Defaults to MoveTo.CENTER.	CENTER
vert_id2_mv	MoveTo	MoveTo instance indicating the position to consider for the second vertebra. Defaults to MoveTo.CENTER.	CENTER
project_2d	bool	If True, computes the 2D projection of the angle. Defaults to False.	required

Returns:

Type	Description
float | None	float | None: The computed angle in degrees. Returns None if either vertebra ID is invalid.

Raises:

Type	Description
NotImplementedError	If the direction direction is not one of the recognized values.

Example

To compute the coplanar angle between two vertebrae with IDs 20 and 21:

compute_angel_between_two_points_(poi, 20, 21, "R")

To compute the lordosis angle between two vertebrae with IDs 20 and 21:

compute_angel_between_two_points_(poi, 20, 21, "P")

Source code in TPTBox/spine/spinestats/angles.py

def compute_angel_between_two_points_(
    poi: POI,
    vert_id1: Vertebra_Instance | int | None,
    vert_id2: Vertebra_Instance | int | None,
    direction: DIRECTIONS,
    vert_id1_mv: MoveTo = MoveTo.CENTER,
    vert_id2_mv: MoveTo = MoveTo.CENTER,
    project_2D=False,
    use_ivd_direction=False,
) -> float | None:
    """Compute the 2D or 3D angle between two anatomical landmarks.

    Useful for calculating coplanar angles, lordosis, and kyphosis depending on
    the direction specified.

    Args:
        poi (POI): An object representing a point of interest that supports indexing
            with point IDs and returns 3D coordinates. Must have methods `reorient_()`
            and `rescale_()` to prepare the data.
        vert_id1 (Vertebra_Instance | int | None): The identifier for the first point of interest.
        vert_id2 (Vertebra_Instance | int | None): The identifier for the second point of interest.
        direction (DIRECTIONS): The direction in which to compute the angle. Possible values are:
            - "P" for Posterior, used for calculating lordosis and kyphosis.
            - "A" for Anterior.
            - "R" for Right, used for calculating coplanar angles.
            - "L" for Left.
            - "S" for Superior.
            - "I" for Inferior.
        vert_id1_mv (MoveTo, optional): MoveTo instance indicating the position to consider for the first vertebra. Defaults to MoveTo.CENTER.
        vert_id2_mv (MoveTo, optional): MoveTo instance indicating the position to consider for the second vertebra. Defaults to MoveTo.CENTER.
        project_2d (bool, optional): If True, computes the 2D projection of the angle. Defaults to False.

    Returns:
        float | None: The computed angle in degrees. Returns None if either vertebra ID is invalid.

    Raises:
        NotImplementedError: If the direction `direction` is not one of the recognized values.

    Example:
        To compute the coplanar angle between two vertebrae with IDs 20 and 21:

        >>> compute_angel_between_two_points_(poi, 20, 21, "R")

        To compute the lordosis angle between two vertebrae with IDs 20 and 21:

        >>> compute_angel_between_two_points_(poi, 20, 21, "P")
    """
    if vert_id1 is None or vert_id2 is None:
        return None
    # assert project_2D, "project_2D == True"
    id1: int = vert_id1.value if isinstance(vert_id1, Enum) else vert_id1
    id2: int = vert_id2.value if isinstance(vert_id2, Enum) else vert_id2
    if (id1, 50) not in poi or (id2, 50) not in poi:
        return None
    assert id1 != id2, id1

    # Ensure id1 is less than id2
    if id1 > id2:
        id1, id2 = id2, id1
    # Reorient and rescale the POI data
    poi.reorient_().rescale_()
    recompute_use_ivd_direction = False
    # Determine direction-specific settings
    location2 = None
    if direction in ["P", "A"]:
        # Note: inf - The value of the direction does noting
        location = Location.Vertebra_Direction_Posterior
        inv = 1 if direction == "P" else -1
    elif direction in ["R", "L"]:
        location = Location.Vertebra_Direction_Right
        inv = 1 if direction == "R" else -1
        if use_ivd_direction:
            location = Location.Vertebra_Disc_Inferior if id1 > IVD_MORE_ACCURATE else Location.Vertebra_Direction_Right
            location2 = Location.Vertebra_Disc_Inferior if id2 > IVD_MORE_ACCURATE else Location.Vertebra_Direction_Right
            recompute_use_ivd_direction = True
    elif direction in ["S", "I"]:
        if use_ivd_direction:
            location = Location.Vertebra_Disc_Inferior if id1 > IVD_MORE_ACCURATE else Location.Vertebra_Direction_Inferior
            location2 = Location.Vertebra_Disc_Inferior if id2 > IVD_MORE_ACCURATE else Location.Vertebra_Direction_Inferior
        else:
            location = Location.Vertebra_Direction_Inferior
        inv = 1 if direction == "I" else -1
    else:
        raise NotImplementedError(f"Direction '{direction}' is not recognized.")
        # Calculate normals for the vertebrae
    if location2 is None:
        location2 = location
    norm1_vert = _get_norm(poi, id1, vert_id1_mv, location, inv=inv)
    norm2_vert = _get_norm(poi, id2, vert_id2_mv, location2, inv=inv)
    if norm1_vert is None or norm2_vert is None:
        return None
    if recompute_use_ivd_direction:
        # Compute right from Post (Vert) + Inferior (Disc)
        if Location.Vertebra_Disc_Inferior == location:
            norm1_post = _get_norm(poi, id1, vert_id1_mv, Location.Vertebra_Direction_Posterior)
            if norm1_post is None:
                return None

            norm1_vert = np.cross(norm1_vert, norm1_post)
        if Location.Vertebra_Disc_Inferior == location2:
            norm2_post = _get_norm(poi, id2, vert_id2_mv, Location.Vertebra_Direction_Posterior)
            if norm2_post is None:
                return None
            norm2_vert = np.cross(norm2_vert, norm2_post)

    assert norm1_vert is not None
    assert norm2_vert is not None
    # if _debug_plot is not None:
    #    nii_old = _debug_plot.copy()
    #    _debug_plot.reorient_().rescale_()
    #    _debug_plot = _debug_plot * 0
    #    plot_ray(poi[id1, subreg], norm1_vert, _debug_plot, inplace=True, value=1)
    #    plot_ray(poi[id2, subreg], norm2_vert, _debug_plot, inplace=True, value=2)
    # Calculate the 3D angle between the normal
    if not project_2D:
        angle_3D = angle_between(norm1_vert, norm2_vert)  # noqa: N806
        return angle_3D / np.pi * 180
    p1 = vert_id1_mv.get_point(id1, poi)
    p2 = vert_id2_mv.get_point(id2, poi)

    if direction in ["S", "I"]:
        a = unit_vector(np.array(p1) - np.array(poi[id1, Location.Vertebra_Direction_Right]))
        b = unit_vector(np.array(p2) - np.array(poi[id2, Location.Vertebra_Direction_Right]))
        norm_down = (a + b) / 2
    else:
        norm_down = unit_vector(np.array(p2) - np.array(p1))
    # Compute cross products to isolate relevant components
    norm_to_remove = np.cross(norm_down, (norm1_vert + norm2_vert))
    norm_keep_other = np.cross(norm_to_remove, norm_down)
    # Bring into a space where we can remove not considered component (like anterior/posterior is ignored for copangles)
    # Transform into a 2D space for angle calculation
    to_space, from_space = get_to_space(norm_keep_other, norm_down, norm_to_remove)
    norm1_vert = to_space @ norm1_vert
    norm2_vert = to_space @ norm2_vert
    norm1_vert[2] = 0
    norm2_vert[2] = 0
    norm1_vert = from_space @ norm1_vert
    norm2_vert = from_space @ norm2_vert
    # Calculate the 2D angle between the transformed normals
    angle_2D = angle_between(norm1_vert, norm2_vert)  # noqa: N806
    # if _debug_plot is not None:
    #    plot_ray(poi[id1, subreg], norm1_vert, _debug_plot, inplace=True, value=3)
    #    plot_ray(poi[id2, subreg], norm2_vert, _debug_plot, inplace=True, value=4)

    #    plot_ray(poi[24, 50], norm_down, _debug_plot, inplace=True, value=5)
    #    plot_ray(poi[24, 50], norm_to_removed, _debug_plot, inplace=True, value=6)
    #    plot_ray(poi[24, 50], norm_keep_other, _debug_plot, inplace=True, value=7)
    #    _debug_plot.dilate_msk_().resample_from_to_(nii_old).save(OURPATH)
    return angle_2D / np.pi * 180

compute_lordosis_and_kyphosis

compute_lordosis_and_kyphosis(poi: POI, project_2D=True) -> dict[str, float | None]

Calculates the angles of cervical lordosis, thoracic kyphosis, and lumbar lordosis based on the given points of interest (POI).

This function determines the angles formed by specific vertebrae along the spine, which are indicative of spinal curvatures. The angles are calculated for three key regions: cervical, thoracic, and lumbar, representing lordosis and kyphosis.

Parameters:

Name	Type	Description	Default
poi	POI	The points of interest object containing 3D coordinates for various vertebrae. It must include the vertebra direction information for proper calculation. (Location.Vertebra_Direction_Posterior)	required
project_2d	bool	If True, the calculation is done in 2D projection; otherwise, in 3D. Defaults to False.	required

Returns:

Name	Type	Description
dict	dict[str, float | None]	A dictionary containing the following key-value pairs: - "cervical_lordosis": The angle of cervical lordosis, calculated between C2 and C7. - "thoracic_kyphosis": The angle of thoracic kyphosis, calculated between T1 and the last thoracic vertebra. - "lumbar_lordosis": The angle of lumbar lordosis, calculated between L1 and the last lumbar vertebra.

Raises:

Type	Description
AssertionError	If the required vertebra direction information is not present in the POI.

Notes

• It is essential that the poi contains the posterior vertebra direction for accurate angle calculations.
• Thoracic kyphosis is calculated from T1 to the last thoracic vertebra identified in the POI.
• Lumbar lordosis is calculated from L1 to the last lumbar vertebra identified in the POI.

Example

To compute the spinal angles for a given POI object:

angles = compute_lordosis_and_kyphosis(poi, project_2d=True) print(angles)

Source code in TPTBox/spine/spinestats/angles.py

def compute_lordosis_and_kyphosis(poi: POI, project_2D=True) -> dict[str, float | None]:
    """Calculates the angles of cervical lordosis, thoracic kyphosis, and lumbar lordosis based on the given points of interest (POI).

    This function determines the angles formed by specific vertebrae along the spine, which are indicative of spinal curvatures.
    The angles are calculated for three key regions: cervical, thoracic, and lumbar, representing lordosis and kyphosis.

    Args:
        poi (POI): The points of interest object containing 3D coordinates for various vertebrae. It must include
            the vertebra direction information for proper calculation. (Location.Vertebra_Direction_Posterior)
        project_2d (bool): If True, the calculation is done in 2D projection; otherwise, in 3D. Defaults to False.

    Returns:
        dict: A dictionary containing the following key-value pairs:
            - "cervical_lordosis": The angle of cervical lordosis, calculated between C2 and C7.
            - "thoracic_kyphosis": The angle of thoracic kyphosis, calculated between T1 and the last thoracic vertebra.
            - "lumbar_lordosis": The angle of lumbar lordosis, calculated between L1 and the last lumbar vertebra.

    Raises:
        AssertionError: If the required vertebra direction information is not present in the POI.

    Notes:
        - It is essential that the `poi` contains the posterior vertebra direction for accurate angle calculations.
        - Thoracic kyphosis is calculated from T1 to the last thoracic vertebra identified in the POI.
        - Lumbar lordosis is calculated from L1 to the last lumbar vertebra identified in the POI.

    Example:
        To compute the spinal angles for a given POI object:

        >>> angles = compute_lordosis_and_kyphosis(poi, project_2d=True)
        >>> print(angles)
        {'cervical_lordosis': 30.5, 'thoracic_kyphosis': 35.0, 'lumbar_lordosis': 45.2}
    """
    assert Location.Vertebra_Direction_Posterior.value in poi.keys_subregion(), (
        "You need to compute the Direction in the Poi (Location.Vertebra_Direction_Posterior)"
    )
    last_t = _get_last_thoracic(poi)
    last_l = _get_last_lumbar(poi)
    poi = poi.copy()
    cervical = compute_angel_between_two_points_(
        poi,
        Vertebra_Instance.C2,
        Vertebra_Instance.C7,
        "P",
        MoveTo.TOP,
        MoveTo.BOTTOM,
        project_2D,
    )
    thoracic = compute_angel_between_two_points_(poi, Vertebra_Instance.T1, last_t, "P", MoveTo.TOP, MoveTo.BOTTOM, project_2D)
    lumbar = compute_angel_between_two_points_(poi, Vertebra_Instance.L1, last_l, "P", MoveTo.TOP, MoveTo.BOTTOM, project_2D)
    return {
        "cervical_lordosis": cervical,
        "thoracic_kyphosis": thoracic,
        "lumbar_lordosis": lumbar,
    }

compute_max_cobb_angle

compute_max_cobb_angle(poi: POI, vertebrae_list=None, vert_id1_mv: MoveTo = MoveTo.TOP, vert_id2_mv: MoveTo = MoveTo.BOTTOM, project_2D=True, use_ivd_direction=False) -> tuple[float, int | None, int | None, int | None]

Calculates the maximum Cobb angle from a list of vertebrae using the points of interest (POI).

The Cobb angle is a measure commonly used to quantify the degree of spinal curvature, particularly for scoliosis. This function identifies the maximum Cobb angle by comparing angles between pairs of vertebrae in the specified list. You must compute have computed pois in the following structures: Version 1: poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right]) ivd position will be interpolated. Version 2: poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc]) ivd (Vertebra_Disc) will be computed by the segmentation Version 3 (best): poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc,Location.Vertebra_Disc_Superior]) ivd (Vertebra_Disc) will be computed by the segmentation + use_ivd_direction = True will use the disc direction and will note if there is a large shift between vertebra without rotation.

Parameters:

Name	Type	Description	Default
poi	POI	The points of interest object containing 3D coordinates for various vertebrae.	required
vertebrae_list	list	A list of vertebra instances to consider for Cobb angle calculation. If not provided, defaults to all cervical, thoracic, and lumbar vertebrae.	None
vert_id1_mv	MoveTo	Enum indicating the move direction for the first vertebra (default is MoveTo.TOP).	TOP
vert_id2_mv	MoveTo	Enum indicating the move direction for the second vertebra (default is MoveTo.BOTTOM).	BOTTOM
project_2d	bool	If True, the calculation is done in 2D projection; otherwise, in 3D. Defaults to False.	required

Returns:

Name	Type	Description
tuple	tuple[float, int | None, int | None, int | None]	A tuple containing the following elements: - max_angle (float): The maximum Cobb angle identified between any pair of vertebrae. - from_vert (int or None): The vertebra ID at which the maximum angle originates. - to_vert (int or None): The vertebra ID at which the maximum angle terminates. - apex (int or None): The vertebra ID that is the apex of the maximum Cobb angle curvature.

Raises:

Type	Description
AssertionError	If the necessary direction data for computation is not present in the POI.

Notes

• The function iterates through pairs of vertebrae to calculate the angles between them.
• It uses the compute_angel_between_two_points_ function to determine the angle between two vertebrae.
• The apex is determined by finding the vertebra with the largest cosine distance to the calculated apex vector.
• Ensure that the vertebrae are provided in a correct anatomical order to avoid inaccurate results.

Example

To compute the maximum Cobb angle for a given POI object:

max_angle, from_vert, to_vert, apex = compute_max_cobb_angle(poi, project_2d=True) print(f"Max Angle: {max_angle}, From: {from_vert}, To: {to_vert}, Apex: {apex}") Max Angle: 35.6, From: 3, To: 12, Apex: 7

Source code in TPTBox/spine/spinestats/angles.py

def compute_max_cobb_angle(
    poi: POI,
    vertebrae_list=None,
    vert_id1_mv: MoveTo = MoveTo.TOP,
    vert_id2_mv: MoveTo = MoveTo.BOTTOM,
    project_2D=True,
    use_ivd_direction=False,
) -> tuple[float, int | None, int | None, int | None]:
    """Calculates the maximum Cobb angle from a list of vertebrae using the points of interest (POI).

    The Cobb angle is a measure commonly used to quantify the degree of spinal curvature, particularly for scoliosis.
    This function identifies the maximum Cobb angle by comparing angles between pairs of vertebrae in the specified list.
    You must compute have computed pois in the following structures:
    Version 1:
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right])
        ivd position will be interpolated.
    Version 2:
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc])
        ivd (Vertebra_Disc) will be computed by the segmentation
    Version 3 (best):
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc,Location.Vertebra_Disc_Superior])
        ivd (Vertebra_Disc) will be computed by the segmentation
        + use_ivd_direction = True will use the disc direction and will note if there is a large shift between vertebra without rotation.

    Args:
        poi (POI): The points of interest object containing 3D coordinates for various vertebrae.
        vertebrae_list (list, optional): A list of vertebra instances to consider for Cobb angle calculation.
            If not provided, defaults to all cervical, thoracic, and lumbar vertebrae.
        vert_id1_mv (MoveTo): Enum indicating the move direction for the first vertebra (default is MoveTo.TOP).
        vert_id2_mv (MoveTo): Enum indicating the move direction for the second vertebra (default is MoveTo.BOTTOM).
        project_2d (bool): If True, the calculation is done in 2D projection; otherwise, in 3D. Defaults to False.

    Returns:
        tuple: A tuple containing the following elements:
            - max_angle (float): The maximum Cobb angle identified between any pair of vertebrae.
            - from_vert (int or None): The vertebra ID at which the maximum angle originates.
            - to_vert (int or None): The vertebra ID at which the maximum angle terminates.
            - apex (int or None): The vertebra ID that is the apex of the maximum Cobb angle curvature.

    Raises:
        AssertionError: If the necessary direction data for computation is not present in the POI.

    Notes:
        - The function iterates through pairs of vertebrae to calculate the angles between them.
        - It uses the `compute_angel_between_two_points_` function to determine the angle between two vertebrae.
        - The apex is determined by finding the vertebra with the largest cosine distance to the calculated apex vector.
        - Ensure that the vertebrae are provided in a correct anatomical order to avoid inaccurate results.

    Example:
        To compute the maximum Cobb angle for a given POI object:

        >>> max_angle, from_vert, to_vert, apex = compute_max_cobb_angle(poi, project_2d=True)
        >>> print(f"Max Angle: {max_angle}, From: {from_vert}, To: {to_vert}, Apex: {apex}")
        Max Angle: 35.6, From: 3, To: 12, Apex: 7
    """
    max_angle = 0
    from_vert: int | None = None
    to_vert: int | None = None
    if vertebrae_list is None:
        # Define the range of vertebrae to consider (e.g., from T1 to L5)
        vertebrae_list = list(Vertebra_Instance.cervical()) + list(Vertebra_Instance.thoracic()) + list(Vertebra_Instance.lumbar())
        vertebrae_list = vertebrae_list[vertebrae_list.index(VERT_START_COBB) :]
    # Iterate through pairs of adjacent vertebrae
    for i in range(len(vertebrae_list)):
        vert_id1: int = vertebrae_list[i].value
        for i2 in range(i):
            vert_id2: int = vertebrae_list[i2].value
            # Ensure that both vertebrae are present in the POI
            if vert_id1 in poi.keys_region() and vert_id2 in poi.keys_region():
                # Compute cobblanar angles for both right and left directions
                angle = compute_angel_between_two_points_(
                    poi,
                    vert_id1,
                    vert_id2,
                    "R",
                    vert_id1_mv,
                    vert_id2_mv,
                    project_2D,
                    use_ivd_direction=use_ivd_direction,
                )
                if angle is None:
                    continue
                # Update max_angle if a larger angle is found
                if max_angle < angle:
                    max_angle, from_vert, to_vert = angle, vert_id2, vert_id1
    apex: int | None = None
    cos_dis = 0
    if from_vert is not None and to_vert is not None:
        a = _get_norm(poi, from_vert, vert_id1_mv, Location.Vertebra_Direction_Right, 1)
        b = _get_norm(poi, to_vert, vert_id2_mv, Location.Vertebra_Direction_Right, 1)
        assert a is not None
        assert b is not None
        apex_v = (a + b) / 2
        for i in vertebrae_list[vertebrae_list.index(Vertebra_Instance(from_vert)) : vertebrae_list.index(Vertebra_Instance(to_vert)) + 1]:
            try:
                a = _get_norm(poi, i, vert_id2_mv, Location.Vertebra_Direction_Right, 1)
                if a is None:
                    continue
            except KeyError:
                continue
            cos_new = cosine_distance(a, apex_v)
            if cos_dis < cos_new:
                cos_dis = cos_new
                apex = i.value
    return max_angle, from_vert, to_vert, apex

compute_max_cobb_angle_multi

compute_max_cobb_angle_multi(poi: POI, vertebrae_list=None, threshold_deg=10, out_list=None, vert_id1_mv: MoveTo = MoveTo.TOP, vert_id2_mv: MoveTo = MoveTo.BOTTOM, use_ivd_direction=False, project_2D=True) -> list[tuple[float, int, int, int | None]]

Identifies multiple Cobb angles along the spine that exceed a given threshold.

This function calculates Cobb angles for a list of vertebrae and recursively finds multiple spinal curvatures that are large enough, as determined by a threshold angle. It is useful for detecting and evaluating scoliosis or other spinal deformities with multiple curves. You must compute have computed pois in the following structures: Version 1: poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right]) ivd position will be interpolated. Version 2: poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc]) ivd (Vertebra_Disc) will be computed by the segmentation Version 3 (best): poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc,Location.Vertebra_Disc_Inferior]) ivd (Vertebra_Disc) will be computed by the segmentation + use_ivd_direction = True will use the disc direction and will note if there is a large shift between vertebra without rotation.

Parameters:

Name	Type	Description	Default
poi	POI	The points of interest object containing 3D coordinates for various vertebrae.	required
vertebrae_list	list	A list of vertebra instances to consider for Cobb angle calculation. If not provided, defaults to all cervical, thoracic, and lumbar vertebrae.	None
threshold_deg	float	The angle threshold in degrees. Only curves with angles greater than this value will be recorded.	10
out_list	list	A list to store the results of the identified Cobb angles. If not provided, an empty list will be created and used.	None
vert_id1_mv	MoveTo	Enum indicating the move direction for the first vertebra (default is MoveTo.TOP).	TOP
vert_id2_mv	MoveTo	Enum indicating the move direction for the second vertebra (default is MoveTo.BOTTOM).	BOTTOM
use_ivd_direction		Uses the IVD direction instead of the Vertebra direction for Lumbar and Thorax region.	False

Returns:

Name	Type	Description
list	list[tuple[float, int, int, int | None]]	A list of tuples, each containing: - max_angle (float): The Cobb angle identified between two vertebrae that exceeds the threshold. - from_vert (int): The vertebra ID at which the angle originates. - to_vert (int): The vertebra ID at which the angle terminates. - apex (int or None): The vertebra ID that is the apex of the curvature.

Notes

• The function splits the list of vertebrae and recursively calculates Cobb angles for sublists.
• Only angles that are above the specified threshold are added to the output list.
• It uses compute_max_cobb_angle to find the maximum Cobb angle within a given set of vertebrae.

Example

To compute multiple Cobb angles for a given POI object with a threshold of 15 degrees:

curves = compute_max_cobb_angle_multi(poi, threshold_deg=15) for curve in curves: print(f"Angle: {curve[0]}, From: {curve[1]}, To: {curve[2]}, Apex: {curve[3]}") Angle: 18.2, From: 2, To: 6, Apex: 4 Angle: 12.4, From: 7, To: 10, Apex: 8

Source code in TPTBox/spine/spinestats/angles.py

def compute_max_cobb_angle_multi(
    poi: POI,
    vertebrae_list=None,
    threshold_deg=10,
    out_list=None,
    vert_id1_mv: MoveTo = MoveTo.TOP,
    vert_id2_mv: MoveTo = MoveTo.BOTTOM,
    use_ivd_direction=False,
    project_2D=True,
) -> list[tuple[float, int, int, int | None]]:
    """Identifies multiple Cobb angles along the spine that exceed a given threshold.

    This function calculates Cobb angles for a list of vertebrae and recursively finds multiple
    spinal curvatures that are large enough, as determined by a threshold angle. It is useful for
    detecting and evaluating scoliosis or other spinal deformities with multiple curves.
    You must compute have computed pois in the following structures:
    Version 1:
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right])
        ivd position will be interpolated.
    Version 2:
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc])
        ivd (Vertebra_Disc) will be computed by the segmentation
    Version 3 (best):
        poi = calc_poi_from_subreg_vert(nii, nii_subreg, subreg_id=[Location.Vertebra_Corpus,Location.Vertebra_Direction_Right,Location.Vertebra_Disc,Location.Vertebra_Disc_Inferior])
        ivd (Vertebra_Disc) will be computed by the segmentation
        + use_ivd_direction = True will use the disc direction and will note if there is a large shift between vertebra without rotation.

    Args:
        poi (POI): The points of interest object containing 3D coordinates for various vertebrae.
        vertebrae_list (list, optional): A list of vertebra instances to consider for Cobb angle calculation.
            If not provided, defaults to all cervical, thoracic, and lumbar vertebrae.
        threshold_deg (float): The angle threshold in degrees. Only curves with angles greater than
            this value will be recorded.
        out_list (list, optional): A list to store the results of the identified Cobb angles.
            If not provided, an empty list will be created and used.
        vert_id1_mv (MoveTo): Enum indicating the move direction for the first vertebra (default is MoveTo.TOP).
        vert_id2_mv (MoveTo): Enum indicating the move direction for the second vertebra (default is MoveTo.BOTTOM).
        use_ivd_direction: Uses the IVD direction instead of the Vertebra direction for Lumbar and Thorax region.

    Returns:
        list: A list of tuples, each containing:
            - max_angle (float): The Cobb angle identified between two vertebrae that exceeds the threshold.
            - from_vert (int): The vertebra ID at which the angle originates.
            - to_vert (int): The vertebra ID at which the angle terminates.
            - apex (int or None): The vertebra ID that is the apex of the curvature.

    Notes:
        - The function splits the list of vertebrae and recursively calculates Cobb angles for sublists.
        - Only angles that are above the specified threshold are added to the output list.
        - It uses `compute_max_cobb_angle` to find the maximum Cobb angle within a given set of vertebrae.

    Example:
        To compute multiple Cobb angles for a given POI object with a threshold of 15 degrees:

        >>> curves = compute_max_cobb_angle_multi(poi, threshold_deg=15)
        >>> for curve in curves:
        >>>     print(f"Angle: {curve[0]}, From: {curve[1]}, To: {curve[2]}, Apex: {curve[3]}")
        Angle: 18.2, From: 2, To: 6, Apex: 4
        Angle: 12.4, From: 7, To: 10, Apex: 8
    """
    if out_list is None:
        out_list = []
    if vertebrae_list is None:
        vertebrae_list = list(Vertebra_Instance.cervical()) + list(Vertebra_Instance.thoracic()) + list(Vertebra_Instance.lumbar())
        vertebrae_list = vertebrae_list[vertebrae_list.index(VERT_START_COBB) :]
    if len(vertebrae_list) <= 2:
        return out_list
    max_angle, from_vert, to_vert, apex = compute_max_cobb_angle(
        poi,
        vertebrae_list=vertebrae_list,
        vert_id1_mv=vert_id1_mv,
        vert_id2_mv=vert_id2_mv,
        use_ivd_direction=use_ivd_direction,
        project_2D=project_2D,
    )  # type: ignore
    # split
    if threshold_deg <= max_angle:
        from_vert: int
        assert from_vert is not None
        assert to_vert is not None
        out_list.append((max_angle, from_vert, to_vert, apex))
        above = vertebrae_list[: vertebrae_list.index(Vertebra_Instance(from_vert))]
        below = vertebrae_list[vertebrae_list.index(Vertebra_Instance(to_vert)) :]
        compute_max_cobb_angle_multi(
            poi,
            above,
            vert_id1_mv=vert_id1_mv,
            vert_id2_mv=vert_id2_mv,
            threshold_deg=threshold_deg,
            out_list=out_list,
            use_ivd_direction=use_ivd_direction,
            project_2D=project_2D,
        )
        compute_max_cobb_angle_multi(
            poi,
            below,
            vert_id1_mv=vert_id1_mv,
            vert_id2_mv=vert_id2_mv,
            threshold_deg=threshold_deg,
            out_list=out_list,
            use_ivd_direction=use_ivd_direction,
            project_2D=project_2D,
        )

    return out_list

plot_compute_lordosis_and_kyphosis

plot_compute_lordosis_and_kyphosis(img_path: str | Path | None, poi: POI, img: Image_Reference, seg: Image_Reference | None = None, line_len=100, project_2D=True) -> tuple[dict[str, float | None], Snapshot_Frame]

Plots and computes the angles of lordosis and kyphosis on a spinal image.

This function calculates cervical lordosis, thoracic kyphosis, and lumbar lordosis angles based on the provided Points of Interest (POI) object. It visualizes these angles on the specified image by drawing line segments corresponding to vertebra orientations, and adds annotations for the calculated angles.

Parameters:

Name	Type	Description	Default
img_path	str | Path | None	Path to save the generated image. If None, the image is not saved.	required
poi	POI	The points of interest object containing 3D coordinates for various vertebrae.	required
img	Image_Reference	The reference image on which to plot the angles and lines.	required
seg	Image_Reference | None	The segmentation image reference. Optional, can be None.	None
line_len	int	The length of the lines representing the vertebrae directions (default is 100).	100

Returns:

Name	Type	Description
tuple	tuple[dict[str, float | None], Snapshot_Frame]	A tuple containing: - out2 (dict): A dictionary with the calculated angles of lordosis and kyphosis, including: - "cervical_lordosis" (float): Angle of cervical lordosis. - "thoracic_kyphosis" (float): Angle of thoracic kyphosis. - "lumbar_lordosis" (float): Angle of lumbar lordosis. - snap (Snapshot_Frame): A Snapshot_Frame object containing the plotted image data.

Notes

• Artificial intervertebral discs (IVD) in not present like for CT.
• The lines drawn indicate vertebra orientations with specified lengths and directions.
• The function also generates text annotations for each computed angle, positioned at the midpoint between the respective vertebrae.
• If an image path is provided, the snapshot is saved as an image file.

Example

To compute and visualize lordosis and kyphosis angles on a given image with a POI:

angles, snapshot = plot_compute_lordosis_and_kyphosis("output_path.png", poi, img, seg) print(angles)

Source code in TPTBox/spine/spinestats/angles.py

def plot_compute_lordosis_and_kyphosis(
    img_path: str | Path | None,
    poi: POI,
    img: Image_Reference,
    seg: Image_Reference | None = None,
    line_len=100,
    project_2D=True,
) -> tuple[dict[str, float | None], Snapshot_Frame]:
    """Plots and computes the angles of lordosis and kyphosis on a spinal image.

    This function calculates cervical lordosis, thoracic kyphosis, and lumbar lordosis angles
    based on the provided Points of Interest (POI) object. It visualizes these angles on the
    specified image by drawing line segments corresponding to vertebra orientations, and adds
    annotations for the calculated angles.

    Args:
        img_path (str | Path | None): Path to save the generated image. If None, the image is not saved.
        poi (POI): The points of interest object containing 3D coordinates for various vertebrae.
        img (Image_Reference): The reference image on which to plot the angles and lines.
        seg (Image_Reference | None): The segmentation image reference. Optional, can be None.
        line_len (int): The length of the lines representing the vertebrae directions (default is 100).

    Returns:
        tuple: A tuple containing:
            - out2 (dict): A dictionary with the calculated angles of lordosis and kyphosis, including:
                - "cervical_lordosis" (float): Angle of cervical lordosis.
                - "thoracic_kyphosis" (float): Angle of thoracic kyphosis.
                - "lumbar_lordosis" (float): Angle of lumbar lordosis.
            - snap (Snapshot_Frame): A Snapshot_Frame object containing the plotted image data.

    Notes:
        - Artificial intervertebral discs (IVD) in not present like for CT.
        - The lines drawn indicate vertebra orientations with specified lengths and directions.
        - The function also generates text annotations for each computed angle, positioned at the
          midpoint between the respective vertebrae.
        - If an image path is provided, the snapshot is saved as an image file.

    Example:
        To compute and visualize lordosis and kyphosis angles on a given image with a POI:

        >>> angles, snapshot = plot_compute_lordosis_and_kyphosis("output_path.png", poi, img, seg)
        >>> print(angles)
        {'cervical_lordosis': 34.5, 'thoracic_kyphosis': 42.7, 'lumbar_lordosis': 50.3}
    """
    poi = poi.reorient().rescale_()
    poi = _add_artificial_ivd(poi)
    out = []
    text_out = []
    last_t = _get_last_thoracic(poi)
    last_l = _get_last_lumbar(poi)
    for id1, vert_id1_mv in [
        (Vertebra_Instance.C2, MoveTo.TOP),
        (Vertebra_Instance.T1, MoveTo.TOP),
        (last_t, MoveTo.BOTTOM),
        (last_l, MoveTo.BOTTOM),
    ]:
        vert_id1_mv: MoveTo
        if id1 is None or (id1.value, 50) not in poi:
            continue
        s = vert_id1_mv.get_location(id1, poi)
        a = _get_norm(poi, id1, vert_id1_mv, Location.Vertebra_Direction_Posterior, 1)
        assert a is not None
        out.append((id1.value, s, (a[0] * line_len, a[1] * line_len)))
        out.append((id1.value, s, (-a[0] * line_len * 3, -a[1] * line_len * 3)))
    out2 = compute_lordosis_and_kyphosis(poi, project_2D=project_2D)
    for (name, v), id1, id2 in zip_strict(
        out2.items(), [Vertebra_Instance.C7, last_t, last_l], [Vertebra_Instance.C2, Vertebra_Instance.C7, last_t]
    ):
        if v is None:
            continue
        if id1 is None or id2 is None or (id1.value, 50) not in poi:
            continue
        vert = round((id1.value + id2.value) / 2)
        while (vert, 50) not in poi and vert != 0:
            vert -= 1
        if (vert, 50) not in poi:
            cord = poi[vert, 50]
            text_out.append((vert, (f"{str(name).split('_')[-1]}: {v:.1f}°", 15, cord[1])))

    poi.info["line_segments_sag"] = out + poi.info.get("line_segments_sag", [])
    poi.info["text_sag"] = text_out + poi.info.get("text_sag", [])
    snap = Snapshot_Frame(img, seg, centroids=poi, show_these_subreg_poi=[100])
    if img_path is not None:
        create_snapshot(img_path, [snap])
    return out2, snap

plot_cobb_angle

plot_cobb_angle(img_path: str | Path | None, poi: POI, img: Image_Reference, seg: Image_Reference | None = None, line_len=100, threshold_deg=10, vert_id1_mv: MoveTo = MoveTo.TOP, vert_id2_mv: MoveTo = MoveTo.BOTTOM, use_ivd_direction=False, project_2D=True) -> tuple[list[tuple[float, int, int, int | None]], Snapshot_Frame]

Visualize Cobb angles on a spinal image by plotting the maximum angles across the spine.

Parameters:

Name	Type	Description	Default
img_path	str | Path | None	The file path where the output image with plotted angles should be saved. If None, the image is not saved.	required
poi	POI	An object representing a point of interest that supports indexing with vertebra identifiers and returns 3D coordinates.	required
img	Image_Reference	The image to plot the cobb angles on.	required
seg	Image_Reference | None	Optional segmentation image to be used in conjunction with the main image.	None
line_len	int	The length of the line segments used to visualize the direction of cobb angles.	100
threshold_deg	int	The angle threshold in degrees above which cobb angles are considered for plotting.	10
vert_id1_mv	MoveTo	The MoveTo option for the first vertebra in each angle calculation.	TOP
vert_id2_mv	MoveTo	The MoveTo option for the second vertebra in each angle calculation.	BOTTOM

Returns:

Name	Type	Description
tuple	tuple[list[tuple[float, int, int, int | None]], Snapshot_Frame]	A tuple containing: - List of angle data and segments for plotting. - Snapshot_Frame object containing the final image with cobb angles plotted.

Notes

• The function uses the maximum cobb angle algorithm to find angles and plots them on the sagittal view.
• It assumes vertebrae are labeled using a standard convention (e.g., C1, T1, L1, etc.).
• Only angles exceeding the specified threshold are plotted.

Example

To plot the cobb angles and save the output image:

plot_cobb_angle("output.png", poi, img, seg, line_len=100, threshold_deg=10)

Source code in TPTBox/spine/spinestats/angles.py

def plot_cobb_angle(
    img_path: str | Path | None,
    poi: POI,
    img: Image_Reference,
    seg: Image_Reference | None = None,
    line_len=100,
    threshold_deg=10,
    vert_id1_mv: MoveTo = MoveTo.TOP,
    vert_id2_mv: MoveTo = MoveTo.BOTTOM,
    use_ivd_direction=False,
    project_2D=True,
) -> tuple[list[tuple[float, int, int, int | None]], Snapshot_Frame]:
    """Visualize Cobb angles on a spinal image by plotting the maximum angles across the spine.

    Args:
        img_path (str | Path | None): The file path where the output image with plotted angles should be saved.
            If None, the image is not saved.
        poi (POI): An object representing a point of interest that supports indexing with vertebra identifiers
            and returns 3D coordinates.
        img (Image_Reference): The image to plot the cobb angles on.
        seg (Image_Reference | None): Optional segmentation image to be used in conjunction with the main image.
        line_len (int): The length of the line segments used to visualize the direction of cobb angles.
        threshold_deg (int): The angle threshold in degrees above which cobb angles are considered for plotting.
        vert_id1_mv (MoveTo): The MoveTo option for the first vertebra in each angle calculation.
        vert_id2_mv (MoveTo): The MoveTo option for the second vertebra in each angle calculation.

    Returns:
        tuple: A tuple containing:
            - List of angle data and segments for plotting.
            - Snapshot_Frame object containing the final image with cobb angles plotted.

    Notes:
        - The function uses the maximum cobb angle algorithm to find angles and plots them on the sagittal view.
        - It assumes vertebrae are labeled using a standard convention (e.g., C1, T1, L1, etc.).
        - Only angles exceeding the specified threshold are plotted.

    Example:
        To plot the cobb angles and save the output image:

        >>> plot_cobb_angle("output.png", poi, img, seg, line_len=100, threshold_deg=10)
    """
    poi = poi.reorient().rescale_()
    poi = _add_artificial_ivd(poi)

    out = []
    text_out = []
    copps = compute_max_cobb_angle_multi(
        poi,
        threshold_deg=threshold_deg,
        vert_id1_mv=vert_id1_mv,
        vert_id2_mv=vert_id2_mv,
        use_ivd_direction=use_ivd_direction,
        project_2D=project_2D,
    )
    for max_angle, from_vert, to_vert, apex in copps:
        if from_vert is not None:
            for id1, mv in zip([from_vert, to_vert], [vert_id1_mv, vert_id2_mv]):
                c = mv.get_location(id1, poi)

                if use_ivd_direction and id1 > IVD_MORE_ACCURATE:
                    norm1_post = _get_norm(poi, id1, mv, Location.Vertebra_Direction_Posterior)
                    a = _get_norm(poi, id1, mv, Location.Vertebra_Disc_Inferior)
                    a = np.cross(a, norm1_post)
                else:
                    a = _get_norm(poi, id1, mv, Location.Vertebra_Direction_Right)

                # print(a, id1, mv, c)

                assert a is not None
                out.append((apex, c, (-a[2] * line_len, a[1] * line_len)))
                out.append((apex, c, (a[2] * line_len, -a[1] * line_len)))
                # a = _get_norm(poi, id1, mv, Location.Vertebra_Disc_Inferior)
                # out.append((apex, c, (a[2] * line_len, -a[1] * line_len)))
        if apex is not None:
            cord = poi[apex, 50]
            s = f"copp angle: {max_angle:.1f}° {Vertebra_Instance(from_vert)} - {Vertebra_Instance(to_vert)}"
            text_out.append((apex, (s, 25, cord[1])))
        poi.info["line_segments_cor"] = out + poi.info.get("line_segments_cor", [])
        poi.info["text_cor"] = text_out + poi.info.get("text_cor", [])
    frame = Snapshot_Frame(
        img,
        seg,
        centroids=poi,
        sagittal=False,
        coronal=True,
        show_these_subreg_poi=[100],
    )
    if img_path is not None:
        create_snapshot(img_path, [frame])
    return copps, frame

plot_cobb_and_lordosis_and_kyphosis

plot_cobb_and_lordosis_and_kyphosis(img_path: str | Path | None, poi: POI, img: Image_Reference, seg: Image_Reference | None = None, line_len=100, threshold_deg=10, project_2D=True) -> tuple[list, dict[str, float | None], list[Snapshot_Frame]]

Plots Cobb angles and lordosis/kyphosis angles on a spinal image.

This function calculates and visualizes both the Cobb angles for spinal curvature and the angles of cervical lordosis, thoracic kyphosis, and lumbar lordosis. It overlays these visualizations on the provided spinal image and can save the resulting image to a specified path.

Parameters:

Name	Type	Description	Default
img_path	str | Path | None	Path to save the generated image. If None, the image is not saved.	required
poi	POI	The points of interest object containing 3D coordinates for various vertebrae.	required
img	Image_Reference	The reference image on which to plot the angles and lines.	required
seg	Image_Reference | None	The segmentation image reference. Optional, can be None.	None
line_len	int	The length of the lines representing the vertebrae directions (default is 100).	100
threshold_deg	int	The threshold angle in degrees to identify significant Cobb angles (default is 10).	10

Returns:

Name	Type	Description
tuple	tuple[list, dict[str, float | None], list[Snapshot_Frame]]	A tuple containing: - out_cobb (list): A list of tuples for each significant Cobb angle found, each with: - max_angle (float): The maximum Cobb angle in the segment. - from_vert (int): The vertebra ID at the start of the Cobb angle measurement. - to_vert (int): The vertebra ID at the end of the Cobb angle measurement. - apex (int | None): The vertebra ID of the apex of the curvature. - out_lak (dict): A dictionary with the calculated angles of lordosis and kyphosis, including: - "cervical_lordosis" (float): Angle of cervical lordosis. - "thoracic_kyphosis" (float): Angle of thoracic kyphosis. - "lumbar_lordosis" (float): Angle of lumbar lordosis. - frames (list): A list containing the generated Snapshot_Frame objects for each plot.

Notes

• This function internally calls plot_cobb_angle and plot_compute_lordosis_and_kyphosis to generate the respective plots.
• The function combines both visualizations into a single output image if a path is specified.
• It effectively allows for simultaneous assessment of scoliosis (via Cobb angles) and sagittal plane curvatures (lordosis and kyphosis).

Example

To visualize and save both Cobb angles and lordosis/kyphosis angles:

cobb_angles, lordosis_kyphosis, frames = plot_cobb_and_lordosis_and_kyphosis( ... "output_path.png", poi, img, seg, line_len=150, threshold_deg=15 ... ) print(cobb_angles) [(22.3, 3, 12, 7), (18.5, 13, 17, 15)] print(lordosis_kyphosis)

Source code in TPTBox/spine/spinestats/angles.py

def plot_cobb_and_lordosis_and_kyphosis(
    img_path: str | Path | None,
    poi: POI,
    img: Image_Reference,
    seg: Image_Reference | None = None,
    line_len=100,
    threshold_deg=10,
    project_2D=True,
) -> tuple[list, dict[str, float | None], list[Snapshot_Frame]]:
    """Plots Cobb angles and lordosis/kyphosis angles on a spinal image.

    This function calculates and visualizes both the Cobb angles for spinal curvature and the angles
    of cervical lordosis, thoracic kyphosis, and lumbar lordosis. It overlays these visualizations
    on the provided spinal image and can save the resulting image to a specified path.

    Args:
        img_path (str | Path | None): Path to save the generated image. If None, the image is not saved.
        poi (POI): The points of interest object containing 3D coordinates for various vertebrae.
        img (Image_Reference): The reference image on which to plot the angles and lines.
        seg (Image_Reference | None): The segmentation image reference. Optional, can be None.
        line_len (int): The length of the lines representing the vertebrae directions (default is 100).
        threshold_deg (int): The threshold angle in degrees to identify significant Cobb angles (default is 10).

    Returns:
        tuple: A tuple containing:
            - out_cobb (list): A list of tuples for each significant Cobb angle found, each with:
                - max_angle (float): The maximum Cobb angle in the segment.
                - from_vert (int): The vertebra ID at the start of the Cobb angle measurement.
                - to_vert (int): The vertebra ID at the end of the Cobb angle measurement.
                - apex (int | None): The vertebra ID of the apex of the curvature.
            - out_lak (dict): A dictionary with the calculated angles of lordosis and kyphosis, including:
                - "cervical_lordosis" (float): Angle of cervical lordosis.
                - "thoracic_kyphosis" (float): Angle of thoracic kyphosis.
                - "lumbar_lordosis" (float): Angle of lumbar lordosis.
            - frames (list): A list containing the generated `Snapshot_Frame` objects for each plot.

    Notes:
        - This function internally calls `plot_cobb_angle` and `plot_compute_lordosis_and_kyphosis`
          to generate the respective plots.
        - The function combines both visualizations into a single output image if a path is specified.
        - It effectively allows for simultaneous assessment of scoliosis (via Cobb angles) and sagittal
          plane curvatures (lordosis and kyphosis).

    Example:
        To visualize and save both Cobb angles and lordosis/kyphosis angles:

        >>> cobb_angles, lordosis_kyphosis, frames = plot_cobb_and_lordosis_and_kyphosis(
        ...     "output_path.png", poi, img, seg, line_len=150, threshold_deg=15
        ... )
        >>> print(cobb_angles)
        [(22.3, 3, 12, 7), (18.5, 13, 17, 15)]
        >>> print(lordosis_kyphosis)
        {'cervical_lordosis': 35.2, 'thoracic_kyphosis': 41.5, 'lumbar_lordosis': 48.1}
    """
    out_cobb, frame1 = plot_cobb_angle(
        None,
        poi,
        img,
        seg,
        line_len=line_len,
        threshold_deg=threshold_deg,
        use_ivd_direction=True,
        project_2D=project_2D,
    )
    out_lak, frame2 = plot_compute_lordosis_and_kyphosis(None, poi, img, seg, line_len=line_len, project_2D=project_2D)
    if img_path is not None:
        create_snapshot(img_path, [frame1, frame2])
    return out_cobb, out_lak, [frame1, frame2]

IVD POIs

TPTBox.spine.spinestats.ivd_pois

strategy_calculate_up_vector

strategy_calculate_up_vector(poi: POI, current_vert: NII, vert_id: int, bb: tuple, log: object = _log) -> POI

Compute IVD superior and inferior extreme points via PCA-based normal estimation.

Fits a principal component axis through the IVD label (vert_id + 100) to determine the superior–inferior direction of the disc, then ray-casts to the extreme points along that axis. The results are stored as :attr:~TPTBox.Location.Vertebra_Disc_Inferior and :attr:~TPTBox.Location.Vertebra_Disc_Superior in poi.

Parameters:

Name	Type	Description	Default
poi	POI	POI object that already contains the IVD centroid for vert_id. Updated in place and returned.	required
current_vert	NII	Cropped vertebra NIfTI used for PCA and ray casting.	required
vert_id	int	Integer vertebra region ID whose IVD label is vert_id + 100.	required
bb	tuple	Bounding-box slices (as returned by :meth:~TPTBox.NII.compute_crop) used to convert local coordinates back to the full image space.	required
log	object	Logger instance used for progress / warning messages.	_log

Returns:

Type	Description
POI	The updated :class:~TPTBox.POI object.

Source code in TPTBox/spine/spinestats/ivd_pois.py

def strategy_calculate_up_vector(
    poi: POI,
    current_vert: NII,
    vert_id: int,
    bb: tuple,
    log: object = _log,
) -> POI:
    """Compute IVD superior and inferior extreme points via PCA-based normal estimation.

    Fits a principal component axis through the IVD label (``vert_id + 100``)
    to determine the superior–inferior direction of the disc, then ray-casts to
    the extreme points along that axis.  The results are stored as
    :attr:`~TPTBox.Location.Vertebra_Disc_Inferior` and
    :attr:`~TPTBox.Location.Vertebra_Disc_Superior` in ``poi``.

    Args:
        poi: POI object that already contains the IVD centroid for ``vert_id``.
            Updated in place and returned.
        current_vert: Cropped vertebra NIfTI used for PCA and ray casting.
        vert_id: Integer vertebra region ID whose IVD label is ``vert_id + 100``.
        bb: Bounding-box slices (as returned by
            :meth:`~TPTBox.NII.compute_crop`) used to convert local
            coordinates back to the full image space.
        log: Logger instance used for progress / warning messages.

    Returns:
        The updated :class:`~TPTBox.POI` object.
    """
    center = to_local_np(Location.Vertebra_Disc, bb, poi, vert_id, log)
    if center is None:
        return poi
    try:
        normal_vector = calculate_pca_normal_np(current_vert.rescale().extract_label(vert_id + 100).get_array(), pca_component=2)
    except ValueError:
        return poi
    extreme_point = center + normal_vector * 10

    axis = poi.get_axis("S")
    s_is_poss = int(poi.orientation[axis] == "S")
    below = int(extreme_point[axis] < center[axis])

    ## INFERIOR ##
    if s_is_poss + below == 1:  # xor
        normal_vector *= -1
    # max_distance_ray_cast_convex
    extreme_point = max_distance_ray_cast_convex(current_vert.extract_label(vert_id + 100), center, normal_vector)
    extreme_point_sup = max_distance_ray_cast_convex(current_vert.extract_label(vert_id + 100), center, -normal_vector)
    # extreme_point = center + normal_vector * 10

    assert extreme_point is not None
    assert extreme_point_sup is not None
    poi[vert_id, Location.Vertebra_Disc_Inferior] = tuple(a.start + b for a, b in zip_strict(bb, extreme_point))
    poi[vert_id, Location.Vertebra_Disc_Superior] = tuple(a.start + b for a, b in zip_strict(bb, extreme_point_sup))

    return poi

compute_fake_ivd

compute_fake_ivd(vert_full: NII, subreg_full: NII, poi: POI, dilate: int = 1) -> NII

Add synthetic intervertebral disc (IVD) labels to the vertebra segmentation.

For every adjacent vertebra pair a synthetic IVD mask is computed using morphological operations on the subregion labels. The disc for vertebra i is assigned label 100 + i in the output. The computation is parallelised over vertebrae using a :class:~concurrent.futures.ProcessPoolExecutor.

Parameters:

Name	Type	Description	Default
vert_full	NII	Vertebra segmentation NIfTI; each vertebra has a unique integer label. Modified in place and returned.	required
subreg_full	NII	Subregion segmentation NIfTI containing spine body labels (e.g. 49, 50).	required
poi	POI	POI object used to ensure that :attr:~TPTBox.Location.Vertebra_Corpus centroids are available.	required
dilate	int	Morphological dilation radius applied when building the IVD hull mask.	1

Returns:

Type	Description
NII	The updated vert_full NIfTI with synthetic IVD labels added.

Source code in TPTBox/spine/spinestats/ivd_pois.py

def compute_fake_ivd(vert_full: NII, subreg_full: NII, poi: POI, dilate: int = 1) -> NII:
    """Add synthetic intervertebral disc (IVD) labels to the vertebra segmentation.

    For every adjacent vertebra pair a synthetic IVD mask is computed using
    morphological operations on the subregion labels.  The disc for vertebra
    ``i`` is assigned label ``100 + i`` in the output.  The computation is
    parallelised over vertebrae using a :class:`~concurrent.futures.ProcessPoolExecutor`.

    Args:
        vert_full: Vertebra segmentation NIfTI; each vertebra has a unique
            integer label.  Modified in place and returned.
        subreg_full: Subregion segmentation NIfTI containing spine body labels
            (e.g. 49, 50).
        poi: POI object used to ensure that
            :attr:`~TPTBox.Location.Vertebra_Corpus` centroids are available.
        dilate: Morphological dilation radius applied when building the IVD
            hull mask.

    Returns:
        The updated ``vert_full`` NIfTI with synthetic IVD labels added.
    """
    subreg_ids_required_for_ivd_generation = [
        # Location.Inferior_Articular_Left,
        # Location.Inferior_Articular_Right,
        # Location.Vertebra_Direction_Inferior,
        Location.Vertebra_Corpus,
        # Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Right,
        # Location.Ligament_Attachment_Point_Anterior_Longitudinal_Superior_Left,
        # Location.Ligament_Attachment_Point_Posterior_Longitudinal_Superior_Right,
        # Location.Ligament_Attachment_Point_Posterior_Longitudinal_Superior_Left,
        # Location.Ligament_Attachment_Point_Anterior_Longitudinal_Inferior_Right,
        # Location.Ligament_Attachment_Point_Anterior_Longitudinal_Inferior_Left,
        # Location.Ligament_Attachment_Point_Posterior_Longitudinal_Inferior_Right,
        # Location.Ligament_Attachment_Point_Posterior_Longitudinal_Inferior_Left,
    ]
    _sub = poi.keys_subregion()
    if any(f.value not in _sub for f in subreg_ids_required_for_ivd_generation):
        poi = calc_poi_from_subreg_vert(
            vert_full,
            subreg_full,
            extend_to=poi,
            subreg_id=subreg_ids_required_for_ivd_generation,
        )

    verts_ids: list[int] = vert_full.unique()
    crop = vert_full.compute_crop(dist=2)
    vert_full_cropped = vert_full.apply_crop(crop)
    subreg_full = subreg_full.apply_crop(crop)
    poi_cropped = poi.apply_crop(crop)

    out = vert_full_cropped.get_array()
    with ProcessPoolExecutor() as executor:
        futures = {
            executor.submit(_process_vertebra, i, verts_ids, vert_full_cropped, subreg_full, poi_cropped, dilate): i
            for i in vert_full.unique()
        }

        for future in tqdm(as_completed(futures), total=len(futures), desc="add artifical IVD to seg"):
            result = future.result()
            if result is not None:
                out_c, slices = result
                out_c = out_c.get_array()
                out[slices][out_c != 0] = out_c[out_c != 0]
    vert_full[crop] = out
    return vert_full

calculate_IVD_POI

calculate_IVD_POI(vert: NII, subreg: NII, poi: POI, ivd_location: set[Location] | None = None) -> POI

Compute IVD-related Points of Interest and add them to poi.

If the subregion image does not yet contain IVD labels (label 100), :func:compute_fake_ivd is called first to synthesise them. Centroid POIs at label 100 are then computed for all vertebrae and optionally extended with superior / inferior extreme disc points via :func:strategy_calculate_up_vector.

Parameters:

Name	Type	Description	Default
vert	NII	Vertebra segmentation NIfTI.	required
subreg	NII	Subregion segmentation NIfTI; updated in place when IVD labels are synthesised.	required
poi	POI	POI object to extend. Updated in place and returned.	required
ivd_location	set[Location] | None	Set of :class:~TPTBox.Location values to compute. Defaults to {Location.Vertebra_Disc_Superior, Location.Vertebra_Disc}. Pass an empty set to skip computation entirely.	None

Returns:

Type	Description
POI	The updated :class:~TPTBox.POI object with IVD-related locations
POI	added.

Source code in TPTBox/spine/spinestats/ivd_pois.py

def calculate_IVD_POI(
    vert: NII,
    subreg: NII,
    poi: POI,
    ivd_location: set[Location] | None = None,
) -> POI:
    """Compute IVD-related Points of Interest and add them to ``poi``.

    If the subregion image does not yet contain IVD labels (label 100),
    :func:`compute_fake_ivd` is called first to synthesise them.  Centroid
    POIs at label 100 are then computed for all vertebrae and optionally
    extended with superior / inferior extreme disc points via
    :func:`strategy_calculate_up_vector`.

    Args:
        vert: Vertebra segmentation NIfTI.
        subreg: Subregion segmentation NIfTI; updated in place when IVD labels
            are synthesised.
        poi: POI object to extend.  Updated in place and returned.
        ivd_location: Set of :class:`~TPTBox.Location` values to compute.
            Defaults to
            ``{Location.Vertebra_Disc_Superior, Location.Vertebra_Disc}``.
            Pass an empty set to skip computation entirely.

    Returns:
        The updated :class:`~TPTBox.POI` object with IVD-related locations
        added.
    """
    if ivd_location is None:
        ivd_location = {Location.Vertebra_Disc_Superior, Location.Vertebra_Disc}
    # Note: currently we always need point 100, so it is computed if any Location is in ivd_location

    if len(ivd_location) == 0:
        return poi
    if 100 not in subreg.unique():
        if Location.Vertebra_Direction_Inferior.value not in poi.keys_subregion():
            poi = calc_poi_from_subreg_vert(
                vert,
                subreg,
                extend_to=poi,
                subreg_id=Location.Vertebra_Direction_Inferior.value,
            )
        vert = compute_fake_ivd(vert, subreg, poi=poi)
        subreg[vert >= 100] = 100

    calc_poi_from_subreg_vert(vert, subreg, subreg_id=100, extend_to=poi)
    if Location.Vertebra_Disc_Superior in ivd_location or Location.Vertebra_Disc_Inferior in ivd_location:
        current_vert = vert.copy()
        bb = current_vert.compute_crop()
        current_vert.apply_crop_(bb)
        for idx in vert.unique():
            if idx >= 100:
                break
            strategy_calculate_up_vector(poi, current_vert, idx, bb)
    return poi