Texture Features API

def calculate_all_texture_features(
    disc_array: npt.NDArray[np.floating[Any]],
    mask_array: npt.NDArray[np.floating[Any]],
    n_bins: int,
    distance_mask_array: Optional[npt.NDArray[np.floating[Any]]] = None,
    ngldm_alpha: int = 0,
    glcm_distance: int = 1,
    ngtdm_distance: int = 1,
    ngldm_distance: int = 1,
    families: Optional[str | Iterable[str]] = None,
) -> dict[str, float]:
    """
    Calculate all texture features (GLCM, GLRLM, GLSZM, GLDZM, NGTDM, NGLDM).

    This is a convenience wrapper that computes the texture matrices once and then
    extracts the features. `families` chooses the families; the matrices of the other
    families are not computed.

    Args:
        disc_array: Discretised image array.
        mask_array: Mask array (ROI). Nonzero values indicate ROI membership.
        n_bins: Number of bins.
        distance_mask_array: Optional mask for GLDZM distance calculation.
                             If None, mask_array is used.
        ngldm_alpha: The coarseness parameter α for NGLDM. Two grey levels are considered
            dependent if their absolute difference is ≤ α. Default is 0 (IBSI standard).
        glcm_distance: The distance δ between the two voxels of a GLCM pair along each of
            the 13 directions (IBSI 1). Default 1.
        ngtdm_distance: The Chebyshev distance δ of the NGTDM neighbourhood. Default 1.
        ngldm_distance: The Chebyshev distance δ of the NGLDM neighbourhood. Default 1.
        families: The families to calculate, from "glcm", "glrlm", "glszm", "gldzm",
            "ngtdm" and "ngldm" (one name, or a list). None (default): all six.

    Returns:
        Dictionary of the texture features of the families, in the order above.

    Raises:
        ValueError: If a family name is not one of the six.
        ValueError: a mask voxel lies outside `distance_mask_array`.

    Example:
        ```python
        import numpy as np
        from pictologics.features.texture import calculate_all_texture_features

        rng = np.random.default_rng(0)
        disc_array = rng.integers(1, 9, size=(10, 10, 10)).astype(np.float64)
        mask_array = np.ones((10, 10, 10), dtype=np.uint8)

        features = calculate_all_texture_features(disc_array, mask_array, n_bins=8)
        print(len(features))
        # 95
        print(round(features["contrast_ACUI"], 3))
        # 10.526
        glcm_only = calculate_all_texture_features(disc_array, mask_array, 8, families="glcm")
        ```
    """
    chosen = set(
        _TEXTURE_FAMILIES
        if families is None
        else [families]
        if isinstance(families, str)
        else families
    )
    for name in chosen - set(_TEXTURE_FAMILIES):
        close = difflib.get_close_matches(str(name), _TEXTURE_FAMILIES, n=1)
        hint = f" (did you mean '{close[0]}'?)" if close else ""
        raise ValueError(
            f"Unknown texture family {name!r}{hint}. The families: {', '.join(_TEXTURE_FAMILIES)}."
        )
    results: dict[str, float] = {}

    # Crop once to the ROI bounding box and use the cropped arrays everywhere below.
    # The per-family feature functions scan the mask (ROI voxel counts, GLCM Ng_eff),
    # so passing full-volume arrays would repeat full-volume scans per family.
    # Feature values are unchanged: cropping only removes zero-mask voxels.
    disc_c, mask_c, distmask_c = _maybe_crop_to_bbox(disc_array, mask_array, distance_mask_array)

    # Calculate all matrices once (compact: same features, smaller matrices)
    planar = _planar_axes(disc_array.shape)
    texture_matrices = _texture_matrices(
        disc_c,
        mask_c,
        n_bins,
        distance_mask=distmask_c,
        ngldm_alpha=ngldm_alpha,
        calc_glcm="glcm" in chosen,
        calc_glrlm="glrlm" in chosen,
        calc_ngtdm="ngtdm" in chosen,
        calc_ngldm="ngldm" in chosen,
        calc_glszm="glszm" in chosen,
        calc_gldzm="gldzm" in chosen,
        compact=True,
        planar=planar,
        glcm_distance=glcm_distance,
        ngtdm_distance=ngtdm_distance,
        ngldm_distance=ngldm_distance,
    )
    # The texture voxel count and the GLCM Ng_eff come with the matrices: the feature
    # functions make no pass over the ROI (`roi` is the bool mask of the texture voxels)
    roi = texture_matrices["roi"]
    n_voxels = texture_matrices["n_voxels"]

    if "glcm" in chosen:
        results.update(_glcm_features(texture_matrices["glcm"], n_bins, texture_matrices["ng_eff"]))
    if "glrlm" in chosen:
        results.update(
            _glrlm_features(
                texture_matrices["glrlm"], roi, _directions(planar).size, n_voxels=n_voxels
            )
        )
    if "glszm" in chosen:
        results.update(
            _glszm_features_from_cells(texture_matrices["glszm_cells"], roi, n_voxels=n_voxels)
        )
    if "gldzm" in chosen:
        results.update(_gldzm_features(texture_matrices["gldzm"], roi, n_voxels=n_voxels))
    if "ngtdm" in chosen:
        results.update(
            calculate_ngtdm_features(
                disc_c,
                mask_c,
                n_bins,
                ngtdm_matrices=(texture_matrices["ngtdm_s"], texture_matrices["ngtdm_n"]),
            )
        )
    if "ngldm" in chosen:
        results.update(_ngldm_features(texture_matrices["ngldm"], roi, n_voxels=n_voxels))
    return results

