diff --git a/noxfile.py b/noxfile.py
index 8f6cae0..025d7c7 100644
--- a/noxfile.py
+++ b/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(
diff --git a/pyproject.toml b/pyproject.toml
index a89d1ed..f7b4e31 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,6 +17,9 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+numba = [
+    "numba>=0.59",
+]
 test = [
     "pytest>=8.4.2",
     "pytest-cov",
diff --git a/src/cp_measure/_detect.py b/src/cp_measure/_detect.py
new file mode 100644
index 0000000..9dd8f9f
--- /dev/null
+++ b/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
diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py
index 1f2dc94..f824940 100644
--- a/src/cp_measure/bulk.py
+++ b/src/cp_measure/bulk.py
@@ -45,6 +45,30 @@
 _3D_FEATURES = ("intensity", "sizeshape", "texture", "granularity")
 
 
+def _numba_registries() -> dict[str, dict[str, Callable]]:
+    """Registries for the 'numba' accelerator.
+
+    Composes the numba implementations (``intensity``, ``granularity``) 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_granularity as _numba_granularity,
+        get_intensity as _numba_intensity,
+    )
+
+    return {
+        "core": {
+            **_CORE,
+            "intensity": _numba_intensity,
+            "granularity": _numba_granularity,
+        },
+        "correlation": _CORRELATION,
+    }
+
+
 def _dispatch(name: str) -> dict[str, Callable]:
     from cp_measure import _ACCELERATOR
 
@@ -55,9 +79,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(
diff --git a/src/cp_measure/core/numba/__init__.py b/src/cp_measure/core/numba/__init__.py
new file mode 100644
index 0000000..65b9beb
--- /dev/null
+++ b/src/cp_measure/core/numba/__init__.py
@@ -0,0 +1,15 @@
+"""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`` and ``granularity``; the global "numba"
+accelerator composes them with the numpy implementations of every other feature
+(see ``cp_measure.bulk``).
+"""
+
+from cp_measure.core.numba.measuregranularity import get_granularity
+from cp_measure.core.numba.measureobjectintensity import get_intensity
+
+__all__ = ["get_granularity", "get_intensity"]
diff --git a/src/cp_measure/core/numba/_granularity.py b/src/cp_measure/core/numba/_granularity.py
new file mode 100644
index 0000000..27b3b22
--- /dev/null
+++ b/src/cp_measure/core/numba/_granularity.py
@@ -0,0 +1,358 @@
+"""Numba 2D morphology kernels for the granularity backend.
+
+Single-threaded ``@njit(cache=True)`` kernels, bit-exact with the skimage
+operations the numpy baseline uses:
+
+- :func:`disk_erosion_2d` / :func:`disk_dilation_2d` — greyscale disk min/max via
+  a **row-decomposed van-Herk/Gil-Werman (VHG)** sliding 1-D min/max over the
+  ``2r+1`` disk rows (``O(r·HW)``), matching ``skimage.morphology.erosion/dilation``
+  with a ``disk(radius)`` footprint (used once for background opening).
+- :func:`erosion_4conn_2d` / :func:`dilation_4conn_2d` — 5-tap (disk(1)/cross)
+  min/max, matching ``skimage`` with a ``disk(1)`` footprint (the per-iteration
+  erosion + the dilation used inside reconstruction).
+- :func:`reconstruction_by_dilation_2d` — morphological reconstruction by dilation
+  (4-connectivity) via the Vincent (1993) FIFO-hybrid, matching
+  ``skimage.morphology.reconstruction(seed, mask, footprint=disk(1))``.
+
+Borders match skimage's footprint ops via ``edge`` (clamp) padding applied on the
+host before the kernel; min/max selection makes the result border-mode-insensitive
+for disks (verified bit-exact vs skimage over reflect/symmetric/edge).
+"""
+
+import numpy as np
+from numba import njit
+from numpy.typing import NDArray
+
+
+@njit(cache=True)
+def bilinear_gather(img, y0, x0, fy, fx, H, W):
+    """Bilinear sample of ``img`` at points with precomputed floor/frac indices.
+
+    Reproduces ``scipy.ndimage.map_coordinates(img, (y, x), order=1)`` with the
+    default ``mode="constant"`` (out-of-range corners contribute ``cval=0``):
+    ``y0=floor(y)``, ``fy=y-y0`` (likewise x). Matches map_coordinates to ~3e-16
+    (FP add-order only; ``tasks/poc_granularity_gather.py``) — far inside this
+    backend's ``rtol=1e-6`` contract — but ~9x faster per call. The floor/frac are
+    computed once on the host and reused across the 16 granularity steps (same
+    sample points, different ``img``), replacing 16 ``map_coordinates`` calls that
+    re-derive them every time (~67% of the default-config runtime).
+    """
+    out = np.empty(y0.shape[0])
+    for t in range(y0.shape[0]):
+        iy = y0[t]
+        ix = x0[t]
+        gy = fy[t]
+        gx = fx[t]
+        v00 = img[iy, ix] if (0 <= iy < H and 0 <= ix < W) else 0.0
+        v01 = img[iy, ix + 1] if (0 <= iy < H and 0 <= ix + 1 < W) else 0.0
+        v10 = img[iy + 1, ix] if (0 <= iy + 1 < H and 0 <= ix < W) else 0.0
+        v11 = img[iy + 1, ix + 1] if (0 <= iy + 1 < H and 0 <= ix + 1 < W) else 0.0
+        out[t] = (v00 * (1.0 - gx) + v01 * gx) * (1.0 - gy) + (
+            v10 * (1.0 - gx) + v11 * gx
+        ) * gy
+    return out
+
+
+def _disk_halfwidths(radius: int) -> NDArray[np.int64]:
+    """Per-row half-widths of ``skimage.morphology.disk(radius)``.
+
+    ``hx[dy] = floor(sqrt(radius² - dy²))`` for ``dy`` in ``0..radius`` — the disk
+    row at vertical offset ``±dy`` spans columns ``[-hx[dy], +hx[dy]]``.
+    """
+    dy = np.arange(radius + 1, dtype=np.float64)
+    return np.floor(np.sqrt(radius * radius - dy * dy)).astype(np.int64)
+
+
+@njit(cache=True)
+def _disk_reduce(P, radius, hx, H, W, is_max):
+    """Row-decomposed VHG disk min/max over an ``edge``-padded image ``P``.
+
+    ``P`` has shape ``(H+2r, W+2r)``. For each of the ``2r+1`` disk rows, a 1-D
+    sliding min/max of radius ``hx[|dy|]`` is taken along the padded row via VHG
+    (prefix ``g`` / suffix ``h`` over blocks of size ``2w+1``), then reduced across
+    rows. ``is_max`` selects dilation (max) vs erosion (min).
+    """
+    L = W + 2 * radius
+    out = np.full((H, W), -np.inf if is_max else np.inf)
+
+    g = np.empty(L, np.float64)
+    h = np.empty(L, np.float64)
+    for dy in range(-radius, radius + 1):
+        w = hx[dy] if dy >= 0 else hx[-dy]
+        k = 2 * w + 1
+        h_block_start = (
+            L - 1
+        ) % k  # x % k of the last element (one modulo per row-band)
+        for i in range(H):
+            pr = i + radius + dy
+            # forward prefix-min/max within blocks of size k (counter tracks x % k)
+            c = 0
+            for x in range(L):
+                if c == 0:
+                    g[x] = P[pr, x]
+                else:
+                    a = g[x - 1]
+                    b = P[pr, x]
+                    if is_max:
+                        g[x] = a if a > b else b
+                    else:
+                        g[x] = a if a < b else b
+                c += 1
+                if c == k:
+                    c = 0
+            # backward suffix-min/max within blocks of size k (counter tracks x % k)
+            c = h_block_start
+            for x in range(L - 1, -1, -1):
+                if x == L - 1 or c == k - 1:
+                    h[x] = P[pr, x]
+                else:
+                    a = h[x + 1]
+                    b = P[pr, x]
+                    if is_max:
+                        h[x] = a if a > b else b
+                    else:
+                        h[x] = a if a < b else b
+                c -= 1
+                if c < 0:
+                    c = k - 1
+            # combine windowed result into the output (center c = radius + j)
+            for j in range(W):
+                c = radius + j
+                a = h[c - w]
+                b = g[c + w]
+                if is_max:
+                    v = a if a > b else b
+                    if v > out[i, j]:
+                        out[i, j] = v
+                else:
+                    v = a if a < b else b
+                    if v < out[i, j]:
+                        out[i, j] = v
+    return out
+
+
+def disk_erosion_2d(img: NDArray, radius: int) -> NDArray[np.float64]:
+    """Greyscale disk erosion, bit-exact with ``skimage…erosion(img, disk(radius))``."""
+    a = np.ascontiguousarray(img, dtype=np.float64)
+    if radius <= 0:
+        return a.copy()
+    P = np.pad(a, radius, mode="edge")
+    return _disk_reduce(
+        P, radius, _disk_halfwidths(radius), a.shape[0], a.shape[1], False
+    )
+
+
+def disk_dilation_2d(img: NDArray, radius: int) -> NDArray[np.float64]:
+    """Greyscale disk dilation, bit-exact with ``skimage…dilation(img, disk(radius))``."""
+    a = np.ascontiguousarray(img, dtype=np.float64)
+    if radius <= 0:
+        return a.copy()
+    P = np.pad(a, radius, mode="edge")
+    return _disk_reduce(
+        P, radius, _disk_halfwidths(radius), a.shape[0], a.shape[1], True
+    )
+
+
+@njit(cache=True)
+def _plus_reduce(P, H, W, is_max):
+    """5-tap (disk(1)/cross) min/max over a 1-pixel ``edge``-padded image ``P``."""
+    out = np.empty((H, W), np.float64)
+    for i in range(H):
+        for j in range(W):
+            v = P[i + 1, j + 1]
+            u = P[i, j + 1]
+            d = P[i + 2, j + 1]
+            le = P[i + 1, j]
+            ri = P[i + 1, j + 2]
+            if is_max:
+                if u > v:
+                    v = u
+                if d > v:
+                    v = d
+                if le > v:
+                    v = le
+                if ri > v:
+                    v = ri
+            else:
+                if u < v:
+                    v = u
+                if d < v:
+                    v = d
+                if le < v:
+                    v = le
+                if ri < v:
+                    v = ri
+            out[i, j] = v
+    return out
+
+
+def erosion_4conn_2d(img: NDArray) -> NDArray[np.float64]:
+    """disk(1) erosion, bit-exact with ``skimage…erosion(img, disk(1))``."""
+    a = np.ascontiguousarray(img, dtype=np.float64)
+    P = np.pad(a, 1, mode="edge")
+    return _plus_reduce(P, a.shape[0], a.shape[1], False)
+
+
+def dilation_4conn_2d(img: NDArray) -> NDArray[np.float64]:
+    """disk(1) dilation, bit-exact with ``skimage…dilation(img, disk(1))``."""
+    a = np.ascontiguousarray(img, dtype=np.float64)
+    P = np.pad(a, 1, mode="edge")
+    return _plus_reduce(P, a.shape[0], a.shape[1], True)
+
+
+@njit(cache=True)
+def _geodesic_raster_fwd(out, mask):
+    """One forward (TL→BR) geodesic-dilation raster of ``out`` under ``mask``."""
+    H, W = out.shape
+    for i in range(H):
+        for j in range(W):
+            v = out[i, j]
+            if i > 0 and out[i - 1, j] > v:
+                v = out[i - 1, j]
+            if j > 0 and out[i, j - 1] > v:
+                v = out[i, j - 1]
+            m = mask[i, j]
+            if v > m:
+                v = m
+            out[i, j] = v
+
+
+@njit(cache=True)
+def _geodesic_raster_bwd(out, mask):
+    """One backward (BR→TL) geodesic-dilation raster of ``out`` under ``mask``."""
+    H, W = out.shape
+    for i in range(H - 1, -1, -1):
+        for j in range(W - 1, -1, -1):
+            v = out[i, j]
+            if i < H - 1 and out[i + 1, j] > v:
+                v = out[i + 1, j]
+            if j < W - 1 and out[i, j + 1] > v:
+                v = out[i, j + 1]
+            m = mask[i, j]
+            if v > m:
+                v = m
+            out[i, j] = v
+
+
+@njit(cache=True)
+def reconstruction_by_dilation_2d(seed, mask):
+    """Morphological reconstruction by dilation (4-connectivity), seed under mask.
+
+    Vincent (1993) / Robinson (2004) hybrid: three forward+backward raster
+    geodesic-dilation passes under ``mask``, the last backward pass seeding a FIFO
+    of pixels that can still propagate, then FIFO propagation until it drains. This
+    computes the exact reconstruction in ``O(N)`` (independent of propagation
+    distance), unlike raster-until-convergence which costs ``O(passes·N)`` and
+    degrades to hundreds of passes on full-resolution images. Bit-identical to
+    ``skimage.morphology.reconstruction(seed, mask, footprint=disk(1))`` — min/max
+    selections only. ``seed`` is clamped to ``min(seed, mask)``.
+
+    Three raster pairs (vs a single pair) cut the FIFO seed count ~50-95%, trading a
+    little cheap cache-sequential scanning for far fewer random-access FIFO updates
+    (~1.5x faster reconstruction, bit-identical — the result is independent of how
+    many raster passes precede the FIFO). FIFO entries are packed as ``(i << 16) | j``
+    int32 (decoded by shift/mask, not int-div/mod; valid for sides ≤ 65535). The FIFO
+    is a ring buffer with a per-pixel ``in-queue`` flag, so at most ``N`` entries are
+    live at once and a buffer of ``N + 1`` cannot overflow.
+    """
+    H, W = seed.shape
+    out = np.empty((H, W), np.float64)
+    for i in range(H):
+        for j in range(W):
+            s = seed[i, j]
+            m = mask[i, j]
+            out[i, j] = s if s < m else m
+
+    # 3 raster pairs before FIFO seeding; the first two pairs do not seed
+    _geodesic_raster_fwd(out, mask)
+    _geodesic_raster_bwd(out, mask)
+    _geodesic_raster_fwd(out, mask)
+    _geodesic_raster_bwd(out, mask)
+    _geodesic_raster_fwd(out, mask)
+
+    cap = H * W + 1
+    queue = np.empty(cap, np.int32)  # packed (i << 16) | j
+    inq = np.zeros((H, W), np.uint8)
+    mask16 = np.int32(0xFFFF)
+    head = 0
+    tail = 0
+
+    # final backward raster: propagate from down/right; seed the FIFO with pixels
+    # whose down/right neighbour can still grow from them
+    for i in range(H - 1, -1, -1):
+        for j in range(W - 1, -1, -1):
+            v = out[i, j]
+            if i < H - 1 and out[i + 1, j] > v:
+                v = out[i + 1, j]
+            if j < W - 1 and out[i, j + 1] > v:
+                v = out[i, j + 1]
+            m = mask[i, j]
+            if v > m:
+                v = m
+            out[i, j] = v
+            seed_p = (
+                i < H - 1 and out[i + 1, j] < v and out[i + 1, j] < mask[i + 1, j]
+            ) or (j < W - 1 and out[i, j + 1] < v and out[i, j + 1] < mask[i, j + 1])
+            if seed_p:
+                queue[tail] = np.int32((i << 16) | j)
+                inq[i, j] = 1
+                tail += 1
+                if tail == cap:
+                    tail = 0
+
+    # FIFO propagation to all 4 neighbours until drained
+    while head != tail:
+        code = queue[head]
+        head += 1
+        if head == cap:
+            head = 0
+        i = code >> np.int32(16)
+        j = code & mask16
+        inq[i, j] = 0
+        v = out[i, j]
+        # up, down, left, right
+        if i > 0:
+            mq = mask[i - 1, j]
+            nv = v if v < mq else mq
+            if nv > out[i - 1, j]:
+                out[i - 1, j] = nv
+                if inq[i - 1, j] == 0:
+                    inq[i - 1, j] = 1
+                    queue[tail] = np.int32(((i - 1) << 16) | j)
+                    tail += 1
+                    if tail == cap:
+                        tail = 0
+        if i < H - 1:
+            mq = mask[i + 1, j]
+            nv = v if v < mq else mq
+            if nv > out[i + 1, j]:
+                out[i + 1, j] = nv
+                if inq[i + 1, j] == 0:
+                    inq[i + 1, j] = 1
+                    queue[tail] = np.int32(((i + 1) << 16) | j)
+                    tail += 1
+                    if tail == cap:
+                        tail = 0
+        if j > 0:
+            mq = mask[i, j - 1]
+            nv = v if v < mq else mq
+            if nv > out[i, j - 1]:
+                out[i, j - 1] = nv
+                if inq[i, j - 1] == 0:
+                    inq[i, j - 1] = 1
+                    queue[tail] = np.int32((i << 16) | (j - 1))
+                    tail += 1
+                    if tail == cap:
+                        tail = 0
+        if j < W - 1:
+            mq = mask[i, j + 1]
+            nv = v if v < mq else mq
+            if nv > out[i, j + 1]:
+                out[i, j + 1] = nv
+                if inq[i, j + 1] == 0:
+                    inq[i, j + 1] = 1
+                    queue[tail] = np.int32((i << 16) | (j + 1))
+                    tail += 1
+                    if tail == cap:
+                        tail = 0
+    return out
diff --git a/src/cp_measure/core/numba/measuregranularity.py b/src/cp_measure/core/numba/measuregranularity.py
new file mode 100644
index 0000000..89dd4d4
--- /dev/null
+++ b/src/cp_measure/core/numba/measuregranularity.py
@@ -0,0 +1,196 @@
+"""Numba-accelerated granularity backend (2D).
+
+Mirrors :func:`cp_measure.core.measuregranularity.get_granularity` but swaps the
+skimage morphology (background opening, per-iteration disk(1) erosion, geodesic
+reconstruction) for the bit-exact numba kernels in :mod:`._granularity`, and the
+16 ``scipy.ndimage.mean`` calls for a precomputed ``bincount`` per-object mean
+sampled by a sparse point-query of the reconstructed image (no full-resolution
+upsample). The 16 per-step resamples of the reconstructed image back to the object
+sample points — the default-config hot spot (~67% of runtime) — use the numba
+:func:`._granularity.bilinear_gather` (floor/frac precomputed once and reused
+across steps) instead of ``scipy.ndimage.map_coordinates``; it matches order-1
+``map_coordinates`` to ~3e-16 (≪ this backend's ``rtol=1e-6`` contract). The three
+one-off grid resamples (subsample / background) stay on ``map_coordinates``.
+
+Batch-shaped via the canonical ``(B, Z, Y, X)`` form (see
+:func:`cp_measure.primitives.shapes.to_bzyx`): a single image is ``B == 1``. Each
+``(Z, Y, X)`` element with ``Z == 1`` runs the numba 2D path; ``Z > 1`` (a true
+3D volume) falls back to the numpy baseline.
+"""
+
+import numpy
+import scipy.ndimage
+from numpy.typing import NDArray
+
+from cp_measure.core.measuregranularity import get_granularity as _get_granularity_numpy
+from cp_measure.core.numba._granularity import (
+    bilinear_gather,
+    disk_dilation_2d,
+    disk_erosion_2d,
+    erosion_4conn_2d,
+    reconstruction_by_dilation_2d,
+)
+from cp_measure.primitives.shapes import to_bzyx
+
+
+def get_granularity(
+    masks,
+    pixels,
+    subsample_size: float = 0.25,
+    image_sample_size: float = 0.25,
+    element_size: int = 10,
+    granular_spectrum_length: int = 16,
+):
+    """Granularity spectrum per object; accepts a single image/volume or a batch.
+
+    Returns a single ``{Granularity_k: array}`` 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). Output arrays are indexed densely by label ``1..max(mask)``.
+    """
+    masks_zyx, pixels_zyx, unwrap = to_bzyx(masks, pixels)
+    results = [
+        _granularity_one(
+            m,
+            p,
+            subsample_size,
+            image_sample_size,
+            element_size,
+            granular_spectrum_length,
+        )
+        for m, p in zip(masks_zyx, pixels_zyx)
+    ]
+    return unwrap(results)
+
+
+def _granularity_one(
+    mask_zyx, pixels_zyx, subsample_size, image_sample_size, element_size, ng
+):
+    """Dispatch one ``(Z, Y, X)`` element: numba 2D path, else numpy baseline (3D)."""
+    if mask_zyx.shape[0] == 1:  # Z == 1 -> 2D image
+        return _granularity_2d(
+            mask_zyx[0],
+            pixels_zyx[0],
+            subsample_size,
+            image_sample_size,
+            element_size,
+            ng,
+        )
+    return _get_granularity_numpy(
+        mask_zyx, pixels_zyx, subsample_size, image_sample_size, element_size, ng
+    )
+
+
+def _granularity_2d(
+    orig_mask: NDArray[numpy.integer],
+    orig_pixels: NDArray[numpy.floating],
+    subsample_size: float,
+    image_sample_size: float,
+    element_size: int,
+    ng: int,
+) -> dict[str, NDArray[numpy.floating]]:
+    orig_shape = numpy.array(orig_pixels.shape)
+    new_shape = orig_shape.copy()
+
+    # 1. Subsample image (scipy map_coordinates, as baseline). The baseline also
+    # resamples the mask here but never reads it (the spectrum loop uses orig_mask),
+    # so we skip that dead order-0 resample.
+    if subsample_size < 1:
+        new_shape = (orig_shape * subsample_size).astype(int)
+        i, j = (
+            numpy.mgrid[0 : new_shape[0], 0 : new_shape[1]].astype(float)
+            / subsample_size
+        )
+        pixels = scipy.ndimage.map_coordinates(orig_pixels, (i, j), order=1)
+    else:
+        pixels = orig_pixels.astype(numpy.float64, copy=True)
+
+    # 2. Background subtraction via greyscale opening (numba disk erosion+dilation)
+    if image_sample_size < 1:
+        back_shape = new_shape * image_sample_size
+        i, j = (
+            numpy.mgrid[0 : back_shape[0], 0 : back_shape[1]].astype(float)
+            / image_sample_size
+        )
+        back_pixels = scipy.ndimage.map_coordinates(pixels, (i, j), order=1)
+    else:
+        back_pixels = pixels
+        back_shape = new_shape
+    radius = element_size
+    # back_mask is all-ones (full-frame), so the baseline's masked re-zeroing is a no-op
+    back_pixels = disk_erosion_2d(back_pixels, radius)
+    back_pixels = disk_dilation_2d(back_pixels, radius)
+    if image_sample_size < 1:
+        i, j = numpy.mgrid[0 : new_shape[0], 0 : new_shape[1]].astype(float)
+        i *= float(back_shape[0] - 1) / float(new_shape[0] - 1)
+        j *= float(back_shape[1] - 1) / float(new_shape[1] - 1)
+        back_pixels = scipy.ndimage.map_coordinates(back_pixels, (i, j), order=1)
+    pixels = pixels - back_pixels
+    pixels[pixels < 0] = 0
+
+    # 3. Per-object readback structures (dense labels 1..max, like baseline)
+    flat_mask = orig_mask.ravel()
+    in_obj = flat_mask > 0
+    max_label = int(orig_mask.max()) if orig_mask.size else 0
+    results: dict[str, NDArray[numpy.floating]] = {}
+    if max_label == 0:  # no objects -> baseline returns empty arrays
+        empty = numpy.zeros((0,))
+        for granularity_id in range(1, ng + 1):
+            results[f"Granularity_{granularity_id}"] = empty
+        return results
+
+    flat_pos = numpy.flatnonzero(in_obj)
+    labels_in = flat_mask[flat_pos]
+    counts = numpy.bincount(labels_in, minlength=max_label + 1)[1:].astype(
+        numpy.float64
+    )
+
+    needs_resize = not numpy.array_equal(new_shape, orig_shape)
+    if needs_resize:
+        yy, xx = numpy.unravel_index(flat_pos, tuple(orig_shape))
+        sy = yy * (float(new_shape[0] - 1) / float(orig_shape[0] - 1))
+        sx = xx * (float(new_shape[1] - 1) / float(orig_shape[1] - 1))
+        # The sample points are fixed across the 16 spectrum steps (only `rec`
+        # changes), so the bilinear floor/frac are computed ONCE and reused, instead
+        # of 16 map_coordinates calls re-deriving them (the default-config hot spot).
+        gy0 = numpy.floor(sy).astype(numpy.int64)
+        gx0 = numpy.floor(sx).astype(numpy.int64)
+        gfy = sy - gy0
+        gfx = sx - gx0
+
+    def _label_mean(values_at_in_obj):
+        sums = numpy.bincount(
+            labels_in, weights=values_at_in_obj, minlength=max_label + 1
+        )[1:]
+        with numpy.errstate(invalid="ignore", divide="ignore"):
+            return (
+                sums / counts
+            )  # 0/0 -> nan for absent labels, matching scipy.ndimage.mean
+
+    orig_valid = orig_pixels.ravel()[flat_pos].astype(numpy.float64)
+    current_mean = _label_mean(orig_valid)
+    start_mean = numpy.maximum(current_mean, numpy.finfo(float).eps)
+
+    # 4. Granular spectrum loop (numba disk(1) erosion + reconstruction, cascaded mask)
+    ero = pixels
+    recon_mask = pixels
+    for granularity_id in range(1, ng + 1):
+        ero = erosion_4conn_2d(ero)
+        # reconstruction of an all-zero seed returns zeros cheaply (no FIFO seeds),
+        # so no all-zero early-out is needed (it would only add a full ero.max() scan)
+        rec = reconstruction_by_dilation_2d(ero, recon_mask)
+        recon_mask = rec  # cascade: rec_g <= rec_{g-1} <= pixels (exact)
+
+        if needs_resize:
+            rec_valid = bilinear_gather(
+                rec, gy0, gx0, gfy, gfx, rec.shape[0], rec.shape[1]
+            )
+        else:
+            rec_valid = rec.ravel()[flat_pos]
+
+        new_mean = _label_mean(rec_valid)
+        results[f"Granularity_{granularity_id}"] = (
+            (current_mean - new_mean) * 100 / start_mean
+        )
+        current_mean = new_mean
+
+    return results
diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py
new file mode 100644
index 0000000..bd058c5
--- /dev/null
+++ b/src/cp_measure/core/numba/measureobjectintensity.py
@@ -0,0 +1,199 @@
+"""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._segment_numba import (
+    flatten_numba,
+    inner_boundary,
+    segment_moments,
+    segment_quantiles,
+    segment_resid_sumsq,
+    segment_stats,
+)
+
+
+def get_intensity(
+    masks: NDArray[numpy.integer],
+    pixels: NDArray[numpy.floating],
+    edge_measurements: bool = True,
+) -> dict[str, NDArray[numpy.floating]]:
+    """masks is a labeled array where 0 are background."""
+    orig_ndim = pixels.ndim
+
+    masked_image = pixels
+    if pixels.ndim == 2:
+        masked_image = pixels.reshape(1, *pixels.shape)
+        if masks.ndim == 2:
+            masks = masks.reshape(1, *masks.shape)
+    elif pixels.ndim == 3 and masks.ndim == 2:  # 3D image, 2D mask
+        masks = masks.reshape(1, *masks.shape)
+
+    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(masked_image),
+        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 = masked_image[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
+    }
diff --git a/src/cp_measure/primitives/__init__.py b/src/cp_measure/primitives/__init__.py
new file mode 100644
index 0000000..3619525
--- /dev/null
+++ b/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.
+"""
diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py
new file mode 100644
index 0000000..ffb9148
--- /dev/null
+++ b/src/cp_measure/primitives/_segment_numba.py
@@ -0,0 +1,228 @@
+"""Numba segment kernels (single-threaded, cached).
+
+These are the numba implementation of the segment-reduce / segment-quantile
+primitives. They loop over the flat ``(values, seg0, coords)`` arrays produced
+by :mod:`cp_measure.primitives.segment` — no image shape, no batch axis — so one
+kernel set covers 2D, 3D, and (future) batched inputs unchanged.
+
+All kernels are ``@njit(cache=True)`` and serial: no ``prange``/``nogil``.
+Parallelism is the job of the (future) batch layer over images, not the kernel.
+"""
+
+import numpy as np
+from numba import njit
+
+
+@njit(cache=True)
+def flatten_numba(masks, pixels, lut):
+    """Flatten a labeled (Z, Y, X) image to ``(values, seg0, xc, yc, zc)``.
+
+    Two grid scans (count, then fill) replace the numpy ``(masks>0)&isfinite``
+    mask + ``numpy.nonzero`` + fancy-index gathers; coordinates are the loop
+    indices. Background and non-finite pixels are dropped, in C (raster) order.
+    ``masks`` and ``pixels`` must be C-contiguous; ``pixels`` may be any float
+    dtype (kept values are upcast into the float64 ``values`` output).
+    """
+    Z, Y, X = masks.shape
+    M = 0
+    for z in range(Z):
+        for y in range(Y):
+            for x in range(X):
+                if masks[z, y, x] > 0 and np.isfinite(pixels[z, y, x]):
+                    M += 1
+    values = np.empty(M, np.float64)
+    seg0 = np.empty(M, np.int64)
+    xc = np.empty(M, np.float64)
+    yc = np.empty(M, np.float64)
+    zc = np.empty(M, np.float64)
+    i = 0
+    for z in range(Z):
+        for y in range(Y):
+            for x in range(X):
+                L = masks[z, y, x]
+                if L <= 0:
+                    continue
+                v = pixels[z, y, x]
+                if not np.isfinite(v):
+                    continue
+                values[i] = v
+                seg0[i] = lut[L]
+                xc[i] = x
+                yc[i] = y
+                zc[i] = z
+                i += 1
+    return values, seg0, xc, yc, zc
+
+
+@njit(cache=True)
+def segment_moments(values, seg0, xc, yc, zc, n):
+    """One pass over the flat pixels accumulating, per segment:
+
+    count, sum, min, max, the max-pixel coordinates, and the six centroid
+    cross-sums (Sx, Sy, Sz, Sx*I, Sy*I, Sz*I).
+
+    The max position uses ``>=`` so the LAST pixel (in the flat raster order the
+    host builds the arrays in) attaining the maximum wins. On real, continuous
+    data the max is unique, so this is bit-identical to the numpy backend's
+    ``scipy.ndimage.maximum_position``; on exact-value ties it is a deterministic
+    rule rather than scipy's quicksort-dependent (version-unstable) pick.
+    """
+    count = np.zeros(n, np.int64)
+    sumI = np.zeros(n, np.float64)
+    minI = np.full(n, np.inf)
+    maxI = np.full(n, -np.inf)
+    mx = np.zeros(n, np.float64)
+    my = np.zeros(n, np.float64)
+    mz = np.zeros(n, np.float64)
+    sx = np.zeros(n, np.float64)
+    sy = np.zeros(n, np.float64)
+    sz = np.zeros(n, np.float64)
+    sxI = np.zeros(n, np.float64)
+    syI = np.zeros(n, np.float64)
+    szI = np.zeros(n, np.float64)
+    M = values.shape[0]
+    for i in range(M):
+        k = seg0[i]
+        v = values[i]
+        x = xc[i]
+        y = yc[i]
+        z = zc[i]
+        count[k] += 1
+        sumI[k] += v
+        if v < minI[k]:
+            minI[k] = v
+        if v >= maxI[k]:  # >= -> keep LAST max in raster order
+            maxI[k] = v
+            mx[k] = x
+            my[k] = y
+            mz[k] = z
+        sx[k] += x
+        sy[k] += y
+        sz[k] += z
+        sxI[k] += x * v
+        syI[k] += y * v
+        szI[k] += z * v
+    return count, sumI, minI, maxI, mx, my, mz, sx, sy, sz, sxI, syI, szI
+
+
+@njit(cache=True)
+def segment_stats(values, seg0, n):
+    """One pass accumulating per-segment count, sum, min, max only.
+
+    The lightweight reduce for callers (e.g. edge measurements) that need
+    neither the max position nor the centroid cross-sums.
+    """
+    count = np.zeros(n, np.int64)
+    sumI = np.zeros(n, np.float64)
+    minI = np.full(n, np.inf)
+    maxI = np.full(n, -np.inf)
+    M = values.shape[0]
+    for i in range(M):
+        k = seg0[i]
+        v = values[i]
+        count[k] += 1
+        sumI[k] += v
+        if v < minI[k]:
+            minI[k] = v
+        if v > maxI[k]:
+            maxI[k] = v
+    return count, sumI, minI, maxI
+
+
+@njit(cache=True)
+def inner_boundary(masks2d):
+    """Labeled inner boundary of a 2D label plane (0 = not a boundary pixel).
+
+    A foreground pixel is a boundary pixel iff any in-bounds 4-neighbour has a
+    different label; out-of-bounds neighbours are ignored. Bit-identical to
+    ``skimage.segmentation.find_boundaries(mode="inner")`` on a single 2D plane,
+    ~12-27x faster (one pass, no morphology). ``masks2d`` must be C-contiguous.
+    """
+    H, W = masks2d.shape
+    out = np.zeros((H, W), masks2d.dtype)
+    for r in range(H):
+        for c in range(W):
+            L = masks2d[r, c]
+            if L <= 0:
+                continue
+            if (
+                (r > 0 and masks2d[r - 1, c] != L)
+                or (r < H - 1 and masks2d[r + 1, c] != L)
+                or (c > 0 and masks2d[r, c - 1] != L)
+                or (c < W - 1 and masks2d[r, c + 1] != L)
+            ):
+                out[r, c] = L
+    return out
+
+
+@njit(cache=True)
+def segment_resid_sumsq(values, seg0, n, mean):
+    """Second pass: per-segment sum of squared residuals (for population std)."""
+    ss = np.zeros(n, np.float64)
+    M = values.shape[0]
+    for i in range(M):
+        k = seg0[i]
+        d = values[i] - mean[k]
+        ss[k] += d * d
+    return ss
+
+
+@njit(cache=True)
+def _interp(sorted_seg, n, frac):
+    """Linear interpolation at position ``n * frac`` within a sorted segment.
+
+    Matches the reference's ``cumsum(areas)``-offset interpolation: the local
+    position is ``count * frac``; clamp the upper neighbour at the last element.
+    """
+    pos = n * frac
+    lo = int(pos)
+    if lo > n - 1:
+        lo = n - 1
+    hi = lo + 1
+    if hi > n - 1:
+        hi = n - 1
+    f = pos - lo
+    return sorted_seg[lo] * (1.0 - f) + sorted_seg[hi] * f
+
+
+@njit(cache=True)
+def segment_quantiles(values, seg0, counts, n, mad_frac):
+    """Per-segment quartiles/median + MAD via scatter-into-segments + sort.
+
+    Builds CSR offsets from ``counts``, scatters values into one flat buffer,
+    sorts each segment slice in place, and interpolates. MAD reuses the sorted
+    absolute deviations from the median at ``mad_frac`` (= 1 / original ndim).
+    """
+    starts = np.zeros(n, np.int64)
+    acc = 0
+    for k in range(n):
+        starts[k] = acc
+        acc += counts[k]
+    total = acc
+
+    buf = np.empty(total, np.float64)
+    cursor = starts.copy()
+    M = values.shape[0]
+    for i in range(M):
+        k = seg0[i]
+        buf[cursor[k]] = values[i]
+        cursor[k] += 1
+
+    lq = np.zeros(n, np.float64)
+    med = np.zeros(n, np.float64)
+    uq = np.zeros(n, np.float64)
+    mad = np.zeros(n, np.float64)
+    for k in range(n):
+        cnt = counts[k]
+        if cnt == 0:
+            continue
+        s = starts[k]
+        seg = buf[s : s + cnt]
+        seg.sort()
+        lq[k] = _interp(seg, cnt, 0.25)
+        med[k] = _interp(seg, cnt, 0.5)
+        uq[k] = _interp(seg, cnt, 0.75)
+        ad = np.abs(seg - med[k])
+        ad.sort()
+        mad[k] = _interp(ad, cnt, mad_frac)
+    return lq, med, uq, mad
diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py
new file mode 100644
index 0000000..e5b5894
--- /dev/null
+++ b/src/cp_measure/primitives/segment.py
@@ -0,0 +1,39 @@
+"""Host-side segment helpers (numpy, backend-agnostic).
+
+A labeled image is reduced to a flat *segment* form — ``values[M]`` intensities,
+``seg0[M]`` 0-based segment indices, and ``xc/yc/zc[M]`` per-pixel coordinates —
+which the segment kernels then reduce. The flattening itself lives in the numba
+layer (:func:`cp_measure.primitives._segment_numba.flatten_numba`); this module
+holds only the numpy label→index lookup that feeds it.
+"""
+
+import numpy
+import scipy.ndimage
+from numpy.typing import NDArray
+
+
+def label_to_idx_lut(
+    masks: NDArray[numpy.integer],
+) -> tuple[NDArray[numpy.int64], int]:
+    """Build a ``label -> 0..n-1`` lookup over the sorted positive labels.
+
+    Returns ``(lut, n)`` where ``lut[label]`` is the segment index (and ``-1``
+    for absent labels / background) and ``n`` is the object count. Output arrays
+    are indexed by this segment rank, which equals ``label - 1`` when labels are
+    the contiguous ``1..n`` that cp_measure expects.
+
+    Present labels are enumerated with ``scipy.ndimage.find_objects`` (one pass)
+    rather than ``numpy.unique`` (a full-image sort): same ascending label set,
+    identical LUT, ~3-5x faster on large/3D masks.
+    """
+    if not numpy.issubdtype(masks.dtype, numpy.integer):
+        masks = masks.astype(numpy.intp, copy=False)
+    bboxes = scipy.ndimage.find_objects(masks)
+    labels = numpy.array(
+        [i + 1 for i, sl in enumerate(bboxes) if sl is not None], dtype=numpy.int64
+    )
+    n = int(labels.size)
+    max_label = int(labels[-1]) if n else 0
+    lut = numpy.full(max_label + 1, -1, dtype=numpy.int64)
+    lut[labels] = numpy.arange(n, dtype=numpy.int64)
+    return lut, n
diff --git a/src/cp_measure/primitives/shapes.py b/src/cp_measure/primitives/shapes.py
new file mode 100644
index 0000000..0c94059
--- /dev/null
+++ b/src/cp_measure/primitives/shapes.py
@@ -0,0 +1,83 @@
+"""Canonical ``(B, Z, Y, X)`` input normalisation (numpy, backend-agnostic).
+
+Every optimised feature works on a batch of volumes. Rather than a single dense
+``(B, Z, Y, X)`` array (which cannot hold ragged, differently-sized images), the
+batch is represented as a **list of B ``(Z, Y, X)`` arrays** — equal-size and
+ragged batches share one code path, and a single image is just ``B == 1``.
+
+Normalisation rules (the only public entry, :func:`to_bzyx`):
+
+================================ ============================ ========
+input                            yields                       batch?
+================================ ============================ ========
+2D ``(H, W)`` ndarray            ``[(1, H, W)]``              no
+3D ``(Z, Y, X)`` ndarray         ``[(Z, Y, X)]`` (one volume) no
+4D ``(B, Z, Y, X)`` ndarray      ``[(Z, Y, X)] * B``          yes
+list/tuple of 2D/3D arrays       one ``(Z, Y, X)`` per item   yes
+================================ ============================ ========
+
+A 3D ndarray is therefore ALWAYS one volume, never a batch — this preserves the
+existing single-volume semantics. To pass a batch of 2D images as an array, use
+``(B, 1, H, W)``; for ragged sizes, pass a list. ``unwrap`` then re-shapes the
+per-image results back to a single dict (single input) or the list (batch).
+
+This is a pure batch normaliser: ``masks`` and ``pixels`` must share the same
+batch/ndim structure. It does NOT broadcast a lower-dimensional mask over a
+volume (e.g. a 2D mask applied to a 3D stack) — per-element ndim handling is the
+caller's job (each backend dispatches its own 2D vs 3D path).
+"""
+
+import numpy
+from numpy.typing import NDArray
+
+
+def _to_zyx(arr: NDArray) -> NDArray:
+    """Promote a single 2D/3D image to ``(Z, Y, X)`` (2D gets a unit Z axis)."""
+    a = numpy.asarray(arr)
+    if a.ndim == 2:
+        return a[numpy.newaxis]
+    if a.ndim == 3:
+        return a
+    raise ValueError(f"expected a 2D or 3D image, got ndim={a.ndim}")
+
+
+def to_bzyx(masks, pixels):
+    """Normalise ``(masks, pixels)`` to the canonical batch-of-volumes form.
+
+    Returns ``(masks_zyx, pixels_zyx, unwrap)`` where ``masks_zyx`` and
+    ``pixels_zyx`` are length-``B`` lists of ``(Z, Y, X)`` arrays (one per image),
+    and ``unwrap(results)`` maps a length-``B`` list of per-image results back to
+    a single result (non-batch input) or the list itself (batch input).
+    """
+    masks_is_seq = isinstance(masks, (list, tuple))
+    pixels_is_seq = isinstance(pixels, (list, tuple))
+    if masks_is_seq != pixels_is_seq:
+        raise ValueError("masks and pixels must both be sequences, or both arrays")
+
+    if masks_is_seq:
+        masks_zyx = [_to_zyx(m) for m in masks]
+        pixels_zyx = [_to_zyx(p) for p in pixels]
+        is_batch = True
+    else:
+        m = numpy.asarray(masks)
+        p = numpy.asarray(pixels)
+        if (m.ndim == 4) != (p.ndim == 4):
+            raise ValueError("masks and pixels must both be 4D for a stacked batch")
+        if m.ndim == 4:
+            masks_zyx = list(m)
+            pixels_zyx = list(p)
+            is_batch = True
+        else:
+            masks_zyx = [_to_zyx(m)]
+            pixels_zyx = [_to_zyx(p)]
+            is_batch = False
+
+    if len(masks_zyx) != len(pixels_zyx):
+        raise ValueError(
+            f"batch size mismatch: {len(masks_zyx)} masks vs {len(pixels_zyx)} images"
+        )
+
+    def unwrap(results):
+        return results if is_batch else results[0]
+
+    return masks_zyx, pixels_zyx, unwrap
diff --git a/test/test_backend_correctness.py b/test/test_backend_correctness.py
new file mode 100644
index 0000000..52ab275
--- /dev/null
+++ b/test/test_backend_correctness.py
@@ -0,0 +1,111 @@
+"""Backend correctness harness: numba `intensity` must match the numpy backend.
+
+Compares ``cp_measure.core.numba.get_intensity`` against the reference
+``cp_measure.core.get_intensity`` key-by-key on continuous random pixels (no
+exact ties, so ``Location_MaxIntensity_*`` is unambiguous), for 2D and 3D input
+with edge measurements on and off. Also checks the dispatch wiring: under
+``set_accelerator("numba")`` the core registry composes the numba intensity with
+the numpy implementations of every other feature, and an absent numba backend
+raises rather than silently falling back.
+"""
+
+import numpy as np
+import pytest
+from conftest import (
+    DEPTH_3D,
+    SIZE_2D,
+    SIZE_3D,
+    _stamp_objects_2d,
+    _stamp_objects_3d,
+    get_rng,
+)
+
+import cp_measure
+import cp_measure._detect
+import cp_measure.bulk
+from cp_measure._detect import HAS_NUMBA
+from cp_measure.core.measureobjectintensity import get_intensity as intensity_numpy
+
+requires_numba = pytest.mark.skipif(not HAS_NUMBA, reason="numba not installed")
+
+
+def _mask_pixels_2d():
+    mask = np.zeros((SIZE_2D, SIZE_2D), dtype=np.int32)
+    _stamp_objects_2d(mask, n_objects=3)
+    return mask, get_rng().random((SIZE_2D, SIZE_2D))
+
+
+def _mask_pixels_3d():
+    mask = np.zeros((DEPTH_3D, SIZE_3D, SIZE_3D), dtype=np.int32)
+    _stamp_objects_3d(mask, n_objects=2)
+    return mask, get_rng().random((DEPTH_3D, SIZE_3D, SIZE_3D))
+
+
+def _assert_dicts_match(ref, got):
+    assert set(got) == set(ref), set(got).symmetric_difference(ref)
+    for key in ref:
+        np.testing.assert_allclose(
+            got[key], ref[key], rtol=1e-6, atol=1e-8, err_msg=f"feature {key!r}"
+        )
+
+
+@requires_numba
+@pytest.mark.parametrize("edge", [True, False], ids=["edge", "noedge"])
+@pytest.mark.parametrize("dim", ["2d", "3d"])
+def test_numba_intensity_matches_numpy(dim, edge):
+    from cp_measure.core.numba import get_intensity as intensity_numba
+
+    mask, pixels = _mask_pixels_2d() if dim == "2d" else _mask_pixels_3d()
+    ref = intensity_numpy(mask, pixels, edge_measurements=edge)
+    got = intensity_numba(mask, pixels, edge_measurements=edge)
+    _assert_dicts_match(ref, got)
+
+
+@requires_numba
+@pytest.mark.parametrize("sub", [1.0, 0.25], ids=["fullres", "subsampled"])
+def test_numba_granularity_matches_numpy(sub):
+    from cp_measure.core.measuregranularity import get_granularity as gran_numpy
+    from cp_measure.core.numba import get_granularity as gran_numba
+
+    mask, pixels = _mask_pixels_2d()
+    ref = gran_numpy(mask, pixels, subsample_size=sub, image_sample_size=sub)
+    got = gran_numba(mask, pixels, subsample_size=sub, image_sample_size=sub)
+    # bit-exact morphology; per-object mean is a float reduction -> small tol, nan-aware
+    assert set(got) == set(ref)
+    for key in ref:
+        np.testing.assert_allclose(
+            got[key], ref[key], rtol=1e-6, atol=1e-8, equal_nan=True, err_msg=key
+        )
+
+
+@requires_numba
+def test_set_accelerator_numba_composes_with_numpy():
+    cp_measure.set_accelerator("numba")
+    try:
+        core = cp_measure.bulk.get_core_measurements()
+        assert core["intensity"].__module__ == (
+            "cp_measure.core.numba.measureobjectintensity"
+        )
+        assert core["granularity"].__module__ == (
+            "cp_measure.core.numba.measuregranularity"
+        )
+        # Every other feature stays on the numpy backend.
+        assert core["sizeshape"].__module__ == "cp_measure.core.measureobjectsizeshape"
+        assert core["texture"].__module__ == "cp_measure.core.measuretexture"
+    finally:
+        cp_measure.set_accelerator(None)
+
+    restored = cp_measure.bulk.get_core_measurements()
+    assert restored["intensity"].__module__ == "cp_measure.core.measureobjectintensity"
+    assert restored["granularity"].__module__ == "cp_measure.core.measuregranularity"
+
+
+def test_set_accelerator_numba_absent_raises(monkeypatch):
+    """When find_spec reports numba absent, selecting it raises (no silent fallback)."""
+    monkeypatch.setattr(cp_measure._detect, "HAS_NUMBA", False)
+    cp_measure.set_accelerator("numba")
+    try:
+        with pytest.raises(RuntimeError, match="numba is not installed"):
+            cp_measure.bulk.get_core_measurements()
+    finally:
+        cp_measure.set_accelerator(None)
diff --git a/test/test_granularity_backend.py b/test/test_granularity_backend.py
new file mode 100644
index 0000000..7c156f8
--- /dev/null
+++ b/test/test_granularity_backend.py
@@ -0,0 +1,116 @@
+"""Correctness of the numba granularity backend vs the numpy baseline."""
+
+import numpy
+import pytest
+
+from cp_measure.core.measuregranularity import get_granularity as gran_numpy
+from cp_measure.core.numba.measuregranularity import get_granularity as gran_numba
+
+TOL = dict(rtol=1e-7, atol=1e-10, equal_nan=True)
+
+
+def _scene(H=80, W=80, n=4, seed=0):
+    """A textured 2D image + a labeled mask with `n` rectangular objects."""
+    rng = numpy.random.default_rng(seed)
+    pixels = rng.random((H, W)) + 0.3 * numpy.add.outer(
+        numpy.sin(numpy.arange(H) / 4), numpy.cos(numpy.arange(W) / 5)
+    )
+    mask = numpy.zeros((H, W), dtype=numpy.int32)
+    step = H // (n + 1)
+    for k in range(n):
+        y0 = step * (k + 1) - 6
+        x0 = 8 + 15 * k
+        mask[y0 : y0 + 12, x0 : x0 + 12] = k + 1
+    return mask, pixels
+
+
+def _assert_same(a, b):
+    assert a.keys() == b.keys()
+    for k in a:
+        numpy.testing.assert_allclose(a[k], b[k], **TOL, err_msg=f"key {k}")
+
+
+@pytest.mark.parametrize(
+    "sub,img", [(1.0, 1.0), (0.25, 0.25), (0.5, 0.25), (0.25, 1.0)]
+)
+def test_2d_matches_baseline(sub, img):
+    mask, pixels = _scene()
+    _assert_same(
+        gran_numba(mask, pixels, subsample_size=sub, image_sample_size=img),
+        gran_numpy(mask, pixels, subsample_size=sub, image_sample_size=img),
+    )
+
+
+def test_single_object():
+    mask = numpy.zeros((60, 60), numpy.int32)
+    mask[20:40, 20:40] = 1
+    rng = numpy.random.default_rng(1)
+    pixels = rng.random((60, 60))
+    _assert_same(gran_numba(mask, pixels), gran_numpy(mask, pixels))
+
+
+def test_one_pixel_object():
+    mask = numpy.zeros((40, 40), numpy.int32)
+    mask[10, 10] = 1
+    mask[25, 30] = 2
+    pixels = numpy.random.default_rng(2).random((40, 40))
+    _assert_same(gran_numba(mask, pixels), gran_numpy(mask, pixels))
+
+
+def test_label_gaps_dense_output():
+    # labels 1 and 3 present, 2 absent -> dense output length 3, nan at index 1
+    mask = numpy.zeros((50, 50), numpy.int32)
+    mask[5:15, 5:15] = 1
+    mask[30:40, 30:40] = 3
+    pixels = numpy.random.default_rng(3).random((50, 50))
+    out = gran_numba(mask, pixels)
+    base = gran_numpy(mask, pixels)
+    assert out["Granularity_1"].shape == (3,)
+    _assert_same(out, base)
+
+
+def test_empty_mask():
+    mask = numpy.zeros((40, 40), numpy.int32)
+    pixels = numpy.random.default_rng(4).random((40, 40))
+    out = gran_numba(mask, pixels)
+    base = gran_numpy(mask, pixels)
+    assert out["Granularity_1"].shape == (0,)
+    _assert_same(out, base)
+
+
+def test_single_image_returns_dict_not_list():
+    mask, pixels = _scene()
+    assert isinstance(gran_numba(mask, pixels), dict)
+
+
+def test_batch_4d_matches_per_image_baseline():
+    scenes = [_scene(seed=s) for s in range(3)]
+    masks = numpy.stack([m[numpy.newaxis] for m, _ in scenes])  # (B,1,H,W)
+    pixels = numpy.stack([p[numpy.newaxis] for _, p in scenes])
+    out = gran_numba(masks, pixels)
+    assert isinstance(out, list) and len(out) == 3
+    for (m, p), got in zip(scenes, out):
+        _assert_same(got, gran_numpy(m, p))
+
+
+def test_batch_list_ragged():
+    a = _scene(H=80, W=80, seed=0)
+    b = _scene(H=64, W=96, n=3, seed=1)
+    out = gran_numba([a[0], b[0]], [a[1], b[1]])
+    assert isinstance(out, list) and len(out) == 2
+    _assert_same(out[0], gran_numpy(*a))
+    _assert_same(out[1], gran_numpy(*b))
+
+
+def test_3d_volume_falls_back_to_numpy():
+    # Z>1 dispatches to the numpy baseline. Use fullres: the baseline collapses Z
+    # to int(Z*subsample_size) and errors when that hits 0, so a thin volume at the
+    # 0.25 default cannot be subsampled at all (a baseline limitation, not ours).
+    rng = numpy.random.default_rng(5)
+    pixels = rng.random((3, 40, 40))
+    mask = numpy.zeros((3, 40, 40), numpy.int32)
+    mask[:, 10:25, 10:25] = 1
+    _assert_same(
+        gran_numba(mask, pixels, subsample_size=1.0, image_sample_size=1.0),
+        gran_numpy(mask, pixels, subsample_size=1.0, image_sample_size=1.0),
+    )
diff --git a/test/test_granularity_kernels.py b/test/test_granularity_kernels.py
new file mode 100644
index 0000000..1b0150d
--- /dev/null
+++ b/test/test_granularity_kernels.py
@@ -0,0 +1,151 @@
+"""Bit-exact validation of the numba granularity morphology kernels vs skimage."""
+
+import numpy
+import pytest
+import scipy.ndimage
+import skimage.morphology as M
+
+from cp_measure.core.numba._granularity import (
+    _disk_halfwidths,
+    bilinear_gather,
+    disk_dilation_2d,
+    disk_erosion_2d,
+    dilation_4conn_2d,
+    erosion_4conn_2d,
+    reconstruction_by_dilation_2d,
+)
+
+
+@pytest.mark.parametrize("seed", [0, 1, 7])
+def test_bilinear_gather_matches_map_coordinates(seed):
+    """The gather must match scipy order-1 map_coordinates (mode=constant) to ~eps,
+    including points on the array boundary (where a corner is out of range)."""
+    rng = numpy.random.default_rng(seed)
+    img = rng.random((37, 41))
+    H, W = img.shape
+    npts = 500
+    sy = numpy.concatenate([rng.uniform(0, H - 1, npts), [0.0, H - 1, H - 1]])
+    sx = numpy.concatenate([rng.uniform(0, W - 1, npts), [0.0, W - 1, 0.0]])
+    ref = scipy.ndimage.map_coordinates(img, (sy, sx), order=1)
+    y0 = numpy.floor(sy).astype(numpy.int64)
+    x0 = numpy.floor(sx).astype(numpy.int64)
+    got = bilinear_gather(img, y0, x0, sy - y0, sx - x0, H, W)
+    numpy.testing.assert_allclose(got, ref, rtol=0, atol=1e-12)
+
+
+def _images(seed):
+    rng = numpy.random.default_rng(seed)
+    H, W = 23, 19
+    return {
+        "random": rng.random((H, W)),
+        "int_ties": rng.integers(0, 5, size=(H, W)).astype(numpy.float64),
+        "gradient": numpy.add.outer(numpy.arange(H), numpy.arange(W)).astype(float),
+        "border_hot": _border_hot(H, W),
+        "single_spike": _spike(H, W),
+    }
+
+
+def _border_hot(H, W):
+    a = numpy.zeros((H, W))
+    a[0, :] = a[-1, :] = a[:, 0] = a[:, -1] = 7.0
+    return a
+
+
+def _spike(H, W):
+    a = numpy.zeros((H, W))
+    a[H // 2, W // 2] = 5.0
+    a[2, 3] = 3.0
+    return a
+
+
+@pytest.mark.parametrize("radius", [1, 2, 3, 5, 7, 10])
+@pytest.mark.parametrize(
+    "name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]
+)
+def test_disk_erosion_bit_exact(radius, name):
+    img = _images(radius)[name]
+    expected = M.erosion(img, M.disk(radius))
+    numpy.testing.assert_array_equal(disk_erosion_2d(img, radius), expected)
+
+
+@pytest.mark.parametrize("radius", [1, 2, 3, 5, 7, 10])
+@pytest.mark.parametrize(
+    "name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]
+)
+def test_disk_dilation_bit_exact(radius, name):
+    img = _images(radius)[name]
+    expected = M.dilation(img, M.disk(radius))
+    numpy.testing.assert_array_equal(disk_dilation_2d(img, radius), expected)
+
+
+@pytest.mark.parametrize(
+    "name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]
+)
+def test_4conn_erosion_bit_exact(name):
+    img = _images(0)[name]
+    numpy.testing.assert_array_equal(erosion_4conn_2d(img), M.erosion(img, M.disk(1)))
+
+
+@pytest.mark.parametrize(
+    "name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]
+)
+def test_4conn_dilation_bit_exact(name):
+    img = _images(0)[name]
+    numpy.testing.assert_array_equal(dilation_4conn_2d(img), M.dilation(img, M.disk(1)))
+
+
+def test_halfwidths_match_disk():
+    for r in (1, 2, 3, 5, 10):
+        hx = _disk_halfwidths(r)
+        fp = M.disk(r)
+        # per-row count of set pixels == 2*hx+1
+        for dy in range(r + 1):
+            row = fp[r + dy]
+            assert int(row.sum()) == 2 * int(hx[dy]) + 1
+
+
+@pytest.mark.parametrize("name", ["random", "int_ties", "gradient", "single_spike"])
+def test_reconstruction_bit_exact(name):
+    img = _images(3)[name]
+    # seed must be <= mask; use an erosion of the mask as the seed (the granularity pattern)
+    mask = img
+    seed = M.erosion(mask, M.disk(1))
+    expected = M.reconstruction(seed, mask, footprint=M.disk(1))
+    got = reconstruction_by_dilation_2d(
+        seed.astype(numpy.float64), mask.astype(numpy.float64)
+    )
+    numpy.testing.assert_array_equal(got, expected)
+
+
+def test_reconstruction_cascaded_mask_chain():
+    # rec_g uses rec_{g-1} as mask (cascade); each step must still match skimage
+    rng = numpy.random.default_rng(7)
+    pixels = rng.random((25, 25))
+    ero = pixels.copy()
+    recon_mask = pixels
+    for _ in range(5):
+        ero = erosion_4conn_2d(ero)
+        rec = reconstruction_by_dilation_2d(ero, recon_mask)
+        expected = M.reconstruction(ero, recon_mask, footprint=M.disk(1))
+        numpy.testing.assert_array_equal(rec, expected)
+        recon_mask = rec
+
+
+def test_reconstruction_long_geodesic():
+    """A serpentine channel forces propagation along a long geodesic path — the
+    case the multi-raster/FIFO split is built for. The result is independent of how
+    many raster passes precede the FIFO, so it must stay bit-exact vs skimage."""
+    H, W = 41, 41
+    mask = numpy.zeros((H, W))
+    # carve a snake of open (high) corridors one row apart, connected at alternating ends
+    mask[1:-1, 1:-1] = 9.0
+    for i in range(2, H - 2, 2):
+        if (i // 2) % 2 == 0:
+            mask[i, 1:-2] = 0.0  # wall from the left, gap on the right
+        else:
+            mask[i, 2:-1] = 0.0  # wall from the right, gap on the left
+    seed = numpy.zeros((H, W))
+    seed[1, 1] = 9.0  # single high seed at one end of the snake
+    got = reconstruction_by_dilation_2d(seed, mask)
+    expected = M.reconstruction(seed, mask, footprint=M.disk(1))
+    numpy.testing.assert_array_equal(got, expected)
diff --git a/test/test_primitives_shapes.py b/test/test_primitives_shapes.py
new file mode 100644
index 0000000..0584e41
--- /dev/null
+++ b/test/test_primitives_shapes.py
@@ -0,0 +1,80 @@
+"""Tests for the canonical (B,Z,Y,X) input normalisation helper."""
+
+import numpy
+import pytest
+
+from cp_measure.primitives.shapes import to_bzyx
+
+
+def _img(shape, seed):
+    rng = numpy.random.default_rng(seed)
+    return rng.random(shape)
+
+
+def test_2d_is_single_unit_z():
+    m, p = numpy.ones((4, 5), int), _img((4, 5), 0)
+    masks, pixels, unwrap = to_bzyx(m, p)
+    assert len(masks) == len(pixels) == 1
+    assert masks[0].shape == (1, 4, 5)
+    assert pixels[0].shape == (1, 4, 5)
+    # single input -> unwrap returns the lone result, not a list
+    assert unwrap(["only"]) == "only"
+
+
+def test_3d_is_single_volume_not_batch():
+    m, p = numpy.ones((3, 4, 5), int), _img((3, 4, 5), 1)
+    masks, pixels, unwrap = to_bzyx(m, p)
+    assert len(masks) == 1
+    assert masks[0].shape == (3, 4, 5)  # one volume, Z preserved
+    assert unwrap(["d"]) == "d"
+
+
+def test_4d_is_batch():
+    m, p = numpy.ones((2, 3, 4, 5), int), _img((2, 3, 4, 5), 2)
+    masks, pixels, unwrap = to_bzyx(m, p)
+    assert len(masks) == len(pixels) == 2
+    assert all(a.shape == (3, 4, 5) for a in masks)
+    out = unwrap([{"a": 1}, {"a": 2}])
+    assert isinstance(out, list) and len(out) == 2  # batch -> list
+
+
+def test_list_of_2d_is_batch_unit_z():
+    imgs = [_img((4, 5), i) for i in range(3)]
+    masks = [numpy.ones((4, 5), int) for _ in range(3)]
+    m, p, unwrap = to_bzyx(masks, imgs)
+    assert len(m) == 3
+    assert all(a.shape == (1, 4, 5) for a in m)
+    assert isinstance(unwrap([1, 2, 3]), list)
+
+
+def test_list_ragged_sizes_ok():
+    imgs = [_img((4, 5), 0), _img((7, 3), 1)]
+    masks = [numpy.ones((4, 5), int), numpy.ones((7, 3), int)]
+    m, p, _ = to_bzyx(masks, imgs)
+    assert m[0].shape == (1, 4, 5)
+    assert m[1].shape == (1, 7, 3)
+
+
+def test_list_of_3d_volumes_batch():
+    vols = [_img((2, 4, 5), 0), _img((3, 4, 5), 1)]
+    masks = [numpy.ones((2, 4, 5), int), numpy.ones((3, 4, 5), int)]
+    m, _, _ = to_bzyx(masks, vols)
+    assert m[0].shape == (2, 4, 5) and m[1].shape == (3, 4, 5)
+
+
+@pytest.mark.parametrize(
+    "masks, pixels, match",
+    [
+        (
+            [numpy.ones((4, 5), int)],
+            [_img((4, 5), 0), _img((4, 5), 1)],
+            "batch size mismatch",
+        ),
+        ([numpy.ones((4, 5), int)], _img((4, 5), 0), "both be sequences"),
+        (numpy.ones((2, 3, 4, 5), int), _img((3, 4, 5), 0), "both be 4D"),
+        (numpy.ones((5,), int), numpy.ones((5,), float), "2D or 3D"),
+    ],
+)
+def test_to_bzyx_raises(masks, pixels, match):
+    with pytest.raises(ValueError, match=match):
+        to_bzyx(masks, pixels)
diff --git a/uv.lock b/uv.lock
index da478cd..9b66ae0 100644
--- a/uv.lock
+++ b/uv.lock
@@ -537,6 +537,9 @@ benchmark = [
     { name = "tifffile", version = "2025.5.10", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },
     { name = "tifffile", version = "2025.10.16", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" },
 ]
+numba = [
+    { name = "numba" },
+]
 test = [
     { name = "pytest" },
     { name = "pytest-cov" },
@@ -558,6 +561,7 @@ dev = [
 requires-dist = [
     { name = "centrosome", specifier = ">=1.3.3" },
     { name = "mahotas", specifier = ">=1.4.13,<2.0.0" },
+    { name = "numba", marker = "extra == 'numba'", specifier = ">=0.59" },
     { name = "numpy", specifier = ">=1.22.1" },
     { name = "pyarrow", marker = "extra == 'benchmark'", specifier = ">=14.0" },
     { name = "pytest", marker = "extra == 'test'", specifier = ">=8.4.2" },
@@ -567,7 +571,7 @@ requires-dist = [
     { name = "scipy", specifier = ">=1.9.1" },
     { name = "tifffile", marker = "extra == 'benchmark'", specifier = ">=2025.5.10" },
 ]
-provides-extras = ["test", "benchmark"]
+provides-extras = ["numba", "test", "benchmark"]
 
 [package.metadata.requires-dev]
 dev = [
@@ -1324,6 +1328,38 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/aa/47/7d70414bcdbb3bc1f458a8d10558f00bbfdb24e5a11740fc8197e12c3255/librt-0.9.0-cp314-cp314t-win_arm64.whl", hash = "sha256:a4b25c6c25cac5d0d9d6d6da855195b254e0021e513e0249f0e3b444dc6e0e61", size = 50009, upload-time = "2026-04-09T16:06:07.995Z" },
 ]
 
+[[package]]
+name = "llvmlite"
+version = "0.47.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/01/88/a8952b6d5c21e74cbf158515b779666f692846502623e9e3c39d8e8ba25f/llvmlite-0.47.0.tar.gz", hash = "sha256:62031ce968ec74e95092184d4b0e857e444f8fdff0b8f9213707699570c33ccc", size = 193614, upload-time = "2026-03-31T18:29:53.497Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/f4/f5/a1bde3aa8c43524b0acaf3f72fb3d80a32dd29dbb42d7dc434f84584cdcc/llvmlite-0.47.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:41270b0b1310717f717cf6f2a9c68d3c43bd7905c33f003825aebc361d0d1b17", size = 37232772, upload-time = "2026-03-31T18:28:12.198Z" },
+    { url = "https://files.pythonhosted.org/packages/7c/fb/76d88fc05ee1f9c1a6efe39eb493c4a727e5d1690412469017cd23bcb776/llvmlite-0.47.0-cp310-cp310-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:f9d118bc1dd7623e0e65ca9ac485ec6dd543c3b77bc9928ddc45ebd34e1e30a7", size = 56275179, upload-time = "2026-03-31T18:28:15.725Z" },
+    { url = "https://files.pythonhosted.org/packages/4d/08/29da7f36217abd56a0c389ef9a18bea47960826e691ced1a36c92c6ce93c/llvmlite-0.47.0-cp310-cp310-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:9ea5cfb04a6ab5b18e46be72b41b015975ba5980c4ddb41f1975b83e19031063", size = 55128632, upload-time = "2026-03-31T18:28:19.946Z" },
+    { url = "https://files.pythonhosted.org/packages/df/f8/5e12e9ed447d65f04acf6fcf2d79cded2355640b5131a46cee4c99a5949d/llvmlite-0.47.0-cp310-cp310-win_amd64.whl", hash = "sha256:166b896a2262a2039d5fc52df5ee1659bd1ccd081183df7a2fba1b74702dd5ea", size = 38138402, upload-time = "2026-03-31T18:28:23.327Z" },
+    { url = "https://files.pythonhosted.org/packages/34/0b/b9d1911cfefa61399821dfb37f486d83e0f42630a8d12f7194270c417002/llvmlite-0.47.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:74090f0dcfd6f24ebbef3f21f11e38111c4d7e6919b54c4416e1e357c3446b07", size = 37232770, upload-time = "2026-03-31T18:28:26.765Z" },
+    { url = "https://files.pythonhosted.org/packages/46/27/5799b020e4cdfb25a7c951c06a96397c135efcdc21b78d853bbd9c814c7d/llvmlite-0.47.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:ca14f02e29134e837982497959a8e2193d6035235de1cb41a9cb2bd6da4eedbb", size = 56275177, upload-time = "2026-03-31T18:28:31.01Z" },
+    { url = "https://files.pythonhosted.org/packages/7e/51/48a53fedf01cb1f3f43ef200be17ebf83c8d9a04018d3783c1a226c342c2/llvmlite-0.47.0-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:12a69d4bb05f402f30477e21eeabe81911e7c251cecb192bed82cd83c9db10d8", size = 55128631, upload-time = "2026-03-31T18:28:36.046Z" },
+    { url = "https://files.pythonhosted.org/packages/a2/50/59227d06bdc96e23322713c381af4e77420949d8cd8a042c79e0043096cc/llvmlite-0.47.0-cp311-cp311-win_amd64.whl", hash = "sha256:c37d6eb7aaabfa83ab9c2ff5b5cdb95a5e6830403937b2c588b7490724e05327", size = 38138400, upload-time = "2026-03-31T18:28:40.076Z" },
+    { url = "https://files.pythonhosted.org/packages/fa/48/4b7fe0e34c169fa2f12532916133e0b219d2823b540733651b34fdac509a/llvmlite-0.47.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:306a265f408c259067257a732c8e159284334018b4083a9e35f67d19792b164f", size = 37232769, upload-time = "2026-03-31T18:28:43.735Z" },
+    { url = "https://files.pythonhosted.org/packages/e6/4b/e3f2cd17822cf772a4a51a0a8080b0032e6d37b2dbe8cfb724eac4e31c52/llvmlite-0.47.0-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:5853bf26160857c0c2573415ff4efe01c4c651e59e2c55c2a088740acfee51cd", size = 56275178, upload-time = "2026-03-31T18:28:48.342Z" },
+    { url = "https://files.pythonhosted.org/packages/b6/55/a3b4a543185305a9bdf3d9759d53646ed96e55e7dfd43f53e7a421b8fbae/llvmlite-0.47.0-cp312-cp312-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:003bcf7fa579e14db59c1a1e113f93ab8a06b56a4be31c7f08264d1d4072d077", size = 55128632, upload-time = "2026-03-31T18:28:52.901Z" },
+    { url = "https://files.pythonhosted.org/packages/2f/f5/d281ae0f79378a5a91f308ea9fdb9f9cc068fddd09629edc0725a5a8fde1/llvmlite-0.47.0-cp312-cp312-win_amd64.whl", hash = "sha256:f3079f25bdc24cd9d27c4b2b5e68f5f60c4fdb7e8ad5ee2b9b006007558f9df7", size = 38138692, upload-time = "2026-03-31T18:28:57.147Z" },
+    { url = "https://files.pythonhosted.org/packages/77/6f/4615353e016799f80fa52ccb270a843c413b22361fadda2589b2922fb9b0/llvmlite-0.47.0-cp313-cp313-macosx_12_0_arm64.whl", hash = "sha256:a3c6a735d4e1041808434f9d440faa3d78d9b4af2ee64d05a66f351883b6ceec", size = 37232771, upload-time = "2026-03-31T18:29:01.324Z" },
+    { url = "https://files.pythonhosted.org/packages/31/b8/69f5565f1a280d032525878a86511eebed0645818492feeb169dfb20ae8e/llvmlite-0.47.0-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:2699a74321189e812d476a43d6d7f652f51811e7b5aad9d9bba842a1c7927acb", size = 56275178, upload-time = "2026-03-31T18:29:05.748Z" },
+    { url = "https://files.pythonhosted.org/packages/d6/da/b32cafcb926fb0ce2aa25553bf32cb8764af31438f40e2481df08884c947/llvmlite-0.47.0-cp313-cp313-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:6c6951e2b29930227963e53ee152441f0e14be92e9d4231852102d986c761e40", size = 55128632, upload-time = "2026-03-31T18:29:11.235Z" },
+    { url = "https://files.pythonhosted.org/packages/46/9f/4898b44e4042c60fafcb1162dfb7014f6f15b1ec19bf29cfea6bf26df90d/llvmlite-0.47.0-cp313-cp313-win_amd64.whl", hash = "sha256:c2e9adf8698d813a9a5efb2d4370caf344dbc1e145019851fee6a6f319ba760e", size = 38138695, upload-time = "2026-03-31T18:29:15.43Z" },
+    { url = "https://files.pythonhosted.org/packages/1c/d4/33c8af00f0bf6f552d74f3a054f648af2c5bc6bece97972f3bfadce4f5ec/llvmlite-0.47.0-cp314-cp314-macosx_12_0_arm64.whl", hash = "sha256:de966c626c35c9dff5ae7bf12db25637738d0df83fc370cf793bc94d43d92d14", size = 37232773, upload-time = "2026-03-31T18:29:19.453Z" },
+    { url = "https://files.pythonhosted.org/packages/64/1d/a760e993e0c0ba6db38d46b9f48f6c7dceb8ac838824997fb9e25f97bc04/llvmlite-0.47.0-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:ddbccff2aeaff8670368340a158abefc032fe9b3ccf7d9c496639263d00151aa", size = 56275176, upload-time = "2026-03-31T18:29:24.149Z" },
+    { url = "https://files.pythonhosted.org/packages/84/3b/e679bc3b29127182a7f4aa2d2e9e5bea42adb93fb840484147d59c236299/llvmlite-0.47.0-cp314-cp314-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:d4a7b778a2e144fc64468fb9bf509ac1226c9813a00b4d7afea5d988c4e22fca", size = 55128631, upload-time = "2026-03-31T18:29:29.536Z" },
+    { url = "https://files.pythonhosted.org/packages/be/f7/19e2a09c62809c9e63bbd14ce71fb92c6ff7b7b3045741bb00c781efc3c9/llvmlite-0.47.0-cp314-cp314-win_amd64.whl", hash = "sha256:694e3c2cdc472ed2bd8bd4555ca002eec4310961dd58ef791d508f57b5cc4c94", size = 39153826, upload-time = "2026-03-31T18:29:33.681Z" },
+    { url = "https://files.pythonhosted.org/packages/40/a1/581a8c707b5e80efdbbe1dd94527404d33fe50bceb71f39d5a7e11bd57b7/llvmlite-0.47.0-cp314-cp314t-macosx_12_0_arm64.whl", hash = "sha256:92ec8a169a20b473c1c54d4695e371bde36489fc1efa3688e11e99beba0abf9c", size = 37232772, upload-time = "2026-03-31T18:29:37.952Z" },
+    { url = "https://files.pythonhosted.org/packages/11/03/16090dd6f74ba2b8b922276047f15962fbeea0a75d5601607edb301ba945/llvmlite-0.47.0-cp314-cp314t-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:fa1cbd800edd3b20bc141521f7fd45a6185a5b84109aa6855134e81397ffe72b", size = 56275178, upload-time = "2026-03-31T18:29:42.58Z" },
+    { url = "https://files.pythonhosted.org/packages/f5/cb/0abf1dd4c5286a95ffe0c1d8c67aec06b515894a0dd2ac97f5e27b82ab0b/llvmlite-0.47.0-cp314-cp314t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:f6725179b89f03b17dabe236ff3422cb8291b4c1bf40af152826dfd34e350ae8", size = 55128632, upload-time = "2026-03-31T18:29:46.939Z" },
+    { url = "https://files.pythonhosted.org/packages/4f/79/d3bbab197e86e0ff4f9c07122895b66a3e0d024247fcff7f12c473cb36d9/llvmlite-0.47.0-cp314-cp314t-win_amd64.whl", hash = "sha256:6842cf6f707ec4be3d985a385ad03f72b2d724439e118fcbe99b2929964f0453", size = 39153839, upload-time = "2026-03-31T18:29:51.004Z" },
+]
+
 [[package]]
 name = "mahotas"
 version = "1.4.18"
@@ -1683,6 +1719,42 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/7f/95/4df134a100b5a9a12378d5301b934366686ef6fbdaffcd21211d5654970e/nox-2026.4.10-py3-none-any.whl", hash = "sha256:082c117627590d9b90aa21f86df89b310b07c5842539524203bcb3c719f116c1", size = 75536, upload-time = "2026-04-10T17:42:40.664Z" },
 ]
 
+[[package]]
+name = "numba"
+version = "0.65.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "llvmlite" },
+    { name = "numpy" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/f6/c5/db2ac3685833d626c0dcae6bd2330cd68433e1fd248d15f70998160d3ad7/numba-0.65.1.tar.gz", hash = "sha256:19357146c32fe9ed25059ab915e8465fb13951cf6b0aace3826b76886373ab23", size = 2765600, upload-time = "2026-04-24T02:02:56.551Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/de/1b/3c5a7daf683a95465bf23504bcd1a2d5db8cd5e5e276ca87505d020dffe9/numba-0.65.1-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:9d993ed0a257aa4116e6f553f114004bcfdee540c7276ab8ea48f650d514c452", size = 2680870, upload-time = "2026-04-24T02:02:10.623Z" },
+    { url = "https://files.pythonhosted.org/packages/0f/a4/1831836814018a898e7d252aebe09c0f3ce1f26d145b68264b4ae0be6822/numba-0.65.1-cp310-cp310-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:5f098109f361681e57295f7e84d8ab2426902539a141811de0703ace52826981", size = 3739780, upload-time = "2026-04-24T02:02:13.097Z" },
+    { url = "https://files.pythonhosted.org/packages/9c/1b/a813ddc81def09e257d2b1f67521982ce4b06204a87268796ffc8187271c/numba-0.65.1-cp310-cp310-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:973fd8173f2312815e6b7aaae887c4ce8a817eeff46a4f8840b828305b75bc95", size = 3446722, upload-time = "2026-04-24T02:02:15.083Z" },
+    { url = "https://files.pythonhosted.org/packages/09/52/ee1d8b3becda384fe0552221641e05aa668a35e8a77470db4db7f6475000/numba-0.65.1-cp310-cp310-win_amd64.whl", hash = "sha256:c63aa0c4193694026452da55d0ef9d85156c1a7a333454c103bb30dec81b7bf8", size = 2747539, upload-time = "2026-04-24T02:02:16.79Z" },
+    { url = "https://files.pythonhosted.org/packages/96/b3/650500c2eab4534d98e9166f4298e0f3c69c742afdf24e6eabccd1f16ad8/numba-0.65.1-cp311-cp311-macosx_12_0_arm64.whl", hash = "sha256:7020d74b19cdb8cff16506542fdd510756e28c5e7f3bd0b7f574f0f42272fcd9", size = 2680563, upload-time = "2026-04-24T02:02:18.414Z" },
+    { url = "https://files.pythonhosted.org/packages/44/0b/0615dbedb98f5b32a35a53290fbdc6e22306968109278d7e58df82d7a9f6/numba-0.65.1-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:f80ed83774b5173abd6581cd8d2165d1d38e13d2e5c8155c0c0b421784745420", size = 3745018, upload-time = "2026-04-24T02:02:20.252Z" },
+    { url = "https://files.pythonhosted.org/packages/49/aa/4361698f35bf63bff67dfe6c90493731177f48ede954f77b0588731537bc/numba-0.65.1-cp311-cp311-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:7ed425a43b0a5f9772f2f4e2dd0bbd12eabecae1af0b24efcfd4e053f012aac6", size = 3450962, upload-time = "2026-04-24T02:02:22.449Z" },
+    { url = "https://files.pythonhosted.org/packages/bd/9a/af61ec03b3116c161fd7a06b9e8a265729a8718458333e8ffbb06d9a3978/numba-0.65.1-cp311-cp311-win_amd64.whl", hash = "sha256:df40a5028a975b9ea66f6a2a3f7abbdbd541a863070e34ed367aff21141248e4", size = 2747417, upload-time = "2026-04-24T02:02:24.43Z" },
+    { url = "https://files.pythonhosted.org/packages/57/bc/76f8f8c5cf9adee47fdb7bbb03be8900f76f902d451d7477cf12b845e1de/numba-0.65.1-cp312-cp312-macosx_12_0_arm64.whl", hash = "sha256:ac3f1e77c352dd0ea9712732c2d8f9ca507717435eec5b5013bf138ac33c4a08", size = 2681371, upload-time = "2026-04-24T02:02:26.105Z" },
+    { url = "https://files.pythonhosted.org/packages/69/47/a415af0283e4db0398104c6d1c11c9861a98dc67a7aa442a7769ed5d6196/numba-0.65.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:52bc6f3ceb8fcaff9b2ae26b4c6b1e9fee39db8d355534c0fe4f39a901246b84", size = 3802467, upload-time = "2026-04-24T02:02:27.712Z" },
+    { url = "https://files.pythonhosted.org/packages/46/36/246f73ec99cfeab2f2cb2ce7d4218766cc36a2da418901223f4f4da9c813/numba-0.65.1-cp312-cp312-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:90ca10b3463bae0bd70589726fe3c77d01d6b5fc86bee54bcdf9fb6b47c28977", size = 3502628, upload-time = "2026-04-24T02:02:29.763Z" },
+    { url = "https://files.pythonhosted.org/packages/db/9e/3c679b2ee078425b9e99a91e44f8d132a6830d8ccce5227bc5e9181aeed8/numba-0.65.1-cp312-cp312-win_amd64.whl", hash = "sha256:5971c632be2a2351500431f46213821dba8d02b18a9f7d02fd36bd2743e41a6a", size = 2750611, upload-time = "2026-04-24T02:02:31.477Z" },
+    { url = "https://files.pythonhosted.org/packages/79/37/14a4579049c1eb673afd0de0cb4842982acd55b9ce2643e763db858bcea0/numba-0.65.1-cp313-cp313-macosx_12_0_arm64.whl", hash = "sha256:1735c15c1134a5108b4d6a5c77fc0947924ea066a738dc09a52008c13df9cad3", size = 2681344, upload-time = "2026-04-24T02:02:33.65Z" },
+    { url = "https://files.pythonhosted.org/packages/a0/22/b8d873f6466b20aa563fc9b33acd48dec89a07803ddaa2f1c8ca1cd33126/numba-0.65.1-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:c09f49117ef255e1f1c6dad0c7a1ed39868243862a73be5706793241a3755f1b", size = 3810619, upload-time = "2026-04-24T02:02:36.041Z" },
+    { url = "https://files.pythonhosted.org/packages/62/08/e16a8b5d9a018962ebb5c66be662317cde32b9f5dab08441f90bed5522fb/numba-0.65.1-cp313-cp313-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:594a8680b3fadac99e97e489b1fd89007177e5336713745c3b769528c635a464", size = 3509783, upload-time = "2026-04-24T02:02:38.245Z" },
+    { url = "https://files.pythonhosted.org/packages/fd/a5/03c970d57f4c1741354837353ce39fb5206952ae1dba8922d29c86f64805/numba-0.65.1-cp313-cp313-win_amd64.whl", hash = "sha256:85be74c0d036842699a30058f82fb88fc5ffdc59f7615cab5792ea92914c9b62", size = 2750534, upload-time = "2026-04-24T02:02:39.903Z" },
+    { url = "https://files.pythonhosted.org/packages/4f/2e/8aed9b726d9ba5f11ad287645fd479e88278db3060a25cb1225d730eb2b7/numba-0.65.1-cp314-cp314-macosx_12_0_arm64.whl", hash = "sha256:33f5eb68eb1c843511615d14663ce60258525d6a4c65ab040e2c2b0c4cf17450", size = 2681554, upload-time = "2026-04-24T02:02:41.812Z" },
+    { url = "https://files.pythonhosted.org/packages/87/96/f3eb235fafa82a34e2ab5dd7dc9ffff998ebf5f0bbc23fa56a96aeb44da6/numba-0.65.1-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:71e73029bf53a62cc6afcf96be4bd942290d8b4c55f0a454fb536158115790f7", size = 3779602, upload-time = "2026-04-24T02:02:43.726Z" },
+    { url = "https://files.pythonhosted.org/packages/09/90/b0f09b48752d23640b8284f22aa597737e8adaddc7fbfacc4708b7f73a4c/numba-0.65.1-cp314-cp314-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:3a07635e0be926b9bdbffb09137c230fb13f6ec0e564914ba937cee12ce3eb35", size = 3479532, upload-time = "2026-04-24T02:02:45.427Z" },
+    { url = "https://files.pythonhosted.org/packages/56/46/3f7fc04fb853559e74b210e0b62c19974ec844cefec611f9e535f4da3761/numba-0.65.1-cp314-cp314-win_amd64.whl", hash = "sha256:2a20fcdabdefbdacf88d85caf70c3b18c4bcb7ebb8f82e6a19486383dd26ab63", size = 2752637, upload-time = "2026-04-24T02:02:47.664Z" },
+    { url = "https://files.pythonhosted.org/packages/81/7b/c1a341a9067367778f4152a5f01061cf281fb09582c92c510ec4918cabf6/numba-0.65.1-cp314-cp314t-macosx_12_0_arm64.whl", hash = "sha256:548dd4b3a4508d5062768d1514b2cd7b015f9a25ec7af651c50dee243965e652", size = 2684600, upload-time = "2026-04-24T02:02:49.653Z" },
+    { url = "https://files.pythonhosted.org/packages/03/36/98ddbcf3e4f04a6dd07e1c67249955920579ba4af6bb6868e3088f4ed282/numba-0.65.1-cp314-cp314t-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:78abc28feff2c2ff8307fff3975b6438352759c9acb797ecd6b1fb6e7e39e31d", size = 3817198, upload-time = "2026-04-24T02:02:51.266Z" },
+    { url = "https://files.pythonhosted.org/packages/a3/83/0dad21057ece5a835599f5d24099b091703995e23dbbf894f259e91c010b/numba-0.65.1-cp314-cp314t-manylinux_2_27_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:ee7676cb389555805f9b9a1840cbcd1ea6c8bd5376ab6918e3a29c5ea1dbda20", size = 3533862, upload-time = "2026-04-24T02:02:52.987Z" },
+    { url = "https://files.pythonhosted.org/packages/32/36/8be7118ffd4c8440881046eac3d0982cc5ab42909508cf5d67024d62a2e4/numba-0.65.1-cp314-cp314t-win_amd64.whl", hash = "sha256:20609346e3bd75204950dcbbfe383a8d7dbf4902f442aedbf00f97fef4aa8f38", size = 2758237, upload-time = "2026-04-24T02:02:54.612Z" },
+]
+
 [[package]]
 name = "numpy"
 version = "1.26.4"