Intensity Features API

pictologics.features.intensity

Intensity Feature Extraction Module

This module provides functions for calculating First Order Statistics (Intensity) features from medical images. It implements the Image Biomarker Standardisation Initiative (IBSI) compliant algorithms.

Key Features:

• First Order Statistics: Mean, Variance, Skewness, Kurtosis, Percentiles, etc.
• Intensity Histogram: Features based on discretised intensity histograms.
• Intensity-Volume Histogram (IVH): Volume fractions and intensity fractions, AUC.
• Spatial Intensity: Moran's I and Geary's C (spatial autocorrelation) [Optimized].
• Local Intensity: Local and Global Intensity Peaks.

Optimization:

Uses numba for JIT compilation, with parallel execution for computationally intensive spatial feature calculations.

calculate_intensity_features(values)

Calculate intensity-based features (First Order Statistics) as defined in IBSI 4.1.

Computes 18 statistical features from the intensity values within the ROI: mean, variance, skewness, kurtosis, median, min/max, percentiles (10th, 90th), interquartile range, range, MAD variants, coefficient of variation, energy, RMS.

Parameters:

Name	Type	Description	Default
values	NDArray[floating[Any]]	1D array of intensity values from the ROI (after mask application).	required

Returns:

Type	Description
dict[str, float]	Dictionary mapping feature names (with IBSI codes) to computed values.
dict[str, float]	Empty dict if input is empty.

Example

Calculate features from an ROI:

from pictologics.features.intensity import calculate_intensity_features
from pictologics.preprocessing import apply_mask

# Get values within ROI
roi_values = apply_mask(image, mask)

# Calculate features
features = calculate_intensity_features(roi_values)
print(features["mean_intensity_Q4LE"])

Source code in pictologics/features/intensity.py