def calculate_all_texture_matrices(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    distance_mask: Optional[npt.NDArray[np.floating[Any]]] = None,
    ngldm_alpha: int = 0,
    calc_glcm: bool = True,
    calc_glrlm: bool = True,
    calc_ngtdm: bool = True,
    calc_ngldm: bool = True,
    calc_glszm: bool = True,
    calc_gldzm: bool = True,
    glcm_distance: int = 1,
    ngtdm_distance: int = 1,
    ngldm_distance: int = 1,
) -> dict[str, Any]:
    """
    Calculate all texture matrices (GLCM, GLRLM, GLSZM, GLDZM, NGTDM, NGLDM) in an optimized single pass.

    This function serves as the computational backbone for texture analysis. It computes the raw
    matrices required to extract specific texture features. By aggregating these calculations,
    it minimizes the number of passes over the image data, significantly improving performance.

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
            Values should be integers in the range [1, n_bins].
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the Region of Interest (ROI).
            Must have the same shape as `data`. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels used for discretization (e.g., 16, 32, 64).
            This determines the size of the resulting matrices.
        distance_mask (Optional[npt.NDArray[np.floating[Any]]]): Optional mask used to calculate the distance map for GLDZM.
            If None, `mask` is used. This allows calculating distances based on the morphological mask
            while analyzing intensities from the intensity mask (e.g., after outlier filtering).
        ngldm_alpha (int): The coarseness parameter α for NGLDM calculation. Two grey levels are
            considered dependent if their absolute difference is ≤ α. Default is 0 (exact match),
            which is the IBSI standard. Use α=1 for tolerance of ±1 grey level difference.
        calc_glcm (bool): Whether to compute the GLCM. Default True.
        calc_glrlm (bool): Whether to compute the GLRLM. Default True.
        calc_ngtdm (bool): Whether to compute the NGTDM. Default True.
        calc_ngldm (bool): Whether to compute the NGLDM. Default True.
        calc_glszm (bool): Whether to compute the GLSZM. Default True.
        calc_gldzm (bool): Whether to compute the GLDZM (skipping it also skips the
            distance transform). Default True.
        glcm_distance (int): The distance δ between the two voxels of a GLCM pair along
            each of the 13 directions (IBSI 1). Default 1.
        ngtdm_distance (int): The Chebyshev distance δ of the NGTDM neighbourhood. Default 1
            (the 26 neighbours).
        ngldm_distance (int): The Chebyshev distance δ of the NGLDM neighbourhood. Default 1.

    Returns:
        dict[str, Any]: A dictionary containing the calculated texture matrices (a disabled
            matrix is a zero-filled placeholder of minimal shape):
            - 'glcm' (npt.NDArray[np.floating[Any]]): Grey Level Co-occurrence Matrix. Shape: (n_dirs, n_bins, n_bins).
            - 'glrlm' (npt.NDArray[np.floating[Any]]): Grey Level Run Length Matrix. Shape: (n_dirs, n_bins, max_run_length).
            - 'ngtdm_s' (npt.NDArray[np.floating[Any]]): NGTDM Sum of absolute differences. Shape: (n_bins,).
            - 'ngtdm_n' (npt.NDArray[np.floating[Any]]): NGTDM Number of valid voxels. Shape: (n_bins,).
            - 'ngldm' (npt.NDArray[np.floating[Any]]): Neighbouring Grey Level Dependence Matrix. Shape: (n_bins, n_dependence).
            - 'glszm' (npt.NDArray[np.floating[Any]]): Grey Level Size Zone Matrix. Shape: (n_bins, max_zone_size).
            - 'gldzm' (npt.NDArray[np.floating[Any]]): Grey Level Distance Zone Matrix. Shape: (n_bins, max_distance).

    Raises:
        ValueError: a mask voxel lies outside `distance_mask`.

    Example:
        Calculate all texture matrices:

        ```python
        import numpy as np
        from pictologics.features.texture import calculate_all_texture_matrices

        # Create dummy data
        data = np.random.randint(1, 33, (50, 50, 50))
        mask = np.ones((50, 50, 50))

        # Calculate matrices
        matrices = calculate_all_texture_matrices(data, mask, n_bins=32)
        print(matrices['glcm'].shape)
        # (13, 32, 32)
        ```
    """
    return _texture_matrices(
        data,
        mask,
        n_bins,
        distance_mask,
        ngldm_alpha,
        calc_glcm,
        calc_glrlm,
        calc_ngtdm,
        calc_ngldm,
        calc_glszm,
        calc_gldzm,
        compact=False,
        glcm_distance=glcm_distance,
        ngtdm_distance=ngtdm_distance,
        ngldm_distance=ngldm_distance,
    )

