Texture Features API

pictologics.features.texture

Texture Feature Extraction Module

This module provides a comprehensive suite of functions for calculating 3D texture features from medical images. It implements the Image Biomarker Standardisation Initiative (IBSI) compliant algorithms for various texture matrices.

Key Concepts:

Texture analysis quantifies the spatial arrangement of grey levels in an image. It assumes that the texture (e.g., "smooth", "coarse", "regular") is contained in the spatial relationship between the grey levels of the voxels.

Implemented Matrices:

1. GLCM (Grey Level Co-occurrence Matrix): Counts how often pairs of grey levels occur at a specific distance and direction. Captures: Contrast, homogeneity, correlation.

2. GLRLM (Grey Level Run Length Matrix): Counts the lengths of consecutive runs of the same grey level. Captures: Coarseness, directionality.

3. GLSZM (Grey Level Size Zone Matrix): Counts the size of zones (connected components) of the same grey level. Captures: Regional homogeneity, size distribution of texture elements.

4. GLDZM (Grey Level Distance Zone Matrix): Counts zones based on their distance from the ROI border. Captures: Spatial distribution relative to the boundary.

5. NGTDM (Neighbourhood Grey Tone Difference Matrix): Quantifies the difference between a voxel and its neighbours. Captures: Human perception of texture (coarseness, contrast, busyness).

6. NGLDM (Neighbourhood Grey Level Dependence Matrix): Captures the dependence of grey levels on their neighbours. Captures: Dependence, spatial relationships.

Optimization:

This module uses numba for Just-In-Time (JIT) compilation to achieve high performance. The core calculations are parallelized and optimized for memory usage. - Single-pass calculation: Multiple matrices are computed in a single pass over the image to minimize memory access overhead. - Flattened DFS: Zone-based features (GLSZM, GLDZM) use a memory-efficient Depth-First Search with flattened stack indices.

Usage:

The main entry point is calculate_all_texture_matrices, which computes all raw matrices. Then, specific feature calculation functions (e.g., calculate_glcm_features) can be called using these matrices.

Example

Calculate texture features:

import numpy as np
from pictologics.features.texture import (
    calculate_all_texture_matrices,
    calculate_glcm_features
)

# 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)

# Extract features
glcm_feats = calculate_glcm_features(
    data,
    mask,
    n_bins=32,
    glcm_matrix=matrices['glcm']
)
print(glcm_feats['contrast_ACUI'])

calculate_all_texture_features(disc_array, mask_array, n_bins, distance_mask_array=None, ngldm_alpha=0)

Calculate all texture features (GLCM, GLRLM, GLSZM, GLDZM, NGTDM, NGLDM).

This is a convenience wrapper that computes all texture matrices and then extracts all available features.

Parameters:

Name	Type	Description	Default
disc_array	NDArray[floating[Any]]	Discretised image array.	required
mask_array	NDArray[floating[Any]]	Mask array (ROI). Nonzero values indicate ROI membership.	required
n_bins	int	Number of bins.	required
distance_mask_array	Optional[NDArray[floating[Any]]]	Optional mask for GLDZM distance calculation. If None, mask_array is used.	None
ngldm_alpha	int	The coarseness parameter α for NGLDM. Two grey levels are considered dependent if their absolute difference is ≤ α. Default is 0 (IBSI standard).	0

Returns:

Type	Description
dict[str, float]	Dictionary of all texture features.

Example

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

Source code in pictologics/features/texture.py

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,
) -> dict[str, float]:
    """
    Calculate all texture features (GLCM, GLRLM, GLSZM, GLDZM, NGTDM, NGLDM).

    This is a convenience wrapper that computes all texture matrices and then
    extracts all available features.

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

    Returns:
        Dictionary of all texture features.

    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
        ```
    """
    results = {}

    # 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
    texture_matrices = calculate_all_texture_matrices(
        disc_c,
        mask_c,
        n_bins,
        distance_mask=distmask_c,
        ngldm_alpha=ngldm_alpha,
    )

    # GLCM
    results.update(
        calculate_glcm_features(disc_c, mask_c, n_bins, glcm_matrix=texture_matrices["glcm"])
    )

    # GLRLM
    results.update(
        calculate_glrlm_features(disc_c, mask_c, n_bins, glrlm_matrix=texture_matrices["glrlm"])
    )

    # GLSZM
    results.update(
        calculate_glszm_features(disc_c, mask_c, n_bins, glszm_matrix=texture_matrices["glszm"])
    )

    # GLDZM
    results.update(
        calculate_gldzm_features(
            disc_c,
            mask_c,
            n_bins,
            gldzm_matrix=texture_matrices["gldzm"],
            distance_mask=(distmask_c if distmask_c is not None else mask_c),
        )
    )

    # NGTDM
    results.update(
        calculate_ngtdm_features(
            disc_c,
            mask_c,
            n_bins,
            ngtdm_matrices=(texture_matrices["ngtdm_s"], texture_matrices["ngtdm_n"]),
        )
    )
    # NGLDM
    results.update(
        calculate_ngldm_features(disc_c, mask_c, n_bins, ngldm_matrix=texture_matrices["ngldm"])
    )

    return results

