Intensity Features API

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).
            The sums run in float64, also for float32 values.

    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 {}
    # float32 values (for example filter responses) would sum in float32: with its rounding
    # errors, and in an order that depends on the numpy version. A C-order copy of a strided
    # array: the sums then do not depend on the memory layout.
    values = np.ascontiguousarray(values, dtype=np.float64)

    features: dict[str, float] = {}

    # 4.1.1 Mean intensity (Q4LE). One pass also gives the energy (4.1.17), the minimum and
    # the maximum.
    mean_val, energy, min_val, max_val = _first_pass(values)
    features["mean_intensity_Q4LE"] = float(mean_val)

    # 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. One call gives all five.
    p10, p25, p75, p90, median_val = _order_statistics(values, (min_val, max_val))

    # 4.1.2 Intensity variance (ECT3): the second central moment. The pass that gives it
    # also gives the third and the fourth moments and the sums of the three deviations
    # (4.1.12 to 4.1.14). Equal finite values have no spread: IBSI gives 0 for all six,
    # where the sums give rounding noise (a variance of about 1e-34 for 0.1 values, whose
    # mean is not exact).
    if min_val == max_val and math.isfinite(min_val):
        m2 = m3 = m4 = mad = robust_mad = median_ad = 0.0
    else:
        m2, m3, m4, mad, robust_mad, median_ad = _moments_and_deviations(
            values, float(mean_val), float(median_val), float(p10), float(p90)
        )
    var_val = m2
    features["intensity_variance_ECT3"] = float(var_val)

    # 4.1.3 Intensity skewness (KE2A) and 4.1.4 kurtosis (IPH6). IBSI defines both as 0
    # when the variance is 0. Equal values count as 0 variance here, also infinite ones
    # (their variance is NaN).
    if var_val == 0.0 or min_val == max_val:
        features["intensity_skewness_KE2A"] = 0.0
        features["intensity_kurtosis_IPH6"] = 0.0
    else:
        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

    # 4.1.5 Median intensity (Y12H)
    features["median_intensity_Y12H"] = float(median_val)

    # 4.1.6 Minimum intensity (1GSF)
    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)
    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(mad)

    # 4.1.13 Robust mean absolute deviation (1128)
    features["intensity_robust_mean_absolute_deviation_1128"] = float(robust_mad)

    # 4.1.14 Median absolute deviation (N72L)
    features["intensity_median_absolute_deviation_N72L"] = float(median_ad)

    # 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): a fixed order, the same bits in every buffer
    energy = float(energy)
    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

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

    disc = np.asarray(discretised_values)
    if disc.ndim != 1:
        raise ValueError(f"discretised_values must be a 1-D array; got shape {disc.shape}")
    n = disc.size
    # Integer values (but uint64, which bincount does not read) take every feature from
    # their bin counts, and so do float64 values that are whole numbers in [1, n_bins]:
    # the levels themselves. Other values are read by the kernels.
    counted = disc.dtype.kind in "iu" and disc.dtype != np.uint64

    counts: Optional[npt.NDArray[np.int64]] = None
    if n_bins is not None and (counted or disc.dtype == np.float64):
        # IBSI: histogram over the full discretisation range [1, N_g]
        counts, low, high, whole = _level_counts(disc, n_bins)
        if not whole:
            counts = None
        elif low < 1 or high > n_bins:
            raise ValueError(_range_error(n_bins, int(low), int(high)))
        else:
            counts, hist_origin = counts[1:], 1
    elif counted:
        # Observed value range; shifting also supports negative values for bincount
        hist_origin = int(np.min(disc))
        counts = np.bincount(
            disc.astype(np.int64) - hist_origin, minlength=int(np.max(disc)) - hist_origin + 1
        )
    if counts is not None:
        out = np.empty(len(_HISTOGRAM_FEATURES))
        _counted_statistics_numba(counts, hist_origin, out)
        return dict(zip(_HISTOGRAM_FEATURES, out.tolist(), strict=True))

    # One float64 copy of another type (float32, float16, uint64, bool): the mean, the
    # variance, the moments, the percentiles and the deviations then sum in float64. The
    # copy keeps every NaN, infinite and non-whole value.
    if disc.dtype != np.float64:
        disc = disc.astype(np.float64)
    # Discretised values are whole numbers (a float array holds them too); a value that is
    # not whole, NaN or infinite is a caller error (a pipeline passes int32). The check
    # reads the values once.
    if not _all_whole(np.ascontiguousarray(disc)):
        raise ValueError(
            "discretised_values must be finite whole numbers (bin numbers); a value "
            "that is not whole, NaN or infinite was given."
        )
    features: dict[str, float] = {}
    min_val_i = int(np.min(disc))
    max_val_i = int(np.max(disc))
    if n_bins is not None:
        if min_val_i < 1 or max_val_i > n_bins:
            raise ValueError(_range_error(n_bins, 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) and 4.2.2 variance (CH89)
    mean_disc = float(np.mean(disc))
    var_disc = float(np.var(disc, ddof=0))
    features["mean_discretised_intensity_X6K6"] = float(mean_disc)
    features["discretised_intensity_variance_CH89"] = float(var_disc)

    # 4.2.3 Discretised intensity skewness (88K1) and 4.2.4 kurtosis (C3I7): 0 when the
    # variance is 0 (IBSI)
    if var_disc == 0.0:
        features["discretised_intensity_skewness_88K1"] = 0.0
        features["discretised_intensity_kurtosis_C3I7"] = 0.0
    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_disc==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, median_val = _order_statistics(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"]
    )

    mad = _mean_abs_dev(disc, float(mean_disc))
    robust_mad = _robust_mean_abs_dev(disc, float(p10), float(p90))
    median_ad = _mean_abs_dev(disc, float(median_val))
    features["intensity_histogram_mean_absolute_deviation_D2ZX"] = float(mad)
    features["intensity_histogram_robust_mean_absolute_deviation_WRZB"] = float(robust_mad)
    features["intensity_histogram_median_absolute_deviation_4RNL"] = float(median_ad)

    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

    # Entropy (TLU2) and uniformity (BJ5W) over the occupied bins
    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

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.ascontiguousarray(discretised_values)  # one layout for the numba kernels
    count_below, kth, unique_vals, below_unique = _sorted_view(vals)

    # -------------------------------------------------------------------------
    # 1. Volume Fractions
    # -------------------------------------------------------------------------
    # The range G of the intensity fractions (IBSI 4.3): the target range, else the range
    # of the discretisation. An end that is not given is that of the ROI, as for a
    # half-open re-segmentation range; with FBS bins, the upper end is the center of the
    # last bin (with no min_val, a lower end is a value: no bins convert).
    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
    if t_min is None:
        t_min = float(kth(0))
    if t_max is None:
        t_max = float(kth(N - 1))
        if bin_width is not None and min_val is not None:
            t_max = min_val + (t_max - 0.5) * bin_width
    t_range = t_max - t_min

    def get_volume_fraction_at_intensity_fraction(frac: float) -> float:
        threshold = t_min + frac * t_range
        if bin_width is not None and min_val is not None:
            threshold = np.floor((threshold - min_val) / bin_width) + 1  # the FBS bin
        return float((N - count_below(threshold)) / N)

    features["volume_at_intensity_fraction_0.10_BC2M_10"] = (
        get_volume_fraction_at_intensity_fraction(0.10)
    )
    features["volume_at_intensity_fraction_0.90_BC2M_90"] = (
        get_volume_fraction_at_intensity_fraction(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(vals.dtype, np.integer)
        ):
            target_count = int(np.floor(vol_frac * N))
            if target_count <= 0:
                return float(kth(N - 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(kth(k - 1))
            t = v + 1
            vmax = int(kth(N - 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
                    candidates: _Grid | npt.NDArray[Any] = _Grid(num_steps, g_min, bin_width, 0.5)
                else:
                    candidates = _Grid(num_steps + 1, g_min, bin_width, 0.0)
            else:
                # The values that occur: the first one that meets the condition below is
                # the first such value of all the sorted values
                candidates = unique_vals.astype(np.float64)
        else:
            candidates = unique_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

            count = N - count_below(check_val)

            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: the volume fraction (the values >= I) against the intensity I, at the
    # unique values of the data.
    if len(unique_vals) == 1:
        # One discretised intensity: the AUC is 0 by definition.
        features["area_under_the_ivh_curve_9CMM"] = 0.0
    else:
        # With a physical mapping, a level k stands for the centre of its FBS bin
        # [min + (k - 1) w, min + k w): min + (k - 0.5) w
        if bin_width is not None and min_val is not None:
            intensities_arr = min_val + (unique_vals.astype(np.float64) - 0.5) * bin_width
        else:
            intensities_arr = unique_vals.astype(np.float64)

        # The volume fractions: the values >= each unique value
        counts = N - below_unique
        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

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
    offsets = _sphere_offsets_for_radius(spacing_tuple, _LOCAL_PEAK_RADIUS_MM)

    # Crop to the ROI box plus the reach of the sphere. Every sphere neighbour of an ROI
    # voxel that is in the image is also in the crop, so the local means and the peaks
    # stay the same, and the ROI search reads only the crop.
    if mask_array.size >= _LOCAL_CROP_MIN_SIZE:
        bbox = compute_nonzero_bbox(mask_array)
        if bbox is None:
            return features
        reach = np.abs(offsets).max(axis=0)
        crop = tuple(
            slice(max(box.start - int(r), 0), min(box.stop + int(r), size))
            for box, r, size in zip(bbox, reach, data.shape, strict=True)
        )
        data = data[crop]
        mask_array = mask_array[crop]
    # The local kernels read float64 or float32, the types that the import compiles; an image
    # of another type goes as a float64 copy (the same values). A column-order image (a small
    # one, without the crop) goes as a row-order copy, the layout that the import compiles.
    if data.dtype not in (np.float64, np.float32):
        data = data.astype(np.float64, order="C")
    elif data.flags.f_contiguous and not data.flags.c_contiguous:
        data = np.ascontiguousarray(data)

    # The ROI scan reads a row-order bool array (one kernel signature per image type).
    roi = np.ascontiguousarray(mask_array != 0)
    indices, values, runs = _roi_scan(roi, data)
    n = indices.shape[0]
    if n == 0:
        return features

    peaks = None
    if n * offsets.shape[0] >= _TWO_STAGE_MIN_WORK:
        peaks = _local_peaks_two_stage(data, indices, values, runs, offsets, spacing_tuple, roi)
    if peaks is None:
        # The exact sphere mean of every ROI voxel, in threads from _LOCAL_MEAN_PARALLEL_MIN pairs
        exact_mean = sized(
            _calculate_local_mean_numba,
            _calculate_local_mean_numba_serial,
            n * offsets.shape[0],
            _LOCAL_MEAN_PARALLEL_MIN,
        )
        roi_means = exact_mean(data, indices, offsets)
        # Compute both peaks without allocating ROI intensity arrays.
        peaks = _calculate_local_peaks_numba(data, indices, roi_means)
    features["global_intensity_peak_0F91"] = float(peaks[0])
    features["local_intensity_peak_VJGA"] = float(peaks[1])

    return features

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.
    A time model picks the faster of two methods (see _spatial_plan): a loop over all
    voxel pairs, whose time grows with N² (N = number of ROI voxels), or FFT convolutions,
    whose time grows with the ROI box. Both give the same values to about 1e-13
    (relative). The FFT needs about 32 bytes per point of a grid that is twice the ROI box
    along each axis, and it uses at most half of the memory (shared by the workers of
    run_batch). When the grid needs more, the FFT runs on slabs of the box: a few times
    longer, within the memory.

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

    # The ROI (mask != 0) in its bounding box (an empty mask gives an empty box). Its
    # voxels keep their full-image order, so `intensities` is the same as data[mask != 0].
    bbox = compute_nonzero_bbox(mask_array) or (slice(0, 0),) * 3
    roi = mask_array[bbox] != 0
    intensities = np.ascontiguousarray(data[bbox][roi].astype(np.float64, copy=False))

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

    mean_int = np.mean(intensities)
    denom = _sum_sq_centered(intensities, float(mean_int))

    # The FFT grid is at least 2 n - 1 along each axis of the box (of a slab, along its cut
    # axis), so no pair wraps around.
    method, axis, rows, fshape = _spatial_plan(N, roi.shape)
    if method == "fft":
        numer_moran, numer_geary, sum_weights = _spatial_sums_fft(
            roi, intensities - mean_int, (sx, sy, sz), fshape
        )
    elif method != "pairs":
        diff = np.zeros(roi.shape)
        diff[roi] = intensities - mean_int
        numer_moran, numer_geary, sum_weights = _spatial_sums_slabs(
            roi, diff, (sx, sy, sz), axis, rows, fshape
        )
    else:
        # Full-image voxel indices (the box start added back), so the pair distances do
        # not depend on the crop.
        x_idx, y_idx, z_idx = (
            np.ascontiguousarray((idx + box.start).astype(np.int32))
            for idx, box in zip(np.nonzero(roi), bbox, strict=True)
        )
        pairs = sized(
            _calculate_spatial_features_numba,
            _calculate_spatial_features_numba_serial,
            N,
            _PAIRS_PARALLEL_MIN,
        )
        numer_moran, numer_geary, sum_weights = pairs(
            x_idx, y_idx, z_idx, intensities, float(mean_int), sx, sy, sz
        )

    # Moran's I - N365

    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