def calculate_glcm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    glcm_matrix: Optional[npt.NDArray[np.floating[Any]]] = None,
    glcm_distance: int = 1,
) -> dict[str, float]:
    r"""
    Calculate Grey Level Co-occurrence Matrix (GLCM) features.

    The GLCM describes the second-order statistical distribution of grey levels in the ROI.
    It counts how often pairs of grey levels occur at a specific distance and direction.
    This implementation computes features based on the 3D merged GLCM (averaged over all 13 directions),
    making the features rotationally invariant.

    **IBSI Reference**: Section 3.6 (Grey Level Co-occurrence Based Features).

    **Mathematical Definition**:
    Let $P(i,j)$ be the co-occurrence matrix, where $i$ and $j$ are grey levels.
    The matrix is normalized such that $\sum_{i,j} P(i,j) = 1$.

    **Calculated Features**:
    *   Joint Maximum (GYBY)
    *   Joint Average (60VM)
    *   Joint Variance (UR99)
    *   Joint Entropy (TU9B)
    *   Difference Average (TF7R)
    *   Difference Variance (D3YU)
    *   Difference Entropy (NTRS)
    *   Sum Average (ZGXS)
    *   Sum Variance (OEEB)
    *   Sum Entropy (P6QZ)
    *   Angular Second Moment (8ZQL)
    *   Contrast (ACUI)
    *   Dissimilarity (8S9J)
    *   Inverse Difference (IB1Z)
    *   Normalised Inverse Difference (NDRX)
    *   Inverse Difference Moment (WF0Z)
    *   Normalised Inverse Difference Moment (1QCO)
    *   Inverse Variance (E8JP)
    *   Correlation (NI2N)
    *   Autocorrelation (QWB0)
    *   Cluster Tendency (DG8W)
    *   Cluster Shade (7NFM)
    *   Cluster Prominence (AE86)
    *   Information Correlation 1 (R8DG)
    *   Information Correlation 2 (JN9H)

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        glcm_matrix (Optional[npt.NDArray[np.floating[Any]]]): Pre-calculated GLCM matrix. If provided, `data` and `mask`
            are ignored for matrix calculation, but `data` is still used for `Ng` estimation if needed.
            If None, the matrix is calculated from scratch.
        glcm_distance (int): The distance δ between the two voxels of a pair along each
            of the 13 directions (IBSI 1), for a matrix calculated from scratch. Default 1.

    Returns:
        dict[str, float]: A dictionary of calculated GLCM features, keyed by their name and IBSI code.
            Example keys: 'joint_maximum_GYBY', 'contrast_ACUI', 'correlation_NI2N'.

    Example:
        ```python
        import numpy as np
        from pictologics.features.texture import calculate_glcm_features

        rng = np.random.default_rng(0)
        data = rng.integers(1, 33, (20, 20, 20)).astype(np.float64)
        mask = np.ones((20, 20, 20), dtype=np.uint8)
        features = calculate_glcm_features(data, mask, n_bins=32)
        print(round(features["contrast_ACUI"], 1))
        ```
    """
    if glcm_matrix is None:
        # Standalone path: crop to the ROI bbox, so the kernel runs on the cropped
        # arrays; Ng_eff comes with the matrices.
        planar = _planar_axes(data.shape)
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        matrices = _one_matrix(data, mask, n_bins, "glcm", planar=planar, distance=glcm_distance)
        return _glcm_features(matrices["glcm"], n_bins, matrices["ng_eff"])
    glcm = glcm_matrix
    if glcm.dtype.kind == "u":  # counts: the kernel
        if not glcm.any():
            return {}
        return _glcm_features(glcm, n_bins, _ng_eff(data, mask))

    # Merge (Sum) -> IAZD
    glcm_sum = np.sum(glcm, axis=0)
    glcm_sym = glcm_sum + glcm_sum.T

    # Normalize
    total_sum = np.sum(glcm_sym)
    if total_sum == 0:
        return {}
    glcm_sym, levels, _ = _occupied(glcm_sym, symmetric=True)

    P = glcm_sym / total_sum

    # Ng_eff is the ROI grey-level span
    Ng_eff = _ng_eff(data, mask)

    # The sums over the cells that hold counts, in two serial passes (_glcm_sums_numba,
    # _glcm_mu_sums_numba); levels: the 1-based grey levels of the rows and columns of P
    levels = levels.astype(np.int64, copy=False)
    # The bins hold every level of P (a given matrix can hold more levels than n_bins)
    n_bins = max(n_bins, int(levels[-1]))
    p_diff = np.zeros(n_bins)
    p_sum_full = np.zeros(2 * n_bins + 1)
    (
        joint_max,
        joint_entropy,
        difference_average,
        contrast,
        angular_second_moment,
        inverse_difference,
        normalised_inverse_difference,
        inverse_difference_moment,
        normalised_inverse_difference_moment,
        inverse_variance,
        autocorrelation,
    ) = _glcm_sums_numba(P, levels, float(Ng_eff), p_diff, p_sum_full)
    # The joint average from the integer counts: exact to half a unit in the last place, as
    # the cluster shade (whose terms cancel) is sensitive to it
    row_counts = glcm_sym.sum(axis=1)
    joint_average = int(levels @ row_counts.astype(np.int64)) / int(total_sum)
    mu = joint_average
    joint_variance, correlation_term, cluster_tendency, cluster_shade, cluster_prominence = (
        _glcm_mu_sums_numba(P, levels, mu)
    )

    features = {}
    features["joint_maximum_GYBY"] = joint_max
    features["joint_average_60VM"] = joint_average
    features["joint_variance_UR99"] = joint_variance
    features["joint_entropy_TU9B"] = joint_entropy
    features["difference_average_TF7R"] = difference_average

    # Difference Variance - D3YU and Difference Entropy - NTRS, from the difference bins
    dk = np.arange(n_bins) - difference_average
    features["difference_variance_D3YU"] = np.sum(dk * dk * p_diff)
    mask_pd = p_diff > 0
    features["difference_entropy_NTRS"] = -np.sum(p_diff[mask_pd] * np.log2(p_diff[mask_pd]))

    # Sum Average - ZGXS, Sum Variance - OEEB and Sum Entropy - P6QZ, from the sum bins
    p_sum = p_sum_full[2:]
    k_vals_sum = np.arange(2, 2 * n_bins + 1)
    features["sum_average_ZGXS"] = np.sum(k_vals_sum * p_sum)
    dk = k_vals_sum - features["sum_average_ZGXS"]
    features["sum_variance_OEEB"] = np.sum(dk * dk * p_sum)
    mask_ps = p_sum > 0
    features["sum_entropy_P6QZ"] = -np.sum(p_sum[mask_ps] * np.log2(p_sum[mask_ps]))

    features["angular_second_moment_8ZQL"] = angular_second_moment
    features["contrast_ACUI"] = contrast
    features["dissimilarity_8S9J"] = difference_average
    features["inverse_difference_IB1Z"] = inverse_difference
    features["normalised_inverse_difference_NDRX"] = normalised_inverse_difference
    features["inverse_difference_moment_WF0Z"] = inverse_difference_moment
    features["normalised_inverse_difference_moment_1QCO"] = normalised_inverse_difference_moment
    features["inverse_variance_E8JP"] = inverse_variance
    if joint_variance != 0:
        features["correlation_NI2N"] = correlation_term / joint_variance
    else:
        features["correlation_NI2N"] = 1.0  # IBSI does not define the value for 0 variance
    features["autocorrelation_QWB0"] = autocorrelation
    features["cluster_tendency_DG8W"] = cluster_tendency
    features["cluster_shade_7NFM"] = cluster_shade
    features["cluster_prominence_AE86"] = cluster_prominence

    # Information Correlation 1 - R8DG and 2 - JN9H
    HXY = joint_entropy
    # The marginal from the counts: exact, so it does not depend on the order of the sums
    p_x = row_counts / total_sum
    mask_px = p_x > 0
    HX = -np.sum(p_x[mask_px] * np.log2(p_x[mask_px]))
    # P is symmetric, so p_y = p_x, and Σ_j P_ij = p_i. The two entropies of the IBSI
    # definitions have closed forms (the sums run over the cells that hold counts):
    #   HXY1 = -Σ_ij P_ij log2(p_i p_j) = -Σ_i p_i log2 p_i - Σ_j p_j log2 p_j = 2 HX
    #   HXY2 = -Σ_ij p_i p_j log2(p_i p_j) = 2 HX Σ_j p_j = 2 HX
    if HX != 0:
        features["information_correlation_1_R8DG"] = (HXY - 2 * HX) / HX
    else:
        features["information_correlation_1_R8DG"] = np.nan
    # HXY2 - HXY is 0 or more; rounding can make it slightly negative (P = p p^T)
    features["information_correlation_2_JN9H"] = np.sqrt(1 - np.exp(-2 * max(0.0, 2 * HX - HXY)))
    return features