calculate_all_texture_matrices(data, mask, n_bins, distance_mask=None, ngldm_alpha=0, calc_glcm=True, calc_glrlm=True, calc_ngtdm=True, calc_ngldm=True, calc_glszm=True, calc_gldzm=True)

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.

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels. Values should be integers in the range [1, n_bins].	required
mask	NDArray[floating[Any]]	The 3D mask array defining the Region of Interest (ROI). Must have the same shape as data. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels used for discretization (e.g., 16, 32, 64). This determines the size of the resulting matrices.	required
distance_mask	Optional[NDArray[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).	None
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.	0
calc_glcm	bool	Whether to compute the GLCM. Default True.	True
calc_glrlm	bool	Whether to compute the GLRLM. Default True.	True
calc_ngtdm	bool	Whether to compute the NGTDM. Default True.	True
calc_ngldm	bool	Whether to compute the NGLDM. Default True.	True
calc_glszm	bool	Whether to compute the GLSZM. Default True.	True
calc_gldzm	bool	Whether to compute the GLDZM (skipping it also skips the distance transform). Default True. Disabled matrices are returned as zero-filled placeholders of minimal shape.	True

Returns:

Type

Description

dict[str, Any]

dict[str, Any]: A dictionary containing the calculated texture matrices: - '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).

Example

Calculate all texture matrices:

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)

Source code in pictologics/features/texture.py

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,
) -> 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.
            Disabled matrices are returned as zero-filled placeholders of minimal shape.

    Returns:
        dict[str, Any]: A dictionary containing the calculated texture matrices:
            - '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).

    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)
        ```
    """
    # Crop to ROI bounding box (union with distance_mask when provided) to reduce memory traffic.
    data_c, mask_c, distmask_c = _maybe_crop_to_bbox(data, mask, distance_mask)

    # Fast exit for empty ROI. Checked on the cropped mask so the common non-empty case
    # avoids a second full-volume scan; mask_c can still be empty when only distance_mask
    # has nonzero voxels.
    if not bool(np.any(mask_c != 0)):
        return {
            "glcm": np.zeros((13, n_bins, n_bins), dtype=np.uint64),
            "glrlm": np.zeros((13, n_bins, 1), dtype=np.uint64),
            "ngtdm_s": np.zeros((n_bins,), dtype=np.float64),
            "ngtdm_n": np.zeros((n_bins,), dtype=np.float64),
            "ngldm": np.zeros((n_bins, 27), dtype=np.uint64),
            "glszm": np.zeros((n_bins, 1), dtype=np.uint32),
            "gldzm": np.zeros((n_bins, 1), dtype=np.uint32),
        }

    # Use a compact binary mask representation for kernels. Label values are
    # membership markers, not weights.
    mask_u8 = (mask_c != 0).astype(np.uint8)
    # 1. Local Features (GLCM, GLRLM, NGTDM, NGLDM)
    if calc_glcm or calc_glrlm or calc_ngtdm or calc_ngldm:
        # Pre-cast data to smallest possible int type (0-based)
        # Input data is 1-based, so we subtract 1. Fused into a single pass via
        # out=, writing straight into the target-dtype buffer instead of allocating
        # a same-dtype-as-input intermediate for the subtraction and casting it after.
        if n_bins <= 256:
            data_int = np.empty(data_c.shape, dtype=np.uint8)
        else:
            data_int = np.empty(data_c.shape, dtype=np.int32)
        np.subtract(data_c, 1, out=data_int, casting="unsafe")

        try:
            n_threads = int(numba.config.NUMBA_NUM_THREADS)
        except (ValueError, TypeError):
            n_threads = 1  # Fallback

        glcm, glrlm, ngtdm_s, ngtdm_n, ngldm = _calculate_local_features_numba(
            data_int,
            mask_u8,
            n_bins,
            calc_glcm=calc_glcm,
            calc_glrlm=calc_glrlm,
            calc_ngtdm=calc_ngtdm,
            calc_ngldm=calc_ngldm,
            offsets_26=OFFSETS_26,
            directions_13=DIRECTIONS_13,
            ngldm_alpha=ngldm_alpha,
            n_threads=n_threads,
        )
    else:
        glcm = np.zeros((13, n_bins, n_bins), dtype=np.uint64)
        glrlm = np.zeros((13, n_bins, 1), dtype=np.uint64)
        ngtdm_s = np.zeros((n_bins,), dtype=np.float64)
        ngtdm_n = np.zeros((n_bins,), dtype=np.float64)
        ngldm = np.zeros((n_bins, 27), dtype=np.uint64)

    # 2. Zone Features (GLSZM, GLDZM)
    if calc_glszm or calc_gldzm:
        if calc_gldzm:
            # Pre-calculate distance map for GLDZM
            # Use distance_mask if provided, else mask
            d_mask = (distmask_c != 0).astype(np.uint8) if distmask_c is not None else mask_u8

            # Distance-to-border map, with the image border treated as an edge.
            mask_bool = d_mask > 0
            dist_map = _gldzm_distance_map(mask_bool)
        else:
            # The zone kernel needs a distance array either way; match the dtype AND
            # layout of the real distance map (an int32 strided view of a padded
            # array) so no extra JIT specialization is compiled. Never read when
            # calc_gldzm is False.
            dist_map = np.zeros(tuple(s + 2 for s in data_c.shape), dtype=np.int32)[
                1:-1, 1:-1, 1:-1
            ]

        glszm, gldzm = calculate_zone_features(
            data_c,
            mask_u8,
            dist_map,
            n_bins,
            calc_glszm=calc_glszm,
            calc_gldzm=calc_gldzm,
        )
    else:
        glszm = np.zeros((n_bins, 1), dtype=np.uint32)
        gldzm = np.zeros((n_bins, 1), dtype=np.uint32)

    return {
        "glcm": glcm,
        "glrlm": glrlm,
        "ngtdm_s": ngtdm_s,
        "ngtdm_n": ngtdm_n,
        "ngldm": ngldm,
        "glszm": glszm,
        "gldzm": gldzm,
    }

calculate_glcm_features(data, mask, n_bins, glcm_matrix=None)

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.

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 numpy import typing as npt # ... assuming data and mask defined ... features = calculate_glcm_features(data, mask, n_bins=32) print(features['contrast_ACUI']) ``` 12.5