def calculate_intensity_features(
    values: npt.NDArray[np.floating[Any]],
) -> dict[str, float]:
    """
    Calculate intensity-based features (First Order Statistics) as defined in IBSI 4.1.

    Computes 18 statistical features from the intensity values within the ROI:
    mean, variance, skewness, kurtosis, median, min/max, percentiles (10th, 90th),
    interquartile range, range, MAD variants, coefficient of variation, energy, RMS.

    Args:
        values: 1D array of intensity values from the ROI (after mask application).

    Returns:
        Dictionary mapping feature names (with IBSI codes) to computed values.
        Empty dict if input is empty.

    Example:
        Calculate features from an ROI:

        ```python
        from pictologics.features.intensity import calculate_intensity_features
        from pictologics.preprocessing import apply_mask

        # Get values within ROI
        roi_values = apply_mask(image, mask)

        # Calculate features
        features = calculate_intensity_features(roi_values)
        print(features["mean_intensity_Q4LE"])
        ```
    """
    if len(values) == 0:
        return {}

    features: dict[str, float] = {}

    # 4.1.1 Mean intensity (Q4LE)
    mean_val = np.mean(values)
    features["mean_intensity_Q4LE"] = float(mean_val)

    # 4.1.2 Intensity variance (ECT3)
    var_val = float(np.var(values, ddof=0))
    features["intensity_variance_ECT3"] = float(var_val)

    # 4.1.3 Intensity skewness (KE2A)
    if var_val == 0.0:
        features["intensity_skewness_KE2A"] = np.nan
        features["intensity_kurtosis_IPH6"] = np.nan
    else:
        m2, m3, m4 = _central_moments_2_3_4(values, float(mean_val))
        denom = m2**1.5
        if denom != 0.0:
            features["intensity_skewness_KE2A"] = float(m3 / denom)
            features["intensity_kurtosis_IPH6"] = float((m4 / (m2 * m2)) - 3.0)

    # 4.1.5 Median intensity (Y12H)
    median_val = np.median(values)
    features["median_intensity_Y12H"] = float(median_val)

    # 4.1.6 Minimum intensity (1GSF)
    min_val = np.min(values)
    features["minimum_intensity_1GSF"] = float(min_val)

    p10, p25, p75, p90 = np.percentile(values, [10, 25, 75, 90])
    features["10th_intensity_percentile_QG58"] = float(p10)
    features["90th_intensity_percentile_8DWT"] = float(p90)

    # 4.1.9 Maximum intensity (84IY)
    max_val = np.max(values)
    features["maximum_intensity_84IY"] = float(max_val)

    # 4.1.10 Intensity interquartile range (SALO)
    features["intensity_interquartile_range_SALO"] = float(p75 - p25)

    # 4.1.11 Intensity range (2OJQ)
    features["intensity_range_2OJQ"] = float(max_val - min_val)

    # 4.1.12 Mean absolute deviation (4FUA)
    features["intensity_mean_absolute_deviation_4FUA"] = float(
        _mean_abs_dev(values, float(mean_val))
    )

    # 4.1.13 Robust mean absolute deviation (1128)
    features["intensity_robust_mean_absolute_deviation_1128"] = float(
        _robust_mean_abs_dev(values, float(p10), float(p90))
    )

    # 4.1.14 Median absolute deviation (N72L)
    features["intensity_median_absolute_deviation_N72L"] = float(
        _mean_abs_dev(values, float(median_val))
    )

    # 4.1.15 Coefficient of variation (7TET)
    if mean_val != 0:
        features["intensity_coefficient_of_variation_7TET"] = float(
            np.sqrt(var_val) / mean_val
        )
    else:
        features["intensity_coefficient_of_variation_7TET"] = np.nan

    # 4.1.16 Quartile coefficient of dispersion (9S40)
    if (p75 + p25) != 0:
        features["intensity_quartile_coefficient_of_dispersion_9S40"] = float(
            (p75 - p25) / (p75 + p25)
        )
    else:
        features["intensity_quartile_coefficient_of_dispersion_9S40"] = np.nan

    # 4.1.17 Energy (N8CA)
    energy = float(np.dot(values, values))
    features["intensity_energy_N8CA"] = energy

    # 4.1.18 Root mean square (5ZWQ)
    features["root_mean_square_intensity_5ZWQ"] = float(np.sqrt(energy / len(values)))

    return features

calculate_intensity_histogram_features(discretised_values)

Calculate intensity histogram features as defined in IBSI 4.2.

Computes features from the discretised intensity histogram including mean, variance, skewness, kurtosis, mode, entropy, uniformity, and gradient features.

Parameters:

Name	Type	Description	Default
discretised_values	NDArray[floating[Any]]	1D array of discretised intensity values (after binning).	required

Returns:

Type	Description
dict[str, float]	Dictionary mapping feature names (with IBSI codes) to computed values.
dict[str, float]	Empty dict if input is empty.

Source code in pictologics/features/intensity.py