def calculate_gldzm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    gldzm_matrix: Optional[npt.NDArray[np.floating[Any]]] = None,
    distance_mask: Optional[npt.NDArray[np.floating[Any]]] = None,
) -> dict[str, float]:
    """
    Calculate Grey Level Distance Zone Matrix (GLDZM) features.

    The GLDZM counts the number of zones of linked voxels with the same grey level,
    categorized by the distance of the zone from the ROI border.
    This captures information about the spatial distribution of textures relative to the boundary.

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        gldzm_matrix (Optional[npt.NDArray[np.floating[Any]]]): Pre-calculated GLDZM matrix.
        distance_mask (Optional[npt.NDArray[np.floating[Any]]]): Optional mask used to calculate the distance map.
            If None, `mask` is used. This allows calculating distances based on the morphological mask
            while analyzing intensities from the intensity mask (e.g., after outlier filtering).

    Returns:
        dict[str, float]: A dictionary of calculated GLDZM features.
            Example keys: 'small_distance_emphasis_0GBI', 'zone_distance_entropy_GBDU'.

    Raises:
        ValueError: a mask voxel lies outside `distance_mask`.
    """
    if gldzm_matrix is None:
        # Standalone path: crop to the bbox of mask ∪ distance_mask. The union keeps the
        # distance map identical (the distance region is fully inside the bbox); the
        # voxel count of the zone percentage comes with the matrices.
        planar = _planar_axes(data.shape)
        data, mask, distance_mask = _maybe_crop_to_bbox(data, mask, distance_mask)
        matrices = _one_matrix(
            data, mask, n_bins, "gldzm", distance_mask=distance_mask, planar=planar
        )
        return _gldzm_features(matrices["gldzm"], None, n_voxels=matrices["n_voxels"])
    gldzm = gldzm_matrix
    if gldzm.dtype.kind == "u":  # counts: the kernel
        return _gldzm_features(gldzm, mask)

    N_zones = np.sum(gldzm)
    if N_zones == 0:
        return {}

    gldzm, rows, cols = _occupied(gldzm)
    P = gldzm / N_zones

    # The 1-based numbers of the rows and columns that hold counts
    I, J = np.meshgrid(rows, cols, indexing="ij")  # noqa: E741
    I2 = I**2
    J2 = J**2

    features = {}

    # Small Distance Emphasis (SDE) - 0GBI
    features["small_distance_emphasis_0GBI"] = np.sum(P / J2)

    # Large Distance Emphasis (LDE) - MB4I
    features["large_distance_emphasis_MB4I"] = np.sum(P * J2)

    # Grey Level Non-Uniformity (GLNU) - VFT7
    s_i = np.sum(gldzm, axis=1)
    features["grey_level_non_uniformity_VFT7"] = np.sum(s_i**2) / N_zones

    # Normalised Grey Level Non-Uniformity (GLNN) - 7HP3
    features["normalised_grey_level_non_uniformity_7HP3"] = np.sum(s_i**2) / (N_zones**2)

    # Zone Distance Non-Uniformity (ZDNU) - V294
    s_j = np.sum(gldzm, axis=0)
    features["zone_distance_non_uniformity_V294"] = np.sum(s_j**2) / N_zones

    # Normalised Zone Distance Non-Uniformity (ZDNN) - IATH
    features["normalised_zone_distance_non_uniformity_IATH"] = np.sum(s_j**2) / (N_zones**2)

    # Zone Percentage (ZP) - VIWW
    n_voxels = _roi_voxel_count(mask)
    features["zone_percentage_VIWW"] = N_zones / n_voxels

    # Grey Level Variance (GLV) - QK93
    mu_i = np.sum(I * P)
    features["grey_level_variance_QK93"] = np.sum(((I - mu_i) ** 2) * P)

    # Zone Distance Variance (ZDV) - 7WT1
    mu_j = np.sum(J * P)
    features["zone_distance_variance_7WT1"] = np.sum(((J - mu_j) ** 2) * P)

    # Zone Distance Entropy (ZDE) - GBDU
    mask_p = P > 0
    features["zone_distance_entropy_GBDU"] = -np.sum(P[mask_p] * np.log2(P[mask_p]))

    # Low Grey Level Zone Emphasis (LGLZE) - S1RA
    features["low_grey_level_zone_emphasis_S1RA"] = np.sum(P / I2)

    # High Grey Level Zone Emphasis (HGLZE) - K26C
    features["high_grey_level_zone_emphasis_K26C"] = np.sum(P * I2)

    # Small Distance Low Grey Level Emphasis (SDLGLE) - RUVG
    features["small_distance_low_grey_level_emphasis_RUVG"] = np.sum(P / (I2 * J2))

    # Small Distance High Grey Level Emphasis (SDHGLE) - DKNJ
    features["small_distance_high_grey_level_emphasis_DKNJ"] = np.sum(P * I2 / J2)

    # Large Distance Low Grey Level Emphasis (LDLGLE) - A7WM
    features["large_distance_low_grey_level_emphasis_A7WM"] = np.sum(P * J2 / I2)

    # Large Distance High Grey Level Emphasis (LDHGLE) - KLTH
    features["large_distance_high_grey_level_emphasis_KLTH"] = np.sum(P * I2 * J2)

    return features

