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)
        else:  # pragma: no cover  # unreachable: denom==0 iff var_val==0, caught above
            features["intensity_skewness_KE2A"] = np.nan
            features["intensity_kurtosis_IPH6"] = np.nan

    # IBSI percentiles use the nearest-rank convention (the smallest value with at
    # least p% of the data at or below it), i.e. numpy's 'inverted_cdf'. This
    # reproduces the IBSI benchmark (e.g. P90 = 4), whereas linear interpolation
    # would give an interpolated 4.2. The median (Y12H) is the conventional
    # sample median, so it is computed separately.
    p10, p25, p75, p90 = np.percentile(values, [10, 25, 75, 90], method="inverted_cdf")
    median_val = np.median(values)

    # 4.1.5 Median intensity (Y12H)
    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)

    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, n_bins=None)

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
n_bins	Optional[int]	Total number of bins N_g used at discretisation. When given, the histogram spans the full IBSI range [1, N_g], including bins that are empty because the data does not reach them (this only affects the gradient features). When None, the histogram spans the observed value range.	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.

Raises:

Type	Description
ValueError	If n_bins is given and the values do not lie in [1, n_bins].

Example

import numpy as np
from pictologics.features.intensity import calculate_intensity_histogram_features

discretised_values = np.array([1, 1, 2, 2, 3, 4, 4, 4], dtype=np.float64)
features = calculate_intensity_histogram_features(discretised_values, n_bins=4)
print(features["mean_discretised_intensity_X6K6"])
# 2.625
print(features["intensity_histogram_mode_AMMC"])
# 4.0

Source code in pictologics/features/intensity.py

def calculate_intensity_histogram_features(
    discretised_values: npt.NDArray[np.floating[Any]],
    n_bins: Optional[int] = None,
) -> 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).
        n_bins: Total number of bins N_g used at discretisation. When given, the
            histogram spans the full IBSI range [1, N_g], including bins that are
            empty because the data does not reach them (this only affects the
            gradient features). When None, the histogram spans the observed value
            range.

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

    Raises:
        ValueError: If n_bins is given and the values do not lie in [1, n_bins].

    Example:
        ```python
        import numpy as np
        from pictologics.features.intensity import calculate_intensity_histogram_features

        discretised_values = np.array([1, 1, 2, 2, 3, 4, 4, 4], dtype=np.float64)
        features = calculate_intensity_histogram_features(discretised_values, n_bins=4)
        print(features["mean_discretised_intensity_X6K6"])
        # 2.625
        print(features["intensity_histogram_mode_AMMC"])
        # 4.0
        ```
    """
    if len(discretised_values) == 0:
        return {}

    features: dict[str, float] = {}

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

    min_val_i = int(np.min(disc))
    max_val_i = int(np.max(disc))

    if n_bins is not None:
        # IBSI: histogram over the full discretisation range [1, N_g].
        if min_val_i < 1 or max_val_i > n_bins:
            raise ValueError(
                f"discretised values must lie in [1, n_bins={n_bins}]; "
                f"got range [{min_val_i}, {max_val_i}]"
            )
        hist_origin = 1
        counts_full = np.bincount(disc.astype(np.int64) - 1, minlength=n_bins)
    else:
        # Observed value range; shifting also supports negative values
        # for bincount compatibility.
        hist_origin = min_val_i
        counts_full = np.bincount(
            disc.astype(np.int64) - min_val_i, 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)
        else:  # pragma: no cover  # unreachable: denom==0 iff var_val==0, caught above
            features["discretised_intensity_skewness_88K1"] = np.nan
            features["discretised_intensity_kurtosis_C3I7"] = np.nan

    # IBSI nearest-rank percentiles (see calculate_intensity_features); the
    # discretised median (WIFQ) is the conventional sample median.
    p10, p25, p75, p90 = np.percentile(disc, [10, 25, 75, 90], method="inverted_cdf")
    median_val = np.median(disc)

    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)

    # IBSI: with multiple modes, select the one closest to the mean discretised
    # intensity; if two modes are equidistant from the mean, select the lower one
    # (argmin returns the first, i.e. lowest, candidate on ties).
    mode_candidates = np.flatnonzero(counts_full == np.max(counts_full)) + hist_origin
    best_mode = int(np.argmin(np.abs(mode_candidates - mean_disc)))
    features["intensity_histogram_mode_AMMC"] = float(mode_candidates[best_mode])

    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(hist_origin + 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(hist_origin + 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.

Raises:

Type	Description
ValueError	If both min_val and max_val are given and max_val < min_val.

Example

import numpy as np
from pictologics.features.intensity import calculate_ivh_features

discretised_values = np.arange(1, 11, dtype=np.float64)
features = calculate_ivh_features(
    discretised_values, bin_width=1.0, min_val=0.0, max_val=10.0
)
print(round(features["area_under_the_ivh_curve_9CMM"], 2))
# 4.95

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.

    Raises:
        ValueError: If both min_val and max_val are given and max_val < min_val.

    Example:
        ```python
        import numpy as np
        from pictologics.features.intensity import calculate_ivh_features

        discretised_values = np.arange(1, 11, dtype=np.float64)
        features = calculate_ivh_features(
            discretised_values, bin_width=1.0, min_val=0.0, max_val=10.0
        )
        print(round(features["area_under_the_ivh_curve_9CMM"], 2))
        # 4.95
        ```
    """
    if len(discretised_values) == 0:
        return {}

    if min_val is not None and max_val is not None and max_val < min_val:
        raise ValueError(f"max_val ({max_val}) must be >= min_val ({min_val})")

    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:
                # Degenerate ranges (g_max within half a bin of g_min) would yield
                # an empty candidate grid; keep at least one candidate.
                num_steps = max(int(np.round((g_max - g_min) / bin_width)), 1)
                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.
    # sorted_vals is already sorted, so unique values are found in O(N) via a
    # neighbour-difference mask instead of np.unique's redundant O(N log N) sort.
    if sorted_vals.size == 0:  # pragma: no cover  # unreachable: empty input returns at func top
        unique_vals = sorted_vals
    else:
        keep = np.empty(sorted_vals.shape, dtype=bool)
        keep[0] = True
        np.not_equal(sorted_vals[1:], sorted_vals[:-1], out=keep[1:])
        unique_vals = sorted_vals[keep]
    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)

        # Trapezoidal integration of fraction(I) over I.
        features["area_under_the_ivh_curve_9CMM"] = float(
            np.trapezoid(fractions, intensities_arr)
        )

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