Source code in pictologics/features/texture.py

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,
) -> 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.

        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 numpy import typing as npt
            # ... assuming data and mask defined ...
            features = calculate_glcm_features(data, mask, n_bins=32)
            print(features['contrast_ACUI'])
            ```
            12.5
    """
    if glcm_matrix is None:
        # Standalone path: crop to the ROI bbox and rebind, so the kernel and the
        # Ng_eff scan below run on the cropped arrays instead of the full volume.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        if mask.dtype == np.uint8:
            # Normalize to C-contiguous: the crop above yields a strided view, which
            # would otherwise compile a second (strided) specialization of the kernel.
            mask_u8 = np.ascontiguousarray(mask)
        else:
            mask_u8 = (mask != 0).astype(np.uint8)

        if n_bins <= 256:
            data_int = np.empty(data.shape, dtype=np.uint8)
        else:
            data_int = np.empty(data.shape, dtype=np.int32)
        np.subtract(data, 1, out=data_int, casting="unsafe")

        # Determine n_threads for JIT call
        try:
            n_threads = int(numba.config.NUMBA_NUM_THREADS)
        except (ValueError, TypeError):
            n_threads = 1  # Fallback

        # Call combined kernel to calculate only GLCM
        glcm, _, _, _, _ = _calculate_local_features_numba(
            data_int,
            mask_u8,
            n_bins,
            calc_glcm=True,
            calc_glrlm=False,
            calc_ngtdm=False,
            calc_ngldm=False,
            offsets_26=OFFSETS_26,
            directions_13=DIRECTIONS_13,
            ngldm_alpha=0,
            n_threads=n_threads,
        )
    else:
        glcm = glcm_matrix

    # 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 {}

    P = glcm_sym / total_sum

    # Indices (0-based from np.indices, convert to 1-based for IBSI)
    i_idx, j_idx = np.indices((n_bins, n_bins))
    I = i_idx + 1  # noqa: E741
    J = j_idx + 1

    features = {}

    # Joint Maximum - GYBY
    features["joint_maximum_GYBY"] = np.max(P)

    # Joint Average - 60VM
    features["joint_average_60VM"] = np.sum(I * P)

    # Joint Variance - UR99
    mu = features["joint_average_60VM"]
    features["joint_variance_UR99"] = np.sum(((I - mu) ** 2) * P)

    # Joint Entropy - TU9B
    mask_p = P > 0
    features["joint_entropy_TU9B"] = -np.sum(P[mask_p] * np.log2(P[mask_p]))

    # Difference Average - TF7R
    k_diff = np.abs(I - J)
    features["difference_average_TF7R"] = np.sum(k_diff * P)

    # Optimized using bincount
    k_diff_flat = k_diff.ravel().astype(np.int32)
    P_flat = P.ravel()
    p_diff = np.bincount(k_diff_flat, weights=P_flat)

    mu_diff = features["difference_average_TF7R"]
    k_vals = np.arange(n_bins)
    features["difference_variance_D3YU"] = np.sum(((k_vals - mu_diff) ** 2) * p_diff)

    # Difference Entropy - NTRS
    mask_pd = p_diff > 0
    features["difference_entropy_NTRS"] = -np.sum(p_diff[mask_pd] * np.log2(p_diff[mask_pd]))

    # Sum Average - ZGXS
    k_sum_grid = I + J

    # Optimized using bincount
    k_sum_flat = k_sum_grid.ravel().astype(np.int32)
    # P_flat is already defined in Difference Variance block
    p_sum_full = np.bincount(k_sum_flat, weights=P_flat)

    # Slice from 2.
    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)

    # Sum Variance - OEEB
    mu_sum = features["sum_average_ZGXS"]
    features["sum_variance_OEEB"] = np.sum(((k_vals_sum - mu_sum) ** 2) * p_sum)

    # Sum Entropy - P6QZ
    mask_ps = p_sum > 0
    features["sum_entropy_P6QZ"] = -np.sum(p_sum[mask_ps] * np.log2(p_sum[mask_ps]))

    # Angular Second Moment (Energy) - 8ZQL
    features["angular_second_moment_8ZQL"] = np.sum(P**2)

    # Contrast - ACUI
    sq_diff = (I - J) ** 2
    features["contrast_ACUI"] = np.sum(sq_diff * P)

    # Dissimilarity - 8S9J
    features["dissimilarity_8S9J"] = np.sum(k_diff * P)

    # Inverse Difference - IB1Z
    features["inverse_difference_IB1Z"] = np.sum(P / (1 + k_diff))

    # Ng_eff is the ROI grey-level span; only ROI min/max are needed. Fused single-pass
    # kernel: no bbox rescan and no boolean-gather temporaries.
    roi_span = roi_min_max(data, mask)
    if roi_span is not None:
        Ng_eff = int(roi_span[1] - roi_span[0] + 1)
    else:
        Ng_eff = 1  # Fallback

    # Normalised Inverse Difference - NDRX
    features["normalised_inverse_difference_NDRX"] = np.sum(P / (1 + k_diff / Ng_eff))

    # Inverse Difference Moment - WF0Z
    features["inverse_difference_moment_WF0Z"] = np.sum(P / (1 + sq_diff))

    # Normalised Inverse Difference Moment - 1QCO
    features["normalised_inverse_difference_moment_1QCO"] = np.sum(
        P / (1 + sq_diff / (Ng_eff**2))
    )

    # Inverse Variance - E8JP
    mask_neq = I != J
    features["inverse_variance_E8JP"] = np.sum(P[mask_neq] / ((I[mask_neq] - J[mask_neq]) ** 2))

    # Correlation - NI2N
    term1 = np.sum((I - mu) * (J - mu) * P)
    if features["joint_variance_UR99"] != 0:
        features["correlation_NI2N"] = term1 / features["joint_variance_UR99"]
    else:
        features["correlation_NI2N"] = 1.0  # Or NaN? IBSI doesn't specify for 0 variance.

    # Autocorrelation - QWB0
    features["autocorrelation_QWB0"] = np.sum(I * J * P)

    # Cluster Tendency - DG8W
    sum_diff2mu = I + J - 2 * mu
    features["cluster_tendency_DG8W"] = np.sum((sum_diff2mu**2) * P)

    # Cluster Shade - 7NFM
    features["cluster_shade_7NFM"] = np.sum((sum_diff2mu**3) * P)

    # Cluster Prominence - AE86
    features["cluster_prominence_AE86"] = np.sum((sum_diff2mu**4) * P)

    # Information Correlation 1 - R8DG
    HXY = features["joint_entropy_TU9B"]
    p_x = np.sum(P, axis=1)
    mask_px = p_x > 0
    HX = -np.sum(p_x[mask_px] * np.log2(p_x[mask_px]))

    HXY1 = -np.sum(P[mask_p] * np.log2(p_x[I[mask_p] - 1] * p_x[J[mask_p] - 1]))

    if HX != 0:
        features["information_correlation_1_R8DG"] = (HXY - HXY1) / HX
    else:
        features["information_correlation_1_R8DG"] = np.nan

    # Information Correlation 2 - JN9H
    P_prod = np.outer(p_x, p_x)
    mask_prod = P_prod > 0
    HXY2 = -np.sum(P_prod[mask_prod] * np.log2(P_prod[mask_prod]))

    features["information_correlation_2_JN9H"] = np.sqrt(1 - np.exp(-2 * (HXY2 - HXY)))
    return features

calculate_gldzm_features(data, mask, n_bins, gldzm_matrix=None, distance_mask=None)

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.

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels.	required
mask	NDArray[floating[Any]]	The 3D mask array defining the ROI. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels.	required
gldzm_matrix	Optional[NDArray[floating[Any]]]	Pre-calculated GLDZM matrix.	None
distance_mask	Optional[NDArray[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).	None

Returns:

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

Source code in pictologics/features/texture.py

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'.
    """
    if gldzm_matrix is None:
        # Standalone path: crop to the bbox of mask ∪ distance_mask and rebind. The
        # union keeps the distance map identical (the distance region is fully inside
        # the bbox), and the voxel count below runs on the cropped mask.
        data, mask, distance_mask = _maybe_crop_to_bbox(data, mask, distance_mask)
        mask_u8 = (mask != 0).astype(np.uint8)

        # Calculate distance map
        # Use distance_mask if provided, else mask
        d_mask = distance_mask if distance_mask is not None else mask

        # Distance-to-border map, with the image border treated as an edge.
        mask_bool = d_mask > 0
        dist_map = _gldzm_distance_map(mask_bool)

        _, gldzm = calculate_zone_features(  # type: ignore[arg-type]
            data,
            mask_u8,
            dist_map,
            n_bins,
            calc_glszm=False,
            calc_gldzm=True,
        )
    else:
        gldzm = gldzm_matrix

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

    P = gldzm / N_zones

    n_g, n_d = gldzm.shape
    i_idx, j_idx = np.indices((n_g, n_d))
    I = i_idx + 1  # noqa: E741
    J = j_idx + 1  # Distance
    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

calculate_glrlm_features(data, mask, n_bins, glrlm_matrix=None)

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

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels.	required
mask	NDArray[floating[Any]]	The 3D mask array defining the ROI. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels.	required
glrlm_matrix	Optional[NDArray[floating[Any]]]	Pre-calculated GLRLM matrix.	None

Returns:

Type	Description
dict[str, float]	dict[str, float]: A dictionary of calculated GLRLM features. Example keys: 'short_runs_emphasis_22OV', 'grey_level_non_uniformity_R5YN'.

Source code in pictologics/features/texture.py

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,
) -> 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.

    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 and rebind, so the kernel and the
        # run-percentage voxel count below run on the cropped arrays.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        if mask.dtype == np.uint8:
            # Normalize to C-contiguous: the crop above yields a strided view, which
            # would otherwise compile a second (strided) specialization of the kernel.
            mask_u8 = np.ascontiguousarray(mask)
        else:
            mask_u8 = (mask != 0).astype(np.uint8)

        if n_bins <= 256:
            data_int = np.empty(data.shape, dtype=np.uint8)
        else:
            data_int = np.empty(data.shape, dtype=np.int32)  # pragma: no cover
        np.subtract(data, 1, out=data_int, casting="unsafe")

        # Determine n_threads for JIT call
        try:
            n_threads = int(numba.config.NUMBA_NUM_THREADS)
        except (ValueError, TypeError):
            n_threads = 1  # Fallback

        # Call combined kernel
        _, glrlm, _, _, _ = _calculate_local_features_numba(
            data_int,
            mask_u8,
            n_bins,
            calc_glcm=False,
            calc_glrlm=True,
            calc_ngtdm=False,
            calc_ngldm=False,
            offsets_26=OFFSETS_26,
            directions_13=DIRECTIONS_13,
            ngldm_alpha=0,
            n_threads=n_threads,
        )
    else:
        glrlm = glrlm_matrix

    # 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 {}

    P = glrlm / N_runs

    n_g, n_r = glrlm.shape
    i_idx, j_idx = np.indices((n_g, n_r))
    I = i_idx + 1  # noqa: E741
    J = j_idx + 1
    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)
    n_dirs = 13  # Fixed for 3D
    features["run_percentage_9ZK5"] = N_runs / (n_voxels * n_dirs)

    # 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

calculate_glszm_features(data, mask, n_bins, glszm_matrix=None)

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.

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels.	required
mask	NDArray[floating[Any]]	The 3D mask array defining the ROI. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels.	required
glszm_matrix	Optional[NDArray[floating[Any]]]	Pre-calculated GLSZM matrix.	None

Returns:

Type	Description
dict[str, float]	dict[str, float]: A dictionary of calculated GLSZM features. Example keys: 'small_zone_emphasis_P001', 'zone_percentage_P30P'.

Source code in pictologics/features/texture.py

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 and rebind, so the zone kernel and the
        # zone-percentage voxel count below run on the cropped arrays.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        if mask.dtype == np.uint8:
            # Normalize to C-contiguous: the crop above yields a strided view, which
            # would otherwise compile a second (strided) specialization of the kernel.
            mask_u8 = np.ascontiguousarray(mask)
        else:
            mask_u8 = (mask != 0).astype(np.uint8)

        # We need dist_map for the combined kernel signature, even if unused for
        # GLSZM; match the dtype AND layout of the real distance map (an int32
        # strided view of a padded array) so no extra JIT specialization is compiled.
        dummy_dist: npt.NDArray[Any] = np.zeros(tuple(s + 2 for s in data.shape), dtype=np.int32)[
            1:-1, 1:-1, 1:-1
        ]

        glszm, _ = calculate_zone_features(
            data,
            mask_u8,
            dummy_dist,
            n_bins,
            calc_glszm=True,
            calc_gldzm=False,  # type: ignore[arg-type]
        )
    else:
        glszm = glszm_matrix

    # The GLSZM is extremely sparse in the zone-size dimension (a single large zone
    # can push n_s into the tens of thousands while only a few hundred cells are
    # non-zero). Operating on the non-zero cells avoids materialising the dense
    # (n_g, n_s) index grids and is arithmetically identical (every term carries a
    # factor of P, so zero cells contribute exactly zero).
    n_g, n_s = glszm.shape
    gl_idx, sz_idx = np.nonzero(glszm)
    if gl_idx.size == 0:
        return {}

    c = glszm[gl_idx, sz_idx].astype(np.float64)
    N_zones = c.sum()

    I = (gl_idx + 1).astype(np.float64)  # noqa: E741
    J = (sz_idx + 1).astype(np.float64)  # Zone size
    I2 = I * I
    J2 = J * J
    P = c / N_zones

    features = {}

    # Small Zone Emphasis (SZE) - P001
    features["small_zone_emphasis_P001"] = np.sum(P / J2)

    # Large Zone Emphasis (LZE) - 48P8
    features["large_zone_emphasis_48P8"] = np.sum(P * J2)

    # Grey Level Non-Uniformity (GLNU) - JNSA
    s_i = np.bincount(gl_idx, weights=c)
    sum_si2 = np.sum(s_i**2)
    features["grey_level_non_uniformity_JNSA"] = sum_si2 / N_zones

    # Normalised Grey Level Non-Uniformity (GLNN) - Y1RO
    features["normalised_grey_level_non_uniformity_Y1RO"] = sum_si2 / (N_zones**2)

    # Zone Size Non-Uniformity (ZSNU) - 4JP3
    s_j = np.bincount(sz_idx, weights=c)
    sum_sj2 = np.sum(s_j**2)
    features["zone_size_non_uniformity_4JP3"] = sum_sj2 / N_zones

    # Normalised Zone Size Non-Uniformity (ZSNN) - VB3A
    features["normalised_zone_size_non_uniformity_VB3A"] = sum_sj2 / (N_zones**2)

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

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

    # Zone Size Variance (ZSV) - 3NSA
    mu_j = np.sum(J * P)
    features["zone_size_variance_3NSA"] = np.sum(((J - mu_j) ** 2) * P)

    # Zone Size Entropy (ZSE) - GU8N
    features["zone_size_entropy_GU8N"] = -np.sum(P * np.log2(P))

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

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

    # Small Zone Low Grey Level Emphasis (SZLGLE) - 5RAI
    features["small_zone_low_grey_level_emphasis_5RAI"] = np.sum(P / (I2 * J2))

    # Small Zone High Grey Level Emphasis (SZHGLE) - HW1V
    features["small_zone_high_grey_level_emphasis_HW1V"] = np.sum(P * I2 / J2)

    # Large Zone Low Grey Level Emphasis (LZLGLE) - YH51
    features["large_zone_low_grey_level_emphasis_YH51"] = np.sum(P * J2 / I2)

    # Large Zone High Grey Level Emphasis (LZHGLE) - J17V
    features["large_zone_high_grey_level_emphasis_J17V"] = np.sum(P * I2 * J2)

    return features

calculate_ngldm_features(data, mask, n_bins, ngldm_matrix=None, ngldm_alpha=0)

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 α).

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels.	required
mask	NDArray[floating[Any]]	The 3D mask array defining the ROI. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels.	required
ngldm_matrix	Optional[NDArray[floating[Any]]]	Pre-calculated NGLDM matrix.	None
ngldm_alpha	int	The coarseness parameter α. Two grey levels are considered dependent if their absolute difference is ≤ α. Default is 0 (exact match, IBSI standard).	0

Returns:

Type	Description
dict[str, float]	dict[str, float]: A dictionary of calculated NGLDM features. Example keys: 'low_dependence_emphasis_SODN', 'dependence_count_entropy_FCBV'.

Source code in pictologics/features/texture.py

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,
) -> 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).

    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 and rebind, so the kernel and the
        # dependence-count-percentage voxel count below run on the cropped arrays.
        data, mask, _ = _maybe_crop_to_bbox(data, mask, None)
        if mask.dtype == np.uint8:
            # Normalize to C-contiguous: the crop above yields a strided view, which
            # would otherwise compile a second (strided) specialization of the kernel.
            mask_u8 = np.ascontiguousarray(mask)
        else:
            mask_u8 = (mask != 0).astype(np.uint8)

        if n_bins <= 256:
            data_int = np.empty(data.shape, dtype=np.uint8)
        else:
            data_int = np.empty(data.shape, dtype=np.int32)  # pragma: no cover
        np.subtract(data, 1, out=data_int, casting="unsafe")

        # Determine n_threads for JIT call
        try:
            n_threads = int(numba.config.NUMBA_NUM_THREADS)
        except (ValueError, TypeError):
            n_threads = 1  # Fallback

        _, _, _, _, ngldm = _calculate_local_features_numba(
            data_int,
            mask_u8,
            n_bins,
            calc_glcm=False,
            calc_glrlm=False,
            calc_ngtdm=False,
            calc_ngldm=True,
            offsets_26=OFFSETS_26,
            directions_13=DIRECTIONS_13,
            ngldm_alpha=ngldm_alpha,
            n_threads=n_threads,
        )
    else:
        ngldm = ngldm_matrix

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

    P = ngldm / N_s

    n_g, n_d = ngldm.shape
    i_idx, j_idx = np.indices((n_g, n_d))
    I = i_idx + 1  # noqa: E741
    J = j_idx + 1  # Dependence count
    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

calculate_ngtdm_features(data, mask, n_bins, ngtdm_matrices=None)

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.

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	The 3D image array containing discretised grey levels.	required
mask	NDArray[floating[Any]]	The 3D mask array defining the ROI. Nonzero values indicate ROI membership.	required
n_bins	int	The number of grey levels.	required
ngtdm_matrices	Optional[tuple[NDArray[floating[Any]], NDArray[floating[Any]]]]	Pre-calculated NGTDM matrices (sum of absolute differences s, and count n).	None

Returns:

Type	Description
dict[str, float]	dict[str, float]: A dictionary of calculated NGTDM features. Example keys: 'coarseness_QCDE', 'contrast_65HE', 'busyness_NQ30'.

Source code in pictologics/features/texture.py

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

    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)
        if mask.dtype == np.uint8:
            # Normalize to C-contiguous: the crop above yields a strided view, which
            # would otherwise compile a second (strided) specialization of the kernel.
            mask_u8 = np.ascontiguousarray(mask)
        else:
            mask_u8 = (mask != 0).astype(np.uint8)

        if n_bins <= 256:
            data_int = np.empty(data.shape, dtype=np.uint8)
        else:
            data_int = np.empty(data.shape, dtype=np.int32)  # pragma: no cover
        np.subtract(data, 1, out=data_int, casting="unsafe")

        # Determine n_threads for JIT call
        try:
            n_threads = int(numba.config.NUMBA_NUM_THREADS)
        except (ValueError, TypeError):
            n_threads = 1  # Fallback

        _, _, s, n, _ = _calculate_local_features_numba(
            data_int,
            mask_u8,
            n_bins,
            calc_glcm=False,
            calc_glrlm=False,
            calc_ngtdm=True,
            calc_ngldm=False,
            offsets_26=OFFSETS_26,
            directions_13=DIRECTIONS_13,
            ngldm_alpha=0,
            n_threads=n_threads,
        )
    else:
        s, n = ngtdm_matrices

    # s[i] is sum of absolute differences for grey level i+1
    # n[i] is number of voxels of grey level i+1 with valid neighborhood

    N_vp = np.sum(n)
    if N_vp == 0:
        return {}

    p = n / N_vp

    # Indices
    i_idx = np.arange(n_bins)
    I = i_idx + 1  # noqa: E741

    features = {}

    # Filter for non-zero probabilities (required for Busyness, Complexity, Strength)
    mask_p = p > 0
    p_nz = p[mask_p]
    s_nz = s[mask_p]
    I_nz = I[mask_p]

    # Pairwise grids over the non-zero grey levels, shared by Contrast/Complexity/Strength.
    Pi, Pj = np.meshgrid(p_nz, p_nz, indexing="ij")
    Ii, Ij = np.meshgrid(I_nz, I_nz, indexing="ij")
    Si, Sj = np.meshgrid(s_nz, s_nz, indexing="ij")

    # Coarseness - QCDE
    sum_ps = np.sum(p_nz * s_nz)
    if sum_ps > 1e-10:
        features["coarseness_QCDE"] = 1 / sum_ps
    else:
        features["coarseness_QCDE"] = 1e6

    # Contrast - 65HE
    Ng_p = len(p_nz)

    if Ng_p > 1:
        # Term 1: Dynamic range variance
        term1_sum = np.sum(Pi * Pj * ((Ii - Ij) ** 2))
        term1 = term1_sum / (Ng_p * (Ng_p - 1))

        # Term 2: Intensity change
        sum_s = np.sum(s)
        term2 = sum_s / N_vp

        features["contrast_65HE"] = term1 * term2
    else:
        features["contrast_65HE"] = 0.0

    # Busyness - NQ30
    IPi = I_nz * p_nz

    # Grid
    IPi_grid, IPj_grid = np.meshgrid(IPi, IPi, indexing="ij")
    denom_busyness = np.sum(np.abs(IPi_grid - IPj_grid))

    if denom_busyness > 1e-10:
        features["busyness_NQ30"] = sum_ps / denom_busyness
    else:
        features["busyness_NQ30"] = 0.0

    # Complexity - HDEZ
    denom_comp = Pi + Pj
    term_comp = np.abs(Ii - Ij) * (Pi * Si + Pj * Sj) / denom_comp

    features["complexity_HDEZ"] = (1 / N_vp) * np.sum(term_comp)

    # Strength - 1X9X
    sum_s = np.sum(s)

    term_str = (Pi + Pj) * ((Ii - Ij) ** 2)
    sum_term_str = np.sum(term_str)

    if sum_s > 1e-10:
        features["strength_1X9X"] = sum_term_str / sum_s
    else:
        features["strength_1X9X"] = 0.0

    return features

calculate_zone_features(data, mask, dist_map, n_bins, calc_glszm=True, calc_gldzm=True)

Wrapper for _calculate_zone_features_numba with buffer pooling.

This function manages pre-allocated buffers to reduce memory allocation overhead for repeated calls (e.g., during batch processing).

Parameters:

Name	Type	Description	Default
data	NDArray[floating[Any]]	3D discretized image data.	required
mask	NDArray[floating[Any]]	3D mask array (not modified - copied internally by JIT function).	required
dist_map	NDArray[floating[Any]]	3D distance map for GLDZM.	required
n_bins	int	Number of grey level bins.	required
calc_glszm	bool	Whether to calculate GLSZM.	True
calc_gldzm	bool	Whether to calculate GLDZM.	True

Returns:

Type	Description
tuple[NDArray[floating[Any]], NDArray[floating[Any]]]	Tuple of (glszm, gldzm) matrices.

Example

import numpy as np
from pictologics.features.texture import calculate_zone_features

rng = np.random.default_rng(0)
data = rng.integers(1, 5, size=(8, 8, 8)).astype(np.float64)
mask = np.ones((8, 8, 8), dtype=np.uint8)
dist_map = np.zeros((8, 8, 8), dtype=np.int32)

glszm, gldzm = calculate_zone_features(data, mask, dist_map, n_bins=4, calc_gldzm=False)
print(int(glszm.sum()))
# 13

Source code in pictologics/features/texture.py

def calculate_zone_features(
    data: npt.NDArray[np.floating[Any]],
    mask: npt.NDArray[np.floating[Any]],
    dist_map: npt.NDArray[np.floating[Any]],
    n_bins: int,
    calc_glszm: bool = True,
    calc_gldzm: bool = True,
) -> tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]:
    """
    Wrapper for _calculate_zone_features_numba with buffer pooling.

    This function manages pre-allocated buffers to reduce memory allocation
    overhead for repeated calls (e.g., during batch processing).

    Args:
        data: 3D discretized image data.
        mask: 3D mask array (not modified - copied internally by JIT function).
        dist_map: 3D distance map for GLDZM.
        n_bins: Number of grey level bins.
        calc_glszm: Whether to calculate GLSZM.
        calc_gldzm: Whether to calculate GLDZM.

    Returns:
        Tuple of (glszm, gldzm) matrices.

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

        rng = np.random.default_rng(0)
        data = rng.integers(1, 5, size=(8, 8, 8)).astype(np.float64)
        mask = np.ones((8, 8, 8), dtype=np.uint8)
        dist_map = np.zeros((8, 8, 8), dtype=np.int32)

        glszm, gldzm = calculate_zone_features(data, mask, dist_map, n_bins=4, calc_gldzm=False)
        print(int(glszm.sum()))
        # 13
        ```
    """
    # For zone features, the worst-case number of zones is bounded by ROI voxel count.
    # Sizing buffers to full image volume is extremely costly for sparse ROIs.
    max_zones = int(np.count_nonzero(mask))
    if max_zones < 1:
        max_zones = 1

    # Get pre-allocated buffers from pool
    pool = _ZoneBufferPool.get_instance()
    res_gl, res_size, res_dist, stack = pool.get_buffers(max_zones)

    try:
        n_chunks = int(numba.config.NUMBA_NUM_THREADS)
    except (ValueError, TypeError):
        n_chunks = 1  # Fallback

    # The parallel kernel pays a fixed threading-dispatch cost per call (~1.5 ms at
    # 14 threads); below the measured crossover volume the serial kernel is faster.
    if n_chunks == 1 or mask.size < 1 << 17:
        return cast(
            tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]],
            _calculate_zone_features_serial_numba(
                data,
                mask,
                dist_map,
                n_bins,
                res_gl,
                res_size,
                res_dist,
                stack,
                calc_glszm,
                calc_gldzm,
            ),
        )

    return cast(
        tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]],
        _calculate_zone_features_numba(
            data,
            mask,
            dist_map,
            n_bins,
            res_gl,
            res_size,
            res_dist,
            stack,
            n_chunks,
            calc_glszm,
            calc_gldzm,
        ),
    )