def calculate_glrlm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    glrlm_matrix: Optional[npt.NDArray[np.floating[Any]]] = None,
    n_directions: Optional[int] = None,
) -> dict[str, float]:
    """
    Calculate Grey Level Run Length Matrix (GLRLM) features.

    The GLRLM quantifies grey level runs, which are defined as the length in number of pixels,
    of consecutive pixels that have the same grey level value.
    This implementation computes features based on the 3D merged GLRLM (averaged over all 13 directions).

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        glrlm_matrix (Optional[npt.NDArray[np.floating[Any]]]): Pre-calculated GLRLM matrix.
        n_directions (Optional[int]): The number of directions in a merged (one-table)
            `glrlm_matrix`, for the run percentage: 13 in 3D, 4 for a one-slice image.
            None: 13, or for 13 direction tables the number of them that hold runs.

    Returns:
        dict[str, float]: A dictionary of calculated GLRLM features.
            Example keys: 'short_runs_emphasis_22OV', 'grey_level_non_uniformity_R5YN'.
    """
    if glrlm_matrix is None:
        # Standalone path: crop to the ROI bbox, so the kernel runs on the cropped
        # arrays; the voxel count of the run percentage comes with the matrices.
        planar = _planar_axes(data.shape)
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        matrices = _one_matrix(data, mask, n_bins, "glrlm", planar=planar)
        return _glrlm_features(
            matrices["glrlm"], None, _directions(planar).size, n_voxels=matrices["n_voxels"]
        )
    glrlm = glrlm_matrix
    if n_directions is None and glrlm.shape[0] == 13:
        # The directions of a one-slice image that leave it have no runs at all
        n_directions = int(np.count_nonzero(glrlm.reshape(13, -1).any(axis=1)))
    if glrlm.dtype.kind == "u":  # counts: the kernel
        return _glrlm_features(glrlm, mask, n_directions)

    # Merge (Sum) -> IAZD
    glrlm_sum = np.sum(glrlm, axis=0)

    # Remove length 0 (column 0)
    glrlm = glrlm_sum[:, 1:]

    N_runs = np.sum(glrlm)
    if N_runs == 0:
        return {}

    glrlm, rows, cols = _occupied(glrlm)
    P = glrlm / N_runs

    # The 1-based numbers of the rows and columns that hold counts
    I, J = np.meshgrid(rows, cols, indexing="ij")  # noqa: E741
    I2 = I**2
    J2 = J**2

    features = {}

    # Short Run Emphasis (SRE) - 22OV
    features["short_runs_emphasis_22OV"] = np.sum(P / J2)

    # Long Run Emphasis (LRE) - W4KF
    features["long_runs_emphasis_W4KF"] = np.sum(P * J2)

    # Grey Level Non-Uniformity (GLNU) - R5YN
    s_i = np.sum(glrlm, axis=1)
    features["grey_level_non_uniformity_R5YN"] = np.sum(s_i**2) / N_runs

    # Normalised Grey Level Non-Uniformity (GLNN) - OVBL
    features["normalised_grey_level_non_uniformity_OVBL"] = np.sum(s_i**2) / (N_runs**2)

    # Run Length Non-Uniformity (RLNU) - W92Y
    s_j = np.sum(glrlm, axis=0)
    features["run_length_non_uniformity_W92Y"] = np.sum(s_j**2) / N_runs

    # Normalised Run Length Non-Uniformity (RLNN) - IC23
    features["normalised_run_length_non_uniformity_IC23"] = np.sum(s_j**2) / (N_runs**2)

    # Run Percentage (RP) - 9ZK5
    n_voxels = _roi_voxel_count(mask)
    features["run_percentage_9ZK5"] = N_runs / (n_voxels * (n_directions or 13))

    # Grey Level Variance (GLV) - 8CE5
    mu_i = np.sum(I * P)
    features["grey_level_variance_8CE5"] = np.sum(((I - mu_i) ** 2) * P)

    # Run Length Variance (RLV) - SXLW
    mu_j = np.sum(J * P)
    features["run_length_variance_SXLW"] = np.sum(((J - mu_j) ** 2) * P)

    # Run Entropy (RE) - HJ9O
    mask_p = P > 0
    features["run_entropy_HJ9O"] = -np.sum(P[mask_p] * np.log2(P[mask_p]))

    # Low Grey Level Run Emphasis (LGLRE) - V3SW
    features["low_grey_level_run_emphasis_V3SW"] = np.sum(P / I2)

    # High Grey Level Run Emphasis (HGLRE) - G3QZ
    features["high_grey_level_run_emphasis_G3QZ"] = np.sum(P * I2)

    # Short Run Low Grey Level Emphasis (SRLGLE) - HTZT
    features["short_run_low_grey_level_emphasis_HTZT"] = np.sum(P / (I2 * J2))

    # Short Run High Grey Level Emphasis (SRHGLE) - GD3A
    features["short_run_high_grey_level_emphasis_GD3A"] = np.sum(P * I2 / J2)

    # Long Run Low Grey Level Emphasis (LRLGLE) - IVPO
    features["long_run_low_grey_level_emphasis_IVPO"] = np.sum(P * J2 / I2)

    # Long Run High Grey Level Emphasis (LRHGLE) - 3KUM
    features["long_run_high_grey_level_emphasis_3KUM"] = np.sum(P * I2 * J2)

    return features