Example

import numpy as np
from pictologics.loader import Image
from pictologics.features.intensity import calculate_local_intensity_features

data = np.random.default_rng(0).normal(size=(6, 6, 6))
mask_arr = np.zeros((6, 6, 6), dtype=np.uint8)
mask_arr[2:4, 2:4, 2:4] = 1

image = Image(array=data, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))
mask = Image(array=mask_arr, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))

features = calculate_local_intensity_features(image, mask)
print(sorted(features.keys()))
# ['global_intensity_peak_0F91', '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'.

    Example:
        ```python
        import numpy as np
        from pictologics.loader import Image
        from pictologics.features.intensity import calculate_local_intensity_features

        data = np.random.default_rng(0).normal(size=(6, 6, 6))
        mask_arr = np.zeros((6, 6, 6), dtype=np.uint8)
        mask_arr[2:4, 2:4, 2:4] = 1

        image = Image(array=data, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))
        mask = Image(array=mask_arr, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))

        features = calculate_local_intensity_features(image, mask)
        print(sorted(features.keys()))
        # ['global_intensity_peak_0F91', '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.

Example

import numpy as np
from pictologics.loader import Image
from pictologics.features.intensity import calculate_spatial_intensity_features

data = np.random.default_rng(0).normal(size=(6, 6, 6))
mask_arr = np.zeros((6, 6, 6), dtype=np.uint8)
mask_arr[2:4, 2:4, 2:4] = 1

image = Image(array=data, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))
mask = Image(array=mask_arr, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))

features = calculate_spatial_intensity_features(image, mask)
print(sorted(features.keys()))
# ['gearys_c_measure_NPT7', 'morans_i_index_N365']

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.

    Example:
        ```python
        import numpy as np
        from pictologics.loader import Image
        from pictologics.features.intensity import calculate_spatial_intensity_features

        data = np.random.default_rng(0).normal(size=(6, 6, 6))
        mask_arr = np.zeros((6, 6, 6), dtype=np.uint8)
        mask_arr[2:4, 2:4, 2:4] = 1

        image = Image(array=data, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))
        mask = Image(array=mask_arr, spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0))

        features = calculate_spatial_intensity_features(image, mask)
        print(sorted(features.keys()))
        # ['gearys_c_measure_NPT7', 'morans_i_index_N365']
        ```
    """
    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)
    roi = mask_array > 0
    x_idx, y_idx, z_idx = np.where(roi)

    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[roi].astype(np.float64))

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

    # Calculate terms using Parallelized Numba Function
    numer_moran, numer_geary, 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:
        geary_c = ((N - 1) / (2 * sum_weights)) * (numer_geary / denom)
        features["gearys_c_measure_NPT7"] = float(geary_c)
    else:
        features["gearys_c_measure_NPT7"] = np.nan

    return features