def calculate_intensity_histogram_features(
    discretised_values: npt.NDArray[np.floating[Any]],
) -> dict[str, float]:
    """
    Calculate intensity histogram features as defined in IBSI 4.2.

    Computes features from the discretised intensity histogram including
    mean, variance, skewness, kurtosis, mode, entropy, uniformity, and gradient features.

    Args:
        discretised_values: 1D array of discretised intensity values (after binning).

    Returns:
        Dictionary mapping feature names (with IBSI codes) to computed values.
        Empty dict if input is empty.
    """
    if len(discretised_values) == 0:
        return {}

    features: dict[str, float] = {}

    disc = np.asarray(discretised_values)
    n = disc.size

    # Support negative values by shifting for bincount compatibility
    min_val_i = int(np.min(disc))
    max_val_i = int(np.max(disc))

    shifted = disc.astype(np.int64) - min_val_i
    counts_full = np.bincount(shifted, minlength=(max_val_i - min_val_i + 1))
    total = float(n)
    p = counts_full[counts_full > 0].astype(np.float64) / total

    # 4.2.1 Mean discretised intensity (X6K6)
    mean_disc = float(np.mean(disc))
    features["mean_discretised_intensity_X6K6"] = float(mean_disc)

    # 4.2.2 Discretised intensity variance (CH89)
    var_disc = float(np.var(disc, ddof=0))
    features["discretised_intensity_variance_CH89"] = float(var_disc)

    # 4.2.3 Discretised intensity skewness (88K1)
    if var_disc == 0.0:
        features["discretised_intensity_skewness_88K1"] = np.nan
        features["discretised_intensity_kurtosis_C3I7"] = np.nan
    else:
        m2, m3, m4 = _central_moments_2_3_4(disc, float(mean_disc))
        denom = m2**1.5
        if denom != 0.0:
            features["discretised_intensity_skewness_88K1"] = float(m3 / denom)
            features["discretised_intensity_kurtosis_C3I7"] = float(
                (m4 / (m2 * m2)) - 3.0
            )

    p10, p25, median_val, p75, p90 = np.percentile(disc, [10, 25, 50, 75, 90])

    features["median_discretised_intensity_WIFQ"] = float(median_val)
    features["minimum_discretised_intensity_1PR8"] = float(min_val_i)
    features["10th_discretised_intensity_percentile_1PR"] = float(p10)
    features["90th_discretised_intensity_percentile_GPMT"] = float(p90)
    features["maximum_discretised_intensity_3NCY"] = float(max_val_i)

    mode_index = int(np.argmax(counts_full))
    features["intensity_histogram_mode_AMMC"] = float(min_val_i + mode_index)

    features["discretised_intensity_interquartile_range_WR0O"] = float(p75 - p25)

    features["discretised_intensity_range_5Z3W"] = float(
        features["maximum_discretised_intensity_3NCY"]
        - features["minimum_discretised_intensity_1PR8"]
    )

    features["intensity_histogram_mean_absolute_deviation_D2ZX"] = float(
        _mean_abs_dev(disc, float(mean_disc))
    )

    features["intensity_histogram_robust_mean_absolute_deviation_WRZB"] = float(
        _robust_mean_abs_dev(disc, float(p10), float(p90))
    )

    features["intensity_histogram_median_absolute_deviation_4RNL"] = float(
        _mean_abs_dev(disc, float(median_val))
    )

    if mean_disc != 0:
        features["intensity_histogram_coefficient_of_variation_CWYJ"] = float(
            np.sqrt(var_disc) / mean_disc
        )
    else:
        features["intensity_histogram_coefficient_of_variation_CWYJ"] = np.nan

    if (p75 + p25) != 0:
        features["intensity_histogram_quartile_coefficient_of_dispersion_SLWD"] = float(
            (p75 - p25) / (p75 + p25)
        )
    else:
        features["intensity_histogram_quartile_coefficient_of_dispersion_SLWD"] = np.nan

    # Vectorized entropy/uniformity
    features["discretised_intensity_entropy_TLU2"] = float(-np.sum(p * np.log2(p)))
    features["discretised_intensity_uniformity_BJ5W"] = float(np.sum(p * p))

    hist_counts = counts_full.astype(np.float64)
    if len(hist_counts) < 2:
        features["maximum_histogram_gradient_12CE"] = np.nan
        features["maximum_histogram_gradient_intensity_8E6O"] = np.nan
        features["minimum_histogram_gradient_VQB3"] = np.nan
        features["minimum_histogram_gradient_intensity_RHQZ"] = np.nan
    else:
        gradient = np.gradient(hist_counts)
        features["maximum_histogram_gradient_12CE"] = float(np.max(gradient))
        max_grad_idx = int(np.argmax(gradient))
        features["maximum_histogram_gradient_intensity_8E6O"] = float(
            min_val_i + max_grad_idx
        )
        features["minimum_histogram_gradient_VQB3"] = float(np.min(gradient))
        min_grad_idx = int(np.argmin(gradient))
        features["minimum_histogram_gradient_intensity_RHQZ"] = float(
            min_val_i + min_grad_idx
        )

    return features