def calculate_glszm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    glszm_matrix: Optional[npt.NDArray[np.floating[Any]]] = None,
) -> dict[str, float]:
    """
    Calculate Grey Level Size Zone Matrix (GLSZM) features.

    The GLSZM counts the number of zones (connected components) of linked voxels
    that share the same grey level intensity. A zone is defined as a group of connected voxels
    with the same grey level. This matrix is rotationally invariant by definition.

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        glszm_matrix (Optional[npt.NDArray[np.floating[Any]]]): Pre-calculated GLSZM matrix.

    Returns:
        dict[str, float]: A dictionary of calculated GLSZM features.
            Example keys: 'small_zone_emphasis_P001', 'zone_percentage_P30P'.
    """
    if glszm_matrix is None:
        # Standalone path: crop to the ROI bbox, so the zone kernel runs on the cropped
        # arrays; the voxel count of the zone percentage comes with the matrices.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        # The kernel returns only the non-zero cells: no dense (n_bins, largest zone) matrix.
        matrices = _one_matrix(data, mask, n_bins, "glszm")
        return _glszm_features_from_cells(
            matrices["glszm_cells"], None, n_voxels=matrices["n_voxels"]
        )

    # np.nonzero gives strided index arrays: row-order copies, for the compiled kernel
    gl_idx, sz_idx = (np.ascontiguousarray(i) for i in np.nonzero(glszm_matrix))
    return _glszm_features(gl_idx, sz_idx, glszm_matrix[gl_idx, sz_idx].astype(np.float64), mask)

