refactor(numba): migrate intensity backend onto (B,Z,Y,X)

noxfile.py

@@ -14,7 +14,7 @@ def tests(session: nox.Session) -> None:
             "pytest-cov",
             "pytest-markdown-docs",
             "six",  # centrosome runtime dep, not declared in its metadata
-            ".",
+            ".[numba]",  # exercise the numba backend + correctness harness in CI
         )
     except Exception:
         session.skip(

pyproject.toml

@@ -17,6 +17,9 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+numba = [
+    "numba>=0.59",
+]
 test = [
     "pytest>=8.4.2",
     "pytest-cov",

src/cp_measure/_detect.py

@@ -0,0 +1,14 @@
+"""Backend capability detection.
+
+Detected ONCE at import via ``importlib.util.find_spec`` — availability is
+checked without importing the package or catching ImportErrors. Dispatch reads
+the flag; the resolved backend path is then called directly and unguarded
+(a backend that is flagged present but raises is a real bug and must surface,
+not be papered over by a try/except fallback).
+
+Other backends (jax, ...) add their own flags here as they are wired.
+"""
+
+import importlib.util
+
+HAS_NUMBA: bool = importlib.util.find_spec("numba") is not None

src/cp_measure/bulk.py

@@ -45,6 +45,23 @@
 _3D_FEATURES = ("intensity", "sizeshape", "texture", "granularity")
 
 
+def _numba_registries() -> dict[str, dict[str, Callable]]:
+    """Registries for the 'numba' accelerator.
+
+    Composes the numba implementations (currently ``intensity`` only) with the
+    numpy implementations of every other feature — a single global "numba"
+    selection still yields a full, working feature set, accelerated where a
+    numba backend exists. This is explicit per-function composition, NOT an
+    error-driven fallback.
+    """
+    from cp_measure.core.numba import get_intensity as _numba_intensity
+
+    return {
+        "core": {**_CORE, "intensity": _numba_intensity},
+        "correlation": _CORRELATION,
+    }
+
+
 def _dispatch(name: str) -> dict[str, Callable]:
     from cp_measure import _ACCELERATOR
 
@@ -55,9 +72,14 @@ def _dispatch(name: str) -> dict[str, Callable]:
             f"'jax' accelerator not yet wired for {name} measurements"
         )
     if _ACCELERATOR == "numba":
-        raise NotImplementedError(
-            f"'numba' accelerator not yet wired for {name} measurements"
-        )
+        from cp_measure._detect import HAS_NUMBA
+
+        if not HAS_NUMBA:
+            raise RuntimeError(
+                "accelerator 'numba' selected but numba is not installed; "
+                "you can install it via `pip install cp_measure[numba]`"
+            )
+        return _numba_registries()[name]
     if _ACCELERATOR == "fastest":
         raise NotImplementedError("'fastest' logic not yet implemented")
     raise ValueError(

src/cp_measure/core/numba/__init__.py

@@ -0,0 +1,14 @@
+"""Numba-accelerated backend.
+
+Selected explicitly by import (``from cp_measure.core.numba import get_intensity``)
+or globally via ``cp_measure.set_accelerator("numba")``. Requires the optional
+``numba`` extra; availability is gated by ``cp_measure._detect.HAS_NUMBA``.
+
+This backend currently accelerates ``intensity`` only; the global "numba"
+accelerator composes it with the numpy implementations of every other feature
+(see ``cp_measure.bulk``).
+"""
+
+from cp_measure.core.numba.measureobjectintensity import get_intensity
+
+__all__ = ["get_intensity"]

src/cp_measure/core/numba/measureobjectintensity.py

@@ -0,0 +1,219 @@
+"""Numba-backed MeasureObjectIntensity.
+
+Drop-in for :func:`cp_measure.core.measureobjectintensity.get_intensity`,
+producing the identical dict of features for 2D and 3D input. The per-label
+reductions — including ``Location_MaxIntensity_*`` — run entirely as fused
+single-pass numba kernels over a flat segment representation
+(:mod:`cp_measure.primitives`). The max position uses a deterministic
+``>=``-last rule that is bit-identical to ``scipy.ndimage.maximum_position`` on
+real (tie-free) data; only exact-value ties can differ (scipy's tie pick is
+quicksort-dependent and not stable across numpy versions).
+"""
+
+import numpy
+import skimage.segmentation
+from numpy.typing import NDArray
+
+from cp_measure.core.measureobjectintensity import (
+    C_LOCATION,
+    INTEGRATED_INTENSITY,
+    INTEGRATED_INTENSITY_EDGE,
+    INTENSITY,
+    LOC_CMI_X,
+    LOC_CMI_Y,
+    LOC_CMI_Z,
+    LOC_MAX_X,
+    LOC_MAX_Y,
+    LOC_MAX_Z,
+    LOWER_QUARTILE_INTENSITY,
+    MAD_INTENSITY,
+    MASS_DISPLACEMENT,
+    MAX_INTENSITY,
+    MAX_INTENSITY_EDGE,
+    MEAN_INTENSITY,
+    MEAN_INTENSITY_EDGE,
+    MEDIAN_INTENSITY,
+    MIN_INTENSITY,
+    MIN_INTENSITY_EDGE,
+    STD_INTENSITY,
+    STD_INTENSITY_EDGE,
+    UPPER_QUARTILE_INTENSITY,
+)
+from cp_measure.primitives.segment import label_to_idx_lut
+from cp_measure.primitives.shapes import to_bzyx
+from cp_measure.primitives._segment_numba import (
+    flatten_numba,
+    inner_boundary,
+    segment_moments,
+    segment_quantiles,
+    segment_resid_sumsq,
+    segment_stats,
+)
+
+
+def get_intensity(
+    masks,
+    pixels,
+    edge_measurements: bool = True,
+):
+    """Per-object intensity features; accepts a single image/volume or a batch.
+
+    Returns a single feature dict for a lone 2D image or 3D volume, or a list of
+    such dicts for a batch (4D ``(B, Z, Y, X)`` array or a list of images). See
+    :func:`cp_measure.primitives.shapes.to_bzyx` for the accepted shapes; a single
+    image is the ``B == 1`` case of one batched code path.
+    """
+    masks_zyx, pixels_zyx, unwrap = to_bzyx(masks, pixels)
+    results = [
+        _intensity_volume(m, p, edge_measurements)
+        for m, p in zip(masks_zyx, pixels_zyx)
+    ]
+    return unwrap(results)
+
+
+def _intensity_volume(
+    masks: NDArray[numpy.integer],
+    pixels: NDArray[numpy.floating],
+    edge_measurements: bool = True,
+) -> dict[str, NDArray[numpy.floating]]:
+    """Intensity features for one ``(Z, Y, X)`` volume (a 2D image is ``Z == 1``).
+
+    ``masks`` is a labeled array where 0 is background. A ``(1, Y, X)`` mask paired
+    with a ``(Z, Y, X)`` image is the 2D-mask-on-stack broadcast (preserved from the
+    original backend). The MAD fraction follows the input dimensionality — 1/2 for a
+    2D image (``Z == 1``), 1/3 for a true volume — matching the numpy baseline's
+    ``1 / pixels.ndim``.
+    """
+    orig_ndim = 2 if pixels.shape[0] == 1 else 3
+
+    lut, nobjects = label_to_idx_lut(masks)
+
+    integrated_intensity = numpy.zeros(nobjects)
+    mean_intensity = numpy.zeros(nobjects)
+    std_intensity = numpy.zeros(nobjects)
+    min_intensity = numpy.zeros(nobjects)
+    max_intensity = numpy.zeros(nobjects)
+    mass_displacement = numpy.zeros(nobjects)
+    lower_quartile_intensity = numpy.zeros(nobjects)
+    median_intensity = numpy.zeros(nobjects)
+    mad_intensity = numpy.zeros(nobjects)
+    upper_quartile_intensity = numpy.zeros(nobjects)
+    cmi_x = numpy.zeros(nobjects)
+    cmi_y = numpy.zeros(nobjects)
+    cmi_z = numpy.zeros(nobjects)
+    max_x = numpy.zeros(nobjects)
+    max_y = numpy.zeros(nobjects)
+    max_z = numpy.zeros(nobjects)
+
+    values, seg0, xc, yc, zc = flatten_numba(
+        numpy.ascontiguousarray(masks),
+        numpy.ascontiguousarray(pixels),
+        lut,
+    )
+    has_objects = values.size > 0
+
+    if has_objects:
+        (
+            count,
+            sumI,
+            minI,
+            maxI,
+            max_x,
+            max_y,
+            max_z,
+            sx,
+            sy,
+            sz,
+            sxI,
+            syI,
+            szI,
+        ) = segment_moments(values, seg0, xc, yc, zc, nobjects)
+        cnt = count.astype(numpy.float64)
+        with numpy.errstate(invalid="ignore", divide="ignore"):
+            integrated_intensity = sumI
+            mean_intensity = sumI / cnt
+            ss = segment_resid_sumsq(values, seg0, nobjects, mean_intensity)
+            std_intensity = numpy.sqrt(ss / cnt)
+            min_intensity = minI
+            max_intensity = maxI
+
+            cm_x = sx / cnt
+            cm_y = sy / cnt
+            cm_z = sz / cnt
+            cmi_x = sxI / sumI
+            cmi_y = syI / sumI
+            cmi_z = szI / sumI
+            mass_displacement = numpy.sqrt(
+                (cm_x - cmi_x) ** 2 + (cm_y - cmi_y) ** 2 + (cm_z - cmi_z) ** 2
+            )
+
+        (
+            lower_quartile_intensity,
+            median_intensity,
+            upper_quartile_intensity,
+            mad_intensity,
+        ) = segment_quantiles(values, seg0, count, nobjects, 1.0 / orig_ndim)
+
+    if edge_measurements:
+        integrated_intensity_edge = numpy.zeros(nobjects)
+        mean_intensity_edge = numpy.zeros(nobjects)
+        std_intensity_edge = numpy.zeros(nobjects)
+        min_intensity_edge = numpy.zeros(nobjects)
+        max_intensity_edge = numpy.zeros(nobjects)
+
+        # 2D plane (Z==1): numba inner-boundary kernel, bit-identical to skimage
+        # mode="inner" but ~12-27x faster. True 3D keeps skimage (6-neighbourhood).
+        if masks.shape[0] == 1:
+            emask = inner_boundary(numpy.ascontiguousarray(masks[0]))[numpy.newaxis] > 0
+        else:
+            emask = skimage.segmentation.find_boundaries(masks, mode="inner") > 0
+        e_values = pixels[emask].astype(numpy.float64)
+        e_seg0 = lut[masks[emask]]
+
+        if e_values.size > 0:
+            ecount, esum, emin, emax = segment_stats(e_values, e_seg0, nobjects)
+            edge_obj = ecount > 0
+            ecnt = ecount.astype(numpy.float64)
+            with numpy.errstate(invalid="ignore", divide="ignore"):
+                emean = esum / ecnt
+                ess = segment_resid_sumsq(e_values, e_seg0, nobjects, emean)
+                estd = numpy.sqrt(ess / ecnt)
+            integrated_intensity_edge[edge_obj] = esum[edge_obj]
+            mean_intensity_edge[edge_obj] = emean[edge_obj]
+            std_intensity_edge[edge_obj] = estd[edge_obj]
+            min_intensity_edge[edge_obj] = emin[edge_obj]
+            max_intensity_edge[edge_obj] = emax[edge_obj]
+
+    measurement_names = [
+        (INTENSITY, INTEGRATED_INTENSITY, integrated_intensity),
+        (INTENSITY, MEAN_INTENSITY, mean_intensity),
+        (INTENSITY, STD_INTENSITY, std_intensity),
+        (INTENSITY, MIN_INTENSITY, min_intensity),
+        (INTENSITY, MAX_INTENSITY, max_intensity),
+        (INTENSITY, MASS_DISPLACEMENT, mass_displacement),
+        (INTENSITY, LOWER_QUARTILE_INTENSITY, lower_quartile_intensity),
+        (INTENSITY, MEDIAN_INTENSITY, median_intensity),
+        (INTENSITY, MAD_INTENSITY, mad_intensity),
+        (INTENSITY, UPPER_QUARTILE_INTENSITY, upper_quartile_intensity),
+        (C_LOCATION, LOC_CMI_X, cmi_x),
+        (C_LOCATION, LOC_CMI_Y, cmi_y),
+        (C_LOCATION, LOC_CMI_Z, cmi_z),
+        (C_LOCATION, LOC_MAX_X, max_x),
+        (C_LOCATION, LOC_MAX_Y, max_y),
+        (C_LOCATION, LOC_MAX_Z, max_z),
+    ]
+    if edge_measurements:
+        measurement_names.extend(
+            [
+                (INTENSITY, INTEGRATED_INTENSITY_EDGE, integrated_intensity_edge),
+                (INTENSITY, MEAN_INTENSITY_EDGE, mean_intensity_edge),
+                (INTENSITY, STD_INTENSITY_EDGE, std_intensity_edge),
+                (INTENSITY, MIN_INTENSITY_EDGE, min_intensity_edge),
+                (INTENSITY, MAX_INTENSITY_EDGE, max_intensity_edge),
+            ]
+        )
+
+    return {
+        "{}_{}".format(category, feature_name): measurement
+        for category, feature_name, measurement in measurement_names
+    }

src/cp_measure/primitives/__init__.py

@@ -0,0 +1,17 @@
+"""Shared primitive layer.
+
+Backend-agnostic building blocks that the per-backend feature implementations
+(``cp_measure.core`` = numpy, ``cp_measure.core.numba`` = numba, ...) compose.
+
+A labeled image is reduced to a flat *segment* representation (values + 0-based
+segment index + per-axis coordinates) which the segment kernels reduce. All
+spatial structure (2D vs 3D) and any future batch/image axis are encoded in that
+flat segment index, so a single set of segment kernels covers every case without
+a rewrite. The numpy ``label_to_idx_lut`` (in ``segment``) builds the
+label→index lookup; the flattening + reductions themselves live in
+``_segment_numba``.
+
+This is an internal layer with no curated public API: import its building
+blocks directly from the concrete modules (``primitives.segment``,
+``primitives._segment_numba``) rather than re-exporting them here.
+"""