calculate_ivh_features(discretised_values, bin_width=None, min_val=None, max_val=None, target_range_min=None, target_range_max=None)

Calculate Intensity-Volume Histogram (IVH) features as defined in IBSI 4.3.

Computes volume fractions at intensity thresholds, intensity values at volume fractions, and Area Under the IVH Curve (AUC).

Parameters:

Name	Type	Description	Default
discretised_values	NDArray[floating[Any]]	1D array of discretised intensity values.	required
bin_width	Optional[float]	Optional bin width used in discretisation (for physical units).	None
min_val	Optional[float]	Optional minimum value used in discretisation.	None
max_val	Optional[float]	Optional maximum value used in discretisation.	None
target_range_min	Optional[float]	Optional target range minimum for fraction calculations.	None
target_range_max	Optional[float]	Optional target range maximum for fraction calculations.	None

Returns:

Type	Description
dict[str, float]	Dictionary mapping feature names (with IBSI codes) to computed values.
dict[str, float]	Empty dict if input is empty.

Source code in pictologics/features/intensity.py

def calculate_ivh_features(
    discretised_values: npt.NDArray[np.floating[Any]],
    bin_width: Optional[float] = None,
    min_val: Optional[float] = None,
    max_val: Optional[float] = None,
    target_range_min: Optional[float] = None,
    target_range_max: Optional[float] = None,
) -> dict[str, float]:
    """
    Calculate Intensity-Volume Histogram (IVH) features as defined in IBSI 4.3.

    Computes volume fractions at intensity thresholds, intensity values at volume
    fractions, and Area Under the IVH Curve (AUC).

    Args:
        discretised_values: 1D array of discretised intensity values.
        bin_width: Optional bin width used in discretisation (for physical units).
        min_val: Optional minimum value used in discretisation.
        max_val: Optional maximum value used in discretisation.
        target_range_min: Optional target range minimum for fraction calculations.
        target_range_max: Optional target range maximum for fraction calculations.

    Returns:
        Dictionary mapping feature names (with IBSI codes) to computed values.
        Empty dict if input is empty.
    """
    if len(discretised_values) == 0:
        return {}

    features: dict[str, float] = {}
    N = len(discretised_values)

    vals = np.asarray(discretised_values)
    sorted_vals = np.sort(vals)

    # -------------------------------------------------------------------------
    # 1. Volume Fractions
    # -------------------------------------------------------------------------
    t_min = target_range_min if target_range_min is not None else min_val
    t_max = target_range_max if target_range_max is not None else max_val

    # Fallback to data min/max if still None
    if t_min is None or t_max is None:
        t_min_idx = float(sorted_vals[0])
        t_max_idx = float(sorted_vals[-1])
        val_range_idx = t_max_idx - t_min_idx

        def get_volume_fraction_at_intensity_fraction_indices(frac: float) -> float:
            threshold_idx = t_min_idx + frac * val_range_idx
            idx = int(np.searchsorted(sorted_vals, threshold_idx, side="left"))
            count = N - idx
            return float(count / N)

        features["volume_at_intensity_fraction_0.10_BC2M_10"] = (
            get_volume_fraction_at_intensity_fraction_indices(0.10)
        )
        features["volume_at_intensity_fraction_0.90_BC2M_90"] = (
            get_volume_fraction_at_intensity_fraction_indices(0.90)
        )

    else:
        # We have physical units for the target range.
        t_range = t_max - t_min

        def get_volume_fraction_at_intensity_fraction_physical(frac: float) -> float:
            threshold_val = t_min + frac * t_range  # type: ignore
            if bin_width is not None and min_val is not None:
                # Convert to index based on FBS
                mv = min_val
                bw = bin_width
                threshold_idx = np.floor((threshold_val - mv) / bw) + 1  # type: ignore
            else:
                threshold_idx = threshold_val

            idx = int(np.searchsorted(sorted_vals, threshold_idx, side="left"))
            count = N - idx
            return float(count / N)

        features["volume_at_intensity_fraction_0.10_BC2M_10"] = (
            get_volume_fraction_at_intensity_fraction_physical(0.10)
        )
        features["volume_at_intensity_fraction_0.90_BC2M_90"] = (
            get_volume_fraction_at_intensity_fraction_physical(0.90)
        )

    features[
        "volume_fraction_difference_between_intensity_0.10_and_0.90_fractions_DDTU"
    ] = float(
        features["volume_at_intensity_fraction_0.10_BC2M_10"]
        - features["volume_at_intensity_fraction_0.90_BC2M_90"]
    )

    # -------------------------------------------------------------------------
    # 2. Intensity Fractions
    # -------------------------------------------------------------------------
    def get_intensity_at_volume_fraction(vol_frac: float) -> float:
        # Fast path for standard integer bins (step=1)
        if (
            bin_width is not None
            and min_val is None
            and bin_width > 0
            and float(bin_width) == 1.0
            and np.issubdtype(sorted_vals.dtype, np.integer)
        ):
            target_count = int(np.floor(vol_frac * N))
            if target_count <= 0:
                return float(sorted_vals[-1])

            # Smallest integer threshold t such that count(vals >= t) <= target_count.
            # Let k = N - target_count. We need searchsorted(t) >= k.
            k = N - target_count
            v = int(sorted_vals[k - 1])
            t = v + 1
            vmax = int(sorted_vals[-1])
            if t > vmax:
                t = vmax
            return float(t)

        # Determine candidates
        if bin_width is not None:
            g_min = min_val if min_val is not None else np.min(discretised_values)
            if max_val is not None:
                g_max = max_val
            elif min_val is not None:
                g_max = min_val + np.max(discretised_values) * bin_width
            else:
                g_max = np.max(discretised_values)

            if bin_width > 0:
                num_steps = int(np.round((g_max - g_min) / bin_width))
                if min_val is not None:
                    # Candidates are bin centers
                    idx = np.arange(num_steps, dtype=np.float64)
                    candidates = g_min + (idx + 0.5) * bin_width
                else:
                    idx = np.arange(num_steps + 1, dtype=np.float64)
                    candidates = g_min + idx * bin_width
            else:
                candidates = sorted_vals.astype(np.float64)
        else:
            candidates = sorted_vals

        target_count = int(np.floor(vol_frac * N))

        # Binary search
        low = 0
        high = len(candidates) - 1
        ans_idx = -1

        while low <= high:
            mid = (low + high) // 2
            val = candidates[mid]

            # Convert physical value to index if in discrete mode
            if bin_width is not None and min_val is not None and bin_width > 0:
                check_val = np.floor((val - min_val) / bin_width) + 1
            else:
                check_val = val

            idx = np.searchsorted(sorted_vals, check_val, side="left")
            count = N - idx

            if count <= target_count:
                ans_idx = mid
                high = mid - 1
            else:
                low = mid + 1

        if ans_idx != -1:
            return float(candidates[ans_idx])
        else:
            return float(candidates[-1])

    features["intensity_at_volume_fraction_0.10_GBPN_10"] = (
        get_intensity_at_volume_fraction(0.10)
    )
    features["intensity_at_volume_fraction_0.90_GBPN_90"] = (
        get_intensity_at_volume_fraction(0.90)
    )

    features[
        "intensity_fraction_difference_between_volume_0.10_and_0.90_fractions_CNV2"
    ] = float(
        features["intensity_at_volume_fraction_0.10_GBPN_10"]
        - features["intensity_at_volume_fraction_0.90_GBPN_90"]
    )

    # -------------------------------------------------------------------------
    # 3. Area Under the IVH Curve (AUC)
    # -------------------------------------------------------------------------
    # IVH Curve: Volume Fraction (phi) vs Intensity (I)
    # We construct the curve points from the unique values in the data.
    unique_vals = np.unique(sorted_vals)
    if len(unique_vals) == 1:
        # If there is only one discretised intensity, AUC is 0 by definition.
        features["area_under_the_ivh_curve_9CMM"] = 0.0
    else:
        # P(X >= i)
        # For each unique value, calculate fraction >= value

        # If we have physical mapping, map unique_vals to physical intensities
        if bin_width is not None and min_val is not None:
            # Map index to physical center: min_val + (idx - 0.5) * w ?
            # Standard FBS mapping: index k corresponds to [min + (k-1)w, min + kw)
            # center = min_val + (k - 1 + 0.5) * w
            # k is the value in unique_vals
            intensities_arr = (
                min_val + (unique_vals.astype(np.float64) - 0.5) * bin_width
            )
        else:
            intensities_arr = unique_vals.astype(np.float64)

        # Calculate volume fractions
        # searchsorted returns first index where val fits.
        # Since sorted_vals is sorted, all elements >= val are from searchsorted(val) onwards.
        indices = np.searchsorted(sorted_vals, unique_vals, side="left")
        counts = N - indices
        fractions = counts.astype(np.float64) / float(N)

        # Riemann Sum (Trapezoidal)
        # Integrate fraction(I) over I.
        auc = 0.0
        for k in range(1, len(intensities_arr)):
            i_curr = intensities_arr[k]
            i_prev = intensities_arr[k - 1]
            phi_curr = fractions[k]
            phi_prev = fractions[k - 1]

            # Trapezoid area
            width = i_curr - i_prev
            avg_height = (phi_curr + phi_prev) * 0.5
            auc += width * avg_height

        features["area_under_the_ivh_curve_9CMM"] = float(auc)

    return features