def calculate_ngldm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    ngldm_matrix: Optional[npt.NDArray[np.floating[Any]]] = None,
    ngldm_alpha: int = 0,
    ngldm_distance: int = 1,
) -> dict[str, float]:
    """
    Calculate Neighbourhood Grey Level Dependence Matrix (NGLDM) features.

    The NGLDM captures the dependence of grey levels on their neighbours.
    A "dependence" is defined as a connected voxel having a similar grey level (within a tolerance α).

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        ngldm_matrix (Optional[npt.NDArray[np.floating[Any]]]): Pre-calculated NGLDM matrix.
        ngldm_alpha (int): The coarseness parameter α. Two grey levels are considered dependent
            if their absolute difference is ≤ α. Default is 0 (exact match, IBSI standard).
        ngldm_distance (int): The Chebyshev distance δ of the neighbourhood, for a matrix
            calculated from scratch. Default 1.

    Returns:
        dict[str, float]: A dictionary of calculated NGLDM features.
            Example keys: 'low_dependence_emphasis_SODN', 'dependence_count_entropy_FCBV'.
    """
    if ngldm_matrix is None:
        # Standalone path: crop to the ROI bbox, so the kernel runs on the cropped
        # arrays; the voxel count of the dependence count percentage comes with the matrices.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        matrices = _one_matrix(
            data, mask, n_bins, "ngldm", ngldm_alpha=ngldm_alpha, distance=ngldm_distance
        )
        return _ngldm_features(matrices["ngldm"], None, n_voxels=matrices["n_voxels"])
    ngldm = ngldm_matrix
    if ngldm.dtype.kind == "u":  # counts: the kernel
        return _ngldm_features(ngldm, mask)

    N_s = np.sum(ngldm)
    if N_s == 0:
        return {}

    ngldm, rows, cols = _occupied(ngldm)
    P = ngldm / N_s

    # The 1-based numbers of the rows and columns that hold counts
    I, J = np.meshgrid(rows, cols, indexing="ij")  # noqa: E741
    I2 = I**2
    J2 = J**2

    features = {}

    # Low Dependence Emphasis (LDE) - SODN
    features["low_dependence_emphasis_SODN"] = np.sum(P / J2)

    # High Dependence Emphasis (HDE) - IMOQ
    features["high_dependence_emphasis_IMOQ"] = np.sum(P * J2)

    # Low Grey Level Count Emphasis (LGCE) - TL9H
    features["low_grey_level_count_emphasis_TL9H"] = np.sum(P / I2)

    # High Grey Level Count Emphasis (HGCE) - OAE7
    features["high_grey_level_count_emphasis_OAE7"] = np.sum(P * I2)

    # Low Dependence Low Grey Level Emphasis (LDLGE) - EQ3F
    features["low_dependence_low_grey_level_emphasis_EQ3F"] = np.sum(P / (I2 * J2))

    # Low Dependence High Grey Level Emphasis (LDHGE) - JA6D
    features["low_dependence_high_grey_level_emphasis_JA6D"] = np.sum(P * I2 / J2)

    # High Dependence Low Grey Level Emphasis (HDLGE) - NBZI
    features["high_dependence_low_grey_level_emphasis_NBZI"] = np.sum(P * J2 / I2)

    # High Dependence High Grey Level Emphasis (HDHGE) - 9QMG
    features["high_dependence_high_grey_level_emphasis_9QMG"] = np.sum(P * I2 * J2)

    # Grey Level Non-Uniformity - FP8K
    s_i = np.sum(ngldm, axis=1)
    features["grey_level_non_uniformity_FP8K"] = np.sum(s_i**2) / N_s

    # Normalised Grey Level Non-Uniformity - 5SPA
    features["normalised_grey_level_non_uniformity_5SPA"] = np.sum(s_i**2) / (N_s**2)

    # Dependence Count Non-Uniformity - Z87G
    s_j = np.sum(ngldm, axis=0)
    features["dependence_count_non_uniformity_Z87G"] = np.sum(s_j**2) / N_s

    # Normalised Dependence Count Non-Uniformity - OKJI
    features["normalised_dependence_count_non_uniformity_OKJI"] = np.sum(s_j**2) / (N_s**2)

    # Dependence Count Percentage - 6XV8
    n_voxels = _roi_voxel_count(mask)
    features["dependence_count_percentage_6XV8"] = N_s / n_voxels

    # Grey Level Variance - 1PFV
    mu_i = np.sum(I * P)
    features["grey_level_variance_1PFV"] = np.sum(((I - mu_i) ** 2) * P)

    # Dependence Count Variance - DNX2
    mu_j = np.sum(J * P)
    features["dependence_count_variance_DNX2"] = np.sum(((J - mu_j) ** 2) * P)

    # Dependence Count Entropy - FCBV
    mask_p = P > 0
    features["dependence_count_entropy_FCBV"] = -np.sum(P[mask_p] * np.log2(P[mask_p]))

    # Dependence Count Energy - CAS9
    features["dependence_count_energy_CAS9"] = np.sum(P**2)

    return features

