feat(accelerator): numba zernike + radial_zernikes backend

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,32 @@
 _3D_FEATURES = ("intensity", "sizeshape", "texture", "granularity")
 
 
+def _numba_registries() -> dict[str, dict[str, Callable]]:
+    """Registries for the 'numba' accelerator.
+
+    Composes the numba implementations (``intensity``, ``zernike``,
+    ``radial_zernikes``) 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,
+        get_radial_zernikes as _numba_radial_zernikes,
+        get_zernike as _numba_zernike,
+    )
+
+    return {
+        "core": {
+            **_CORE,
+            "intensity": _numba_intensity,
+            "zernike": _numba_zernike,
+            "radial_zernikes": _numba_radial_zernikes,
+        },
+        "correlation": _CORRELATION,
+    }
+
+
 def _dispatch(name: str) -> dict[str, Callable]:
     from cp_measure import _ACCELERATOR
 
@@ -55,9 +81,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,22 @@
+"""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 accelerates ``intensity``, ``zernike`` and ``radial_zernikes``; the
+global "numba" accelerator composes them with the numpy implementations of every
+other feature (see ``cp_measure.bulk``).
+"""
+
+from cp_measure.core.numba.measureobjectintensity import get_intensity
+from cp_measure.core.numba.measureobjectintensitydistribution import (
+    get_radial_zernikes,
+)
+from cp_measure.core.numba.measureobjectsizeshape import get_zernike
+
+__all__ = [
+    "get_intensity",
+    "get_radial_zernikes",
+    "get_zernike",
+]

src/cp_measure/core/numba/_zernike.py

@@ -0,0 +1,140 @@
+"""Numba Zernike-moment kernel (fused basis-eval + weighted-complex segment-sum).
+
+Both Zernike features — shape (`get_zernike`, weight=1) and intensity-weighted
+(`get_radial_zernikes`, weight=pixel) — reduce to the same per-object
+weighted-complex sum of the Zernike polynomial basis ``V_nm`` over the object's
+pixels. This module fuses the basis evaluation and the segment-sum into one
+single-pass ``@njit`` kernel (no ``(M, K)`` intermediate).
+
+The polynomial *coefficients* come from centrosome (``construct_zernike_lookuptable``,
+degree-only, cheap) via :func:`zernike_coeffs`; the per-pixel evaluation replicates
+``centrosome.zernike.construct_zernike_polynomials`` exactly:
+
+- radial part ``R_nm(r)/r^m`` via Horner over ``r² = xm²+ym²`` using the LUT row,
+  looping exactly ``(n-m)//2 + 1`` terms (forward order, LUT[k,0] = highest power),
+- azimuthal factor ``z^m`` with ``z = ym + 1j*xm`` (ym = normalised row offset,
+  xm = normalised column offset),
+- strict unit-disk cutoff: zero where ``r² > 1`` (keep ``== 1``).
+
+``minimum_enclosing_circle`` (center/radius) and label enumeration stay on the host.
+"""
+
+import centrosome.cpmorphology
+import centrosome.zernike
+import numpy as np
+from numba import njit
+from numpy.typing import NDArray
+
+from cp_measure.primitives.segment import label_to_idx_lut
+
+
+def zernike_coeffs(
+    zernike_indexes: NDArray[np.integer],
+) -> tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]:
+    """Repack centrosome's radial LUT for the Horner kernel.
+
+    Returns ``(lut, nterms, m_arr)`` where ``lut[k, :]`` are the radial-polynomial
+    coefficients of index ``k`` (highest power first), ``nterms[k] = (n-m)//2 + 1``
+    is exactly how many of them to consume in the Horner loop (trailing LUT entries
+    are zero and must NOT be iterated — they would spuriously multiply by r²), and
+    ``m_arr[k]`` is the azimuthal order ``m``.
+    """
+    lut = centrosome.zernike.construct_zernike_lookuptable(zernike_indexes).astype(
+        np.float64
+    )
+    n = zernike_indexes[:, 0].astype(np.int64)
+    m_arr = zernike_indexes[:, 1].astype(np.int64)
+    nterms = ((n - m_arr) // 2 + 1).astype(np.int64)
+    return lut, nterms, m_arr
+
+
+@njit(cache=True)
+def zernike_moments(weights, xm, ym, seg0, lut, nterms, m_arr, n):
+    """Fused per-object weighted-complex Zernike sum (single pass, no (M,K) buffer).
+
+    For each pixel: evaluate every polynomial's ``V_nm`` from ``lut``/``m_arr``
+    (Horner over r² + ``z^m``, ``z = ym + 1j*xm``; skip the pixel when ``r² > 1``)
+    and scatter-add ``weight·Re/Im`` into ``vr``/``vi`` ``(n, K)``. ``seg0[i]`` is the
+    0-based object index of pixel ``i`` (``< 0`` ⇒ skip). Returns ``(vr, vi)``.
+    """
+    K = lut.shape[0]
+    M = xm.shape[0]
+    vr = np.zeros((n, K))
+    vi = np.zeros((n, K))
+    max_m = 0
+    for k in range(K):
+        if m_arr[k] > max_m:
+            max_m = m_arr[k]
+    zr = np.empty(max_m + 1)
+    zi = np.empty(max_m + 1)
+    for i in range(M):
+        seg = seg0[i]
+        if seg < 0:
+            continue
+        x = xm[i]
+        y = ym[i]
+        r2 = x * x + y * y
+        if r2 > 1.0:  # whole basis is 0 outside the unit disk (strict cutoff)
+            continue
+        w = weights[i]
+        zr[0] = 1.0
+        zi[0] = 0.0
+        for mm in range(1, max_m + 1):  # z^mm = z^(mm-1) * (y + i x)
+            zr[mm] = zr[mm - 1] * y - zi[mm - 1] * x
+            zi[mm] = zr[mm - 1] * x + zi[mm - 1] * y
+        for k in range(K):
+            s = 0.0
+            nt = nterms[k]
+            for t in range(nt):  # Horner: R_nm(r)/r^m
+                s = s * r2 + lut[k, t]
+            m = m_arr[k]
+            vr[seg, k] += w * (s * zr[m])
+            vi[seg, k] += w * (s * zi[m])
+    return vr, vi
+
+
+def zernike_moments_per_object(labels, pixels, coeffs):
+    """Per-object weighted-complex Zernike sum shared by both 2D backends.
+
+    Does the part identical to shape and radial Zernikes: enumerate objects,
+    find the per-object minimum enclosing circle, normalise each labeled pixel's
+    coordinates into the unit disk (``ym`` from rows, ``xm`` from cols), and run the
+    fused kernel. ``pixels`` is the intensity image for the weighted (radial) case
+    or ``None`` for the unweighted (shape) case. ``coeffs`` is :func:`zernike_coeffs`.
+
+    Returns ``(vr, vi, radii, seg0)``: ``vr``/``vi`` are ``(n, K)`` Re/Im sums,
+    ``radii`` ``(n,)`` the enclosing-circle radii (shape's ``π·r²`` denominator),
+    ``seg0`` ``(M,)`` the 0-based object index per pixel (radial's pixel-count
+    denominator via ``bincount``). ``n == vr.shape[0]`` is 0 when there are no objects.
+    """
+    lut, nterms, m_arr = coeffs
+    K = lut.shape[0]
+    label_lut, n = label_to_idx_lut(labels)
+    if n == 0:
+        z = np.zeros((0, K))
+        return z, z, np.zeros(0), np.zeros(0, np.int64)
+
+    unique_labels = np.flatnonzero(label_lut >= 0)  # present labels, ascending
+    centers, radii = centrosome.cpmorphology.minimum_enclosing_circle(
+        labels, unique_labels
+    )
+    rows, cols = np.nonzero(labels)
+    seg0 = label_lut[labels[rows, cols]]  # int64 (label_lut is int64), 0-based rank
+    ym = (rows - centers[seg0, 0]) / radii[seg0]  # normalised row offset
+    xm = (cols - centers[seg0, 1]) / radii[seg0]  # normalised column offset
+    weights = (
+        np.ones(seg0.shape[0])
+        if pixels is None
+        else pixels[rows, cols].astype(np.float64)
+    )
+    vr, vi = zernike_moments(
+        weights,
+        np.ascontiguousarray(xm),
+        np.ascontiguousarray(ym),
+        seg0,
+        lut,
+        nterms,
+        m_arr,
+        n,
+    )
+    return vr, vi, radii, seg0