calculate_local_intensity_features(image, mask, *, enabled=True)

Calculate local intensity features: Local and Global Intensity Peak (IBSI 4.5).

Computes intensity peaks using a 1 cm³ spherical neighborhood (radius ~6.2mm). Global peak is the maximum local mean, local peak is the local mean at max intensity.

Parameters:

Name	Type	Description	Default
image	Image	Image object containing intensity data.	required
mask	Image	Image object containing the ROI mask.	required
enabled	bool	If False, returns empty dict immediately (for performance).	True

Returns:

Type	Description
dict[str, float]	Dictionary with 'global_intensity_peak_0F91' and 'local_intensity_peak_VJGA'.

Source code in pictologics/features/intensity.py

def calculate_local_intensity_features(
    image: Image,
    mask: Image,
    *,
    enabled: bool = True,
) -> dict[str, float]:
    """
    Calculate local intensity features: Local and Global Intensity Peak (IBSI 4.5).

    Computes intensity peaks using a 1 cm³ spherical neighborhood (radius ~6.2mm).
    Global peak is the maximum local mean, local peak is the local mean at max intensity.

    Args:
        image: Image object containing intensity data.
        mask: Image object containing the ROI mask.
        enabled: If False, returns empty dict immediately (for performance).

    Returns:
        Dictionary with 'global_intensity_peak_0F91' and 'local_intensity_peak_VJGA'.
    """
    if not enabled:
        return {}

    features: dict[str, float] = {}

    mask_array = mask.array
    data = image.array
    spacing_tuple = (
        float(image.spacing[0]),
        float(image.spacing[1]),
        float(image.spacing[2]),
    )

    # Radius for 1 cm^3 sphere
    radius_mm = 6.2035

    # Get ROI indices
    x_idx, y_idx, z_idx = np.where(mask_array > 0)
    if len(x_idx) == 0:
        return features

    mask_indices = np.ascontiguousarray(
        np.stack([x_idx, y_idx, z_idx], axis=1).astype(np.int32)
    )
    offsets = _sphere_offsets_for_radius(spacing_tuple, radius_mm)

    # Calculate local means only for ROI voxels
    roi_means = _calculate_local_mean_numba(data, mask_indices, offsets)

    # Compute both peaks without allocating ROI intensity arrays.
    global_peak, local_peak = _calculate_local_peaks_numba(
        data, mask_indices, roi_means
    )
    features["global_intensity_peak_0F91"] = float(global_peak)
    features["local_intensity_peak_VJGA"] = float(local_peak)

    return features