def calculate_ngtdm_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    n_bins: int,
    ngtdm_matrices: Optional[
        tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]
    ] = None,
    ngtdm_distance: int = 1,
) -> dict[str, float]:
    """
    Calculate Neighbourhood Grey Tone Difference Matrix (NGTDM) features.

    The NGTDM quantifies the difference between a grey value and the average grey value
    of its neighbours. It captures the coarseness and contrast of the texture.

    Args:
        data (npt.NDArray[np.floating[Any]]): The 3D image array containing discretised grey levels.
        mask (npt.NDArray[np.floating[Any]]): The 3D mask array defining the ROI. Nonzero values indicate ROI membership.
        n_bins (int): The number of grey levels.
        ngtdm_matrices (Optional[tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]]): Pre-calculated NGTDM matrices
            (sum of absolute differences `s`, and count `n`).
        ngtdm_distance (int): The Chebyshev distance δ of the neighbourhood, for matrices
            calculated from scratch. Default 1 (the 26 neighbours).

    Returns:
        dict[str, float]: A dictionary of calculated NGTDM features.
            Example keys: 'coarseness_QCDE', 'contrast_65HE', 'busyness_NQ30'.
    """
    if ngtdm_matrices is None:
        # Standalone path: crop to the ROI bbox so the kernel runs on the cropped arrays.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        matrices = _one_matrix(data, mask, n_bins, "ngtdm", distance=ngtdm_distance)
        s, n = matrices["ngtdm_s"], matrices["ngtdm_n"]
    else:
        s, n = ngtdm_matrices

    # s[i] and n[i]: the sum of the absolute differences and the count of the voxels of
    # grey level i + 1 (with a valid neighbourhood). The kernel computes the features with
    # the bits of the numpy expressions; the closed forms of the sums over the pairs of
    # levels are in its comment.
    out = np.empty(5)
    n_vp = _ngtdm_features_numba(
        np.ascontiguousarray(s, dtype=np.float64), np.ascontiguousarray(n, dtype=np.float64), out
    )
    if n_vp == 0.0:
        return {}
    return {
        "coarseness_QCDE": float(out[0]),
        "contrast_65HE": float(out[1]),
        "busyness_NQ30": float(out[2]),
        "complexity_HDEZ": float(out[3]),
        "strength_1X9X": float(out[4]),
    }