calculate_spatial_intensity_features(image, mask, *, enabled=True)

Calculate spatial intensity features: Moran's I and Geary's C (IBSI 4.4).

These features measure spatial autocorrelation of intensity values within the ROI. Computationally intensive (O(N²) where N = number of ROI voxels).

Parameters:

Name	Type	Description	Default
image	Image	Image object containing intensity data.	required
mask	Image	Image object containing the ROI mask.	required
enabled	bool	If False, returns empty dict immediately (for performance).	True

Returns:

Type	Description
dict[str, float]	Dictionary with 'morans_i_index_N365' and 'gearys_c_measure_NPT7'.
dict[str, float]	Returns NaN values if ROI has fewer than 2 voxels or constant intensity.

Source code in pictologics/features/intensity.py

def calculate_spatial_intensity_features(
    image: Image,
    mask: Image,
    *,
    enabled: bool = True,
) -> dict[str, float]:
    """
    Calculate spatial intensity features: Moran's I and Geary's C (IBSI 4.4).

    These features measure spatial autocorrelation of intensity values within the ROI.
    Computationally intensive (O(N²) where N = number of ROI voxels).

    Args:
        image: Image object containing intensity data.
        mask: Image object containing the ROI mask.
        enabled: If False, returns empty dict immediately (for performance).

    Returns:
        Dictionary with 'morans_i_index_N365' and 'gearys_c_measure_NPT7'.
        Returns NaN values if ROI has fewer than 2 voxels or constant intensity.
    """
    if not enabled:
        return {}

    features: dict[str, float] = {}

    mask_array = mask.array
    data = image.array
    sx, sy, sz = (
        float(image.spacing[0]),
        float(image.spacing[1]),
        float(image.spacing[2]),
    )

    # Get ROI indices (X, Y, Z)
    x_idx, y_idx, z_idx = np.where(mask_array > 0)

    if len(x_idx) < 2:
        features["morans_i_index_N365"] = np.nan
        features["gearys_c_measure_NPT7"] = np.nan
        return features

    xi = np.ascontiguousarray(x_idx.astype(np.int32))
    yi = np.ascontiguousarray(y_idx.astype(np.int32))
    zi = np.ascontiguousarray(z_idx.astype(np.int32))

    intensities = np.ascontiguousarray(data[mask_array > 0].astype(np.float64))

    N = len(intensities)
    mean_int = np.mean(intensities)

    # Calculate terms using Parallelized Numba Function
    numer_moran, numer_geary_1, numer_geary_2, sum_weights = (
        _calculate_spatial_features_numba(
            xi, yi, zi, intensities, float(mean_int), sx, sy, sz
        )
    )

    # Moran's I - N365
    denom = _sum_sq_centered(intensities, float(mean_int))

    if denom != 0 and sum_weights != 0:
        moran_i = (N / sum_weights) * (numer_moran / denom)
        features["morans_i_index_N365"] = float(moran_i)
    else:
        features["morans_i_index_N365"] = np.nan

    # Geary's C - NPT7
    if denom != 0 and sum_weights != 0:
        numer = 2 * numer_geary_1 - 2 * numer_geary_2
        geary_c = ((N - 1) / (2 * sum_weights)) * (numer / denom)
        features["gearys_c_measure_NPT7"] = float(geary_c)
    else:
        features["gearys_c_measure_NPT7"] = np.nan

    return features