From 913e3d388d7225e664e902096951d8010409feba Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 04:47:44 +0200 Subject: [PATCH 01/15] feat(accelerator): wire numba intensity backend on the #49 dispatch seam MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit First real accelerator end-to-end on top of the merged #49 dispatch: `set_accelerator("numba")` now routes `intensity` to a numba implementation and composes it with the numpy backend for every other feature. - _detect.py: capability flags (HAS_NUMBA/HAS_JAX/HAS_JAX_GPU) via find_spec, resolved once at import. No try/except — an absent backend is never attempted, a present-but-broken one raises. - primitives/: shared host segment layer. flatten_labeled reduces a labeled (Z,Y,X) image to flat (values, seg0, coords); a single kernel set then covers 2D, 3D and future batches with no image/batch axis baked in. max_position is a host scipy.ndimage.maximum_position call for bit-exact parity with the numpy backend's tie-break. - primitives/_segment_numba.py: @njit(cache=True), single-threaded kernels — fused single-pass moments + centroid cross-sums, residual-sumsq std, CSR per-segment quantiles/MAD. - core/numba/: import-selected backend (`from cp_measure.core.numba import get_intensity`); identical dict contract, 2D and 3D. - bulk._dispatch: "numba" composes numba intensity + numpy rest; raises if numba is not installed (no silent fallback). - numba is an optional extra ([numba]); the default install stays numba-free. CI tests install .[numba] and run the correctness harness. test/test_backend_correctness.py asserts numba == numpy (2D/3D, edge on/off, rtol=1e-6), the dispatch composition, and the absent-numba raise path. Co-Authored-By: Claude Opus 4.8 (1M context) --- noxfile.py | 2 +- pyproject.toml | 3 + src/cp_measure/_detect.py | 24 +++ src/cp_measure/bulk.py | 28 ++- src/cp_measure/core/numba/__init__.py | 14 ++ .../core/numba/measureobjectintensity.py | 185 ++++++++++++++++++ src/cp_measure/primitives/__init__.py | 22 +++ src/cp_measure/primitives/_segment_numba.py | 126 ++++++++++++ src/cp_measure/primitives/segment.py | 109 +++++++++++ test/test_backend_correctness.py | 90 +++++++++ 10 files changed, 599 insertions(+), 4 deletions(-) create mode 100644 src/cp_measure/_detect.py create mode 100644 src/cp_measure/core/numba/__init__.py create mode 100644 src/cp_measure/core/numba/measureobjectintensity.py create mode 100644 src/cp_measure/primitives/__init__.py create mode 100644 src/cp_measure/primitives/_segment_numba.py create mode 100644 src/cp_measure/primitives/segment.py create mode 100644 test/test_backend_correctness.py 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..d9ca785 --- /dev/null +++ b/src/cp_measure/_detect.py @@ -0,0 +1,24 @@ +"""Backend capability detection. + +Detected ONCE at import via ``importlib.util.find_spec`` — availability is +checked without importing the package or catching ImportErrors. Dispatch reads +these flags; 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). +""" + +import importlib.util + +HAS_NUMBA: bool = importlib.util.find_spec("numba") is not None +HAS_JAX: bool = importlib.util.find_spec("jax") is not None + + +def _detect_jax_gpu() -> bool: + if not HAS_JAX: + return False + import jax + + return jax.default_backend() != "cpu" + + +HAS_JAX_GPU: bool = _detect_jax_gpu() diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py index 1f2dc94..1824431 100644 --- a/src/cp_measure/bulk.py +++ b/src/cp_measure/bulk.py @@ -45,6 +45,23 @@ _3D_FEATURES = ("intensity", "sizeshape", "texture", "granularity") +def _numba_registries() -> dict[str, dict[str, Callable]]: + """Registries for the 'numba' accelerator. + + Composes the numba implementations (currently ``intensity`` only) with the + numpy implementations of every other feature — a single global "numba" + selection still yields a full, working feature set, accelerated where a + numba backend exists. This is explicit per-function composition, NOT an + error-driven fallback. + """ + from cp_measure.core.numba import get_intensity as _numba_intensity + + return { + "core": {**_CORE, "intensity": _numba_intensity}, + "correlation": _CORRELATION, + } + + def _dispatch(name: str) -> dict[str, Callable]: from cp_measure import _ACCELERATOR @@ -55,9 +72,14 @@ def _dispatch(name: str) -> dict[str, Callable]: f"'jax' accelerator not yet wired for {name} measurements" ) if _ACCELERATOR == "numba": - raise NotImplementedError( - f"'numba' accelerator not yet wired for {name} measurements" - ) + from cp_measure._detect import HAS_NUMBA + + if not HAS_NUMBA: + raise RuntimeError( + "accelerator 'numba' selected but numba is not installed; " + "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..671d723 --- /dev/null +++ b/src/cp_measure/core/numba/__init__.py @@ -0,0 +1,14 @@ +"""Numba-accelerated backend. + +Selected explicitly by import (``from cp_measure.core.numba import get_intensity``) +or globally via ``cp_measure.set_accelerator("numba")``. Requires the optional +``numba`` extra; availability is gated by ``cp_measure._detect.HAS_NUMBA``. + +This backend currently accelerates ``intensity`` only; the global "numba" +accelerator composes it with the numpy implementations of every other feature +(see ``cp_measure.bulk``). +""" + +from cp_measure.core.numba.measureobjectintensity import get_intensity + +__all__ = ["get_intensity"] diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py new file mode 100644 index 0000000..64580f1 --- /dev/null +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -0,0 +1,185 @@ +"""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 run as fused single-pass numba kernels over a flat segment +representation (:mod:`cp_measure.primitives`); ``Location_MaxIntensity_*`` is +computed on the host via ``scipy.ndimage.maximum_position`` to reproduce the +numpy backend's tie-break exactly. +""" + +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 import ( + flatten_labeled, + label_to_idx_lut, + max_position_per_object, +) +from cp_measure.primitives._segment_numba import ( + segment_moments, + segment_quantiles, + segment_resid_sumsq, +) + + +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, labels, 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_labeled(masks, masked_image, lut) + has_objects = values.size > 0 + + if has_objects: + (count, sumI, minI, maxI, 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) + + max_x, max_y, max_z = max_position_per_object( + values, seg0, xc, yc, zc, labels + ) + + 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) + + outlines = skimage.segmentation.find_boundaries(masks, mode="inner") + emask = outlines > 0 + e_values = masked_image[emask].astype(numpy.float64) + e_seg0 = lut[masks[emask]] + + if e_values.size > 0: + zeros = numpy.zeros(e_values.size) + (ecount, esum, emin, emax, *_rest) = segment_moments( + e_values, e_seg0, zeros, zeros, zeros, 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..673b714 --- /dev/null +++ b/src/cp_measure/primitives/__init__.py @@ -0,0 +1,22 @@ +"""Shared primitive layer. + +Backend-agnostic building blocks that the per-backend feature implementations +(``cp_measure.core`` = numpy, ``cp_measure.core.numba`` = numba, ...) compose. + +The host helpers here flatten a labeled image into a 1-D *segment* representation +(values + 0-based segment index + per-axis coordinates). 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. +""" + +from cp_measure.primitives.segment import ( + flatten_labeled, + label_to_idx_lut, + max_position_per_object, +) + +__all__ = [ + "flatten_labeled", + "label_to_idx_lut", + "max_position_per_object", +] diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py new file mode 100644 index 0000000..4c236a3 --- /dev/null +++ b/src/cp_measure/primitives/_segment_numba.py @@ -0,0 +1,126 @@ +"""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 segment_moments(values, seg0, xc, yc, zc, n): + """One pass over the flat pixels accumulating, per segment: + + count, sum, min, max, and the six centroid cross-sums + (Sx, Sy, Sz, Sx*I, Sy*I, Sz*I). max position is intentionally NOT computed + here (it is done on the host via scipy for exact parity). + """ + count = np.zeros(n, np.int64) + sumI = np.zeros(n, np.float64) + minI = np.full(n, np.inf) + maxI = np.full(n, -np.inf) + 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] + count[k] += 1 + sumI[k] += v + if v < minI[k]: + minI[k] = v + if v > maxI[k]: + maxI[k] = v + x = xc[i] + y = yc[i] + z = zc[i] + 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, sx, sy, sz, sxI, syI, szI + + +@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..c7dcdf8 --- /dev/null +++ b/src/cp_measure/primitives/segment.py @@ -0,0 +1,109 @@ +"""Host-side segment helpers (numpy/scipy, backend-agnostic). + +A labeled image is reduced to a flat *segment* form: + +* ``values[M]`` finite labeled pixel intensities (float64) +* ``seg0[M]`` 0-based segment index = rank of the pixel's label among the + sorted positive labels (``label_to_idx``) +* ``xc/yc/zc[M]``per-pixel coordinates (z = 0 plane for 2D input) + +Segment kernels (numpy or numba) then loop over these flat M-length arrays with +no notion of image shape or batch axis. 2D vs 3D differ only in which coordinate +arrays are populated; a batch differs only in how ``seg0`` offsets are assigned. +""" + +import numpy +import scipy.ndimage +from numpy.typing import NDArray + + +def label_to_idx_lut( + masks: NDArray[numpy.integer], +) -> tuple[NDArray[numpy.int64], NDArray[numpy.integer], int]: + """Build a ``label -> 0..n-1`` lookup over the sorted positive labels. + + Returns ``(lut, labels_sorted, n)`` where ``lut[label]`` is the segment + index (and ``-1`` for absent labels / background), ``labels_sorted`` are the + ascending positive labels, 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. + """ + unique = numpy.unique(masks) + labels = unique[unique > 0] + 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, labels, n + + +def flatten_labeled( + masks: NDArray[numpy.integer], + pixels: NDArray[numpy.floating], + lut: NDArray[numpy.int64], +) -> tuple[ + NDArray[numpy.float64], + NDArray[numpy.int64], + NDArray[numpy.float64], + NDArray[numpy.float64], + NDArray[numpy.float64], +]: + """Flatten a labeled (Z, Y, X) image to ``(values, seg0, xc, yc, zc)``. + + Non-finite pixels and background are dropped, matching the numpy reference's + ``(masks > 0) & isfinite(pixels)`` mask. Coordinates use ``numpy.mgrid`` over + the full volume exactly as the reference does, so positions line up. + """ + lmask = (masks > 0) & numpy.isfinite(pixels) + values = pixels[lmask].astype(numpy.float64) + seg0 = lut[masks[lmask]] + mesh_z, mesh_y, mesh_x = numpy.mgrid[ + 0 : masks.shape[0], 0 : masks.shape[1], 0 : masks.shape[2] + ] + xc = mesh_x[lmask].astype(numpy.float64) + yc = mesh_y[lmask].astype(numpy.float64) + zc = mesh_z[lmask].astype(numpy.float64) + return values, seg0, xc, yc, zc + + +def max_position_per_object( + values: NDArray[numpy.float64], + seg0: NDArray[numpy.int64], + xc: NDArray[numpy.float64], + yc: NDArray[numpy.float64], + zc: NDArray[numpy.float64], + labels_sorted: NDArray[numpy.integer], +) -> tuple[ + NDArray[numpy.float64], NDArray[numpy.float64], NDArray[numpy.float64] +]: + """Per-object argmax position via ``scipy.ndimage.maximum_position``. + + Reproduces the shipped numpy implementation EXACTLY, including scipy's + labeled tie-break on exact-equal maxima (an implementation artifact that is + neither reliably first nor last). The reference calls ``maximum_position`` + once per object on that object's pixels in raster order; this does the same, + so ``Location_MaxIntensity_*`` is bit-identical to the numpy backend. (We do + NOT use the numba kernel's ``>=``-last argmax here, by design — exact parity + with the existing output is the goal.) + """ + n = int(labels_sorted.size) + max_x = numpy.zeros(n) + max_y = numpy.zeros(n) + max_z = numpy.zeros(n) + for k in range(n): + sel = seg0 == k + vals_k = values[sel] + if vals_k.size == 0: + continue + # Constant labels + single index reproduces the reference's per-object + # labeled call; the position returned indexes into vals_k's raster order. + lbl = int(labels_sorted[k]) + labels_k = numpy.full(vals_k.size, lbl, dtype=numpy.int64) + pos = numpy.array( + scipy.ndimage.maximum_position(vals_k, labels_k, numpy.array([lbl])), + dtype=int, + ).ravel()[0] + max_x[k] = xc[sel][pos] + max_y[k] = yc[sel][pos] + max_z[k] = zc[sel][pos] + return max_x, max_y, max_z diff --git a/test/test_backend_correctness.py b/test/test_backend_correctness.py new file mode 100644 index 0000000..55a8279 --- /dev/null +++ b/test/test_backend_correctness.py @@ -0,0 +1,90 @@ +"""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 +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" + ) + # 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" + + +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) From 8673cf45a62a931cf5387cdff84210f871e64ffb Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 05:01:34 +0200 Subject: [PATCH 02/15] refactor(numba): compute max_position inside the numba kernel MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Move Location_MaxIntensity_* out of the host scipy per-object call and into the fused segment_moments kernel via a deterministic `>=`-last argmax (records the max pixel's coordinates in the same single pass). scipy.ndimage.maximum_position's labeled tie-break is `argsort` (quicksort) + last-write-wins, i.e. an arbitrary tied pixel that is not stable across numpy versions — so there is no stable rule to replicate. On real continuous data the max is unique, so the kernel's `>=`-last result is bit-identical to scipy (the correctness harness confirms 2D/3D, edge on/off); only exact-value ties can differ, and the kernel's rule is the more reproducible of the two. Drops the now-unused max_position_per_object host helper (and its scipy import) from the primitive layer. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../core/numba/measureobjectintensity.py | 35 +++++++++------ src/cp_measure/primitives/__init__.py | 2 - src/cp_measure/primitives/_segment_numba.py | 27 ++++++++---- src/cp_measure/primitives/segment.py | 44 ------------------- 4 files changed, 41 insertions(+), 67 deletions(-) diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index 64580f1..399c24e 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -2,10 +2,12 @@ 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 run as fused single-pass numba kernels over a flat segment -representation (:mod:`cp_measure.primitives`); ``Location_MaxIntensity_*`` is -computed on the host via ``scipy.ndimage.maximum_position`` to reproduce the -numpy backend's tie-break exactly. +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 @@ -40,7 +42,6 @@ from cp_measure.primitives import ( flatten_labeled, label_to_idx_lut, - max_position_per_object, ) from cp_measure.primitives._segment_numba import ( segment_moments, @@ -65,7 +66,7 @@ def get_intensity( elif pixels.ndim == 3 and masks.ndim == 2: # 3D image, 2D mask masks = masks.reshape(1, *masks.shape) - lut, labels, nobjects = label_to_idx_lut(masks) + lut, _labels, nobjects = label_to_idx_lut(masks) integrated_intensity = numpy.zeros(nobjects) mean_intensity = numpy.zeros(nobjects) @@ -88,9 +89,21 @@ def get_intensity( has_objects = values.size > 0 if has_objects: - (count, sumI, minI, maxI, sx, sy, sz, sxI, syI, szI) = segment_moments( - values, seg0, xc, yc, zc, nobjects - ) + ( + 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 @@ -117,10 +130,6 @@ def get_intensity( mad_intensity, ) = segment_quantiles(values, seg0, count, nobjects, 1.0 / orig_ndim) - max_x, max_y, max_z = max_position_per_object( - values, seg0, xc, yc, zc, labels - ) - if edge_measurements: integrated_intensity_edge = numpy.zeros(nobjects) mean_intensity_edge = numpy.zeros(nobjects) diff --git a/src/cp_measure/primitives/__init__.py b/src/cp_measure/primitives/__init__.py index 673b714..b913779 100644 --- a/src/cp_measure/primitives/__init__.py +++ b/src/cp_measure/primitives/__init__.py @@ -12,11 +12,9 @@ from cp_measure.primitives.segment import ( flatten_labeled, label_to_idx_lut, - max_position_per_object, ) __all__ = [ "flatten_labeled", "label_to_idx_lut", - "max_position_per_object", ] diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py index 4c236a3..45669d0 100644 --- a/src/cp_measure/primitives/_segment_numba.py +++ b/src/cp_measure/primitives/_segment_numba.py @@ -17,14 +17,22 @@ def segment_moments(values, seg0, xc, yc, zc, n): """One pass over the flat pixels accumulating, per segment: - count, sum, min, max, and the six centroid cross-sums - (Sx, Sy, Sz, Sx*I, Sy*I, Sz*I). max position is intentionally NOT computed - here (it is done on the host via scipy for exact parity). + 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) @@ -35,22 +43,25 @@ def segment_moments(values, seg0, xc, yc, zc, n): 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]: + if v >= maxI[k]: # >= -> keep LAST max in raster order maxI[k] = v - x = xc[i] - y = yc[i] - z = zc[i] + 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, sx, sy, sz, sxI, syI, szI + return count, sumI, minI, maxI, mx, my, mz, sx, sy, sz, sxI, syI, szI @njit(cache=True) diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py index c7dcdf8..4874ee0 100644 --- a/src/cp_measure/primitives/segment.py +++ b/src/cp_measure/primitives/segment.py @@ -13,7 +13,6 @@ """ import numpy -import scipy.ndimage from numpy.typing import NDArray @@ -64,46 +63,3 @@ def flatten_labeled( yc = mesh_y[lmask].astype(numpy.float64) zc = mesh_z[lmask].astype(numpy.float64) return values, seg0, xc, yc, zc - - -def max_position_per_object( - values: NDArray[numpy.float64], - seg0: NDArray[numpy.int64], - xc: NDArray[numpy.float64], - yc: NDArray[numpy.float64], - zc: NDArray[numpy.float64], - labels_sorted: NDArray[numpy.integer], -) -> tuple[ - NDArray[numpy.float64], NDArray[numpy.float64], NDArray[numpy.float64] -]: - """Per-object argmax position via ``scipy.ndimage.maximum_position``. - - Reproduces the shipped numpy implementation EXACTLY, including scipy's - labeled tie-break on exact-equal maxima (an implementation artifact that is - neither reliably first nor last). The reference calls ``maximum_position`` - once per object on that object's pixels in raster order; this does the same, - so ``Location_MaxIntensity_*`` is bit-identical to the numpy backend. (We do - NOT use the numba kernel's ``>=``-last argmax here, by design — exact parity - with the existing output is the goal.) - """ - n = int(labels_sorted.size) - max_x = numpy.zeros(n) - max_y = numpy.zeros(n) - max_z = numpy.zeros(n) - for k in range(n): - sel = seg0 == k - vals_k = values[sel] - if vals_k.size == 0: - continue - # Constant labels + single index reproduces the reference's per-object - # labeled call; the position returned indexes into vals_k's raster order. - lbl = int(labels_sorted[k]) - labels_k = numpy.full(vals_k.size, lbl, dtype=numpy.int64) - pos = numpy.array( - scipy.ndimage.maximum_position(vals_k, labels_k, numpy.array([lbl])), - dtype=int, - ).ravel()[0] - max_x[k] = xc[sel][pos] - max_y[k] = yc[sel][pos] - max_z[k] = zc[sel][pos] - return max_x, max_y, max_z From 2c23bc5eeec04f6f87b4da14e075d8babb36ba67 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 05:21:39 +0200 Subject: [PATCH 03/15] refactor(numba): tidy primitive layer after review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - flatten_labeled: derive (z,y,x) coords from numpy.nonzero(lmask) instead of materialising three full-volume mgrid arrays then masking them — same coords in the same C order, no per-call O(volume) temporaries. - label_to_idx_lut: drop the unused sorted-labels return value (now just (lut, n)); the max_position-in-kernel refactor removed its only consumer. - add a lighter segment_stats kernel (count/sum/min/max) and use it for the edge path, replacing the segment_moments call that needed throwaway zero coordinate arrays and discarded the centroid cross-sums. No behaviour change; correctness harness + full suite stay green. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../core/numba/measureobjectintensity.py | 8 +++--- src/cp_measure/primitives/_segment_numba.py | 24 ++++++++++++++++++ src/cp_measure/primitives/segment.py | 25 ++++++++----------- 3 files changed, 37 insertions(+), 20 deletions(-) diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index 399c24e..2f4abe3 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -47,6 +47,7 @@ segment_moments, segment_quantiles, segment_resid_sumsq, + segment_stats, ) @@ -66,7 +67,7 @@ def get_intensity( elif pixels.ndim == 3 and masks.ndim == 2: # 3D image, 2D mask masks = masks.reshape(1, *masks.shape) - lut, _labels, nobjects = label_to_idx_lut(masks) + lut, nobjects = label_to_idx_lut(masks) integrated_intensity = numpy.zeros(nobjects) mean_intensity = numpy.zeros(nobjects) @@ -143,10 +144,7 @@ def get_intensity( e_seg0 = lut[masks[emask]] if e_values.size > 0: - zeros = numpy.zeros(e_values.size) - (ecount, esum, emin, emax, *_rest) = segment_moments( - e_values, e_seg0, zeros, zeros, zeros, nobjects - ) + 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"): diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py index 45669d0..0a629a4 100644 --- a/src/cp_measure/primitives/_segment_numba.py +++ b/src/cp_measure/primitives/_segment_numba.py @@ -64,6 +64,30 @@ def segment_moments(values, seg0, xc, yc, zc, n): 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 segment_resid_sumsq(values, seg0, n, mean): """Second pass: per-segment sum of squared residuals (for population std).""" diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py index 4874ee0..1dadc8a 100644 --- a/src/cp_measure/primitives/segment.py +++ b/src/cp_measure/primitives/segment.py @@ -18,14 +18,13 @@ def label_to_idx_lut( masks: NDArray[numpy.integer], -) -> tuple[NDArray[numpy.int64], NDArray[numpy.integer], int]: +) -> tuple[NDArray[numpy.int64], int]: """Build a ``label -> 0..n-1`` lookup over the sorted positive labels. - Returns ``(lut, labels_sorted, n)`` where ``lut[label]`` is the segment - index (and ``-1`` for absent labels / background), ``labels_sorted`` are the - ascending positive labels, 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. + 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. """ unique = numpy.unique(masks) labels = unique[unique > 0] @@ -33,7 +32,7 @@ def label_to_idx_lut( 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, labels, n + return lut, n def flatten_labeled( @@ -50,16 +49,12 @@ def flatten_labeled( """Flatten a labeled (Z, Y, X) image to ``(values, seg0, xc, yc, zc)``. Non-finite pixels and background are dropped, matching the numpy reference's - ``(masks > 0) & isfinite(pixels)`` mask. Coordinates use ``numpy.mgrid`` over - the full volume exactly as the reference does, so positions line up. + ``(masks > 0) & isfinite(pixels)`` mask. ``numpy.nonzero`` yields the (z, y, x) + coordinates of the kept pixels directly, in the same C order as ``pixels[lmask]`` + — no full-volume coordinate grids are materialised. """ lmask = (masks > 0) & numpy.isfinite(pixels) values = pixels[lmask].astype(numpy.float64) seg0 = lut[masks[lmask]] - mesh_z, mesh_y, mesh_x = numpy.mgrid[ - 0 : masks.shape[0], 0 : masks.shape[1], 0 : masks.shape[2] - ] - xc = mesh_x[lmask].astype(numpy.float64) - yc = mesh_y[lmask].astype(numpy.float64) - zc = mesh_z[lmask].astype(numpy.float64) + zc, yc, xc = (c.astype(numpy.float64) for c in numpy.nonzero(lmask)) return values, seg0, xc, yc, zc From 15f4b411a21f94be20ebf72a74000b39fac76b58 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 15:36:36 +0200 Subject: [PATCH 04/15] perf(numba): flatten the labeled image in a numba scan, not numpy MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit flatten_labeled built the flat (values, seg0, coords) arrays with a numpy (masks>0)&isfinite mask + numpy.nonzero + two fancy-index gathers — several full-image passes plus a boolean-array allocation, and the dominant cost of the non-edge path. Replace it with flatten_numba: two grid scans (count, then fill) in a single @njit kernel, coordinates taken from the loop indices. The flat-segment kernels and the rest of the backend are unchanged — only how the flat arrays are built. Measured (single image, non-edge core): flatten step ~4-10x faster (10x at 1024^2), full core ~1.1x (256^2) / ~1.5x (1024^2); the gain grows with image size. Bit-identical output (correctness harness stays green). The numpy flatten_labeled (its only consumer) is removed; primitives/segment.py now holds just the numpy label->index lookup, the numba layer owns the flatten. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../core/numba/measureobjectintensity.py | 12 +++--- src/cp_measure/primitives/__init__.py | 16 ++++--- src/cp_measure/primitives/_segment_numba.py | 40 ++++++++++++++++++ src/cp_measure/primitives/segment.py | 42 +++---------------- 4 files changed, 60 insertions(+), 50 deletions(-) diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index 2f4abe3..59926bc 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -39,11 +39,9 @@ STD_INTENSITY_EDGE, UPPER_QUARTILE_INTENSITY, ) -from cp_measure.primitives import ( - flatten_labeled, - label_to_idx_lut, -) +from cp_measure.primitives import label_to_idx_lut from cp_measure.primitives._segment_numba import ( + flatten_numba, segment_moments, segment_quantiles, segment_resid_sumsq, @@ -86,7 +84,11 @@ def get_intensity( max_y = numpy.zeros(nobjects) max_z = numpy.zeros(nobjects) - values, seg0, xc, yc, zc = flatten_labeled(masks, masked_image, lut) + values, seg0, xc, yc, zc = flatten_numba( + numpy.ascontiguousarray(masks), + numpy.ascontiguousarray(masked_image, dtype=numpy.float64), + lut, + ) has_objects = values.size > 0 if has_objects: diff --git a/src/cp_measure/primitives/__init__.py b/src/cp_measure/primitives/__init__.py index b913779..1f67af8 100644 --- a/src/cp_measure/primitives/__init__.py +++ b/src/cp_measure/primitives/__init__.py @@ -3,18 +3,16 @@ Backend-agnostic building blocks that the per-backend feature implementations (``cp_measure.core`` = numpy, ``cp_measure.core.numba`` = numba, ...) compose. -The host helpers here flatten a labeled image into a 1-D *segment* representation -(values + 0-based segment index + per-axis coordinates). 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. +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`` here builds the label→index lookup; the +flattening + reductions themselves live in ``_segment_numba``. """ -from cp_measure.primitives.segment import ( - flatten_labeled, - label_to_idx_lut, -) +from cp_measure.primitives.segment import label_to_idx_lut __all__ = [ - "flatten_labeled", "label_to_idx_lut", ] diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py index 0a629a4..18453da 100644 --- a/src/cp_measure/primitives/_segment_numba.py +++ b/src/cp_measure/primitives/_segment_numba.py @@ -13,6 +13,46 @@ 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 and ``pixels`` float64. + """ + 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: diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py index 1dadc8a..c671234 100644 --- a/src/cp_measure/primitives/segment.py +++ b/src/cp_measure/primitives/segment.py @@ -1,15 +1,10 @@ -"""Host-side segment helpers (numpy/scipy, backend-agnostic). +"""Host-side segment helpers (numpy, backend-agnostic). -A labeled image is reduced to a flat *segment* form: - -* ``values[M]`` finite labeled pixel intensities (float64) -* ``seg0[M]`` 0-based segment index = rank of the pixel's label among the - sorted positive labels (``label_to_idx``) -* ``xc/yc/zc[M]``per-pixel coordinates (z = 0 plane for 2D input) - -Segment kernels (numpy or numba) then loop over these flat M-length arrays with -no notion of image shape or batch axis. 2D vs 3D differ only in which coordinate -arrays are populated; a batch differs only in how ``seg0`` offsets are assigned. +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 @@ -33,28 +28,3 @@ def label_to_idx_lut( lut = numpy.full(max_label + 1, -1, dtype=numpy.int64) lut[labels] = numpy.arange(n, dtype=numpy.int64) return lut, n - - -def flatten_labeled( - masks: NDArray[numpy.integer], - pixels: NDArray[numpy.floating], - lut: NDArray[numpy.int64], -) -> tuple[ - NDArray[numpy.float64], - NDArray[numpy.int64], - NDArray[numpy.float64], - NDArray[numpy.float64], - NDArray[numpy.float64], -]: - """Flatten a labeled (Z, Y, X) image to ``(values, seg0, xc, yc, zc)``. - - Non-finite pixels and background are dropped, matching the numpy reference's - ``(masks > 0) & isfinite(pixels)`` mask. ``numpy.nonzero`` yields the (z, y, x) - coordinates of the kept pixels directly, in the same C order as ``pixels[lmask]`` - — no full-volume coordinate grids are materialised. - """ - lmask = (masks > 0) & numpy.isfinite(pixels) - values = pixels[lmask].astype(numpy.float64) - seg0 = lut[masks[lmask]] - zc, yc, xc = (c.astype(numpy.float64) for c in numpy.nonzero(lmask)) - return values, seg0, xc, yc, zc From 66953d611cab029797d49487297e2063dc20637c Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 15:48:09 +0200 Subject: [PATCH 05/15] refactor(numba): trim _detect to HAS_NUMBA; avoid float64 image copy - _detect.py: drop the unused HAS_JAX / HAS_JAX_GPU flags. Besides being dead for this PR, HAS_JAX_GPU eagerly imported jax at module load whenever jax was installed, just to set a flag nothing reads. jax detection lands with the jax backend; HAS_NUMBA alone establishes the find_spec pattern. - flatten the image without a forced float64 copy: pass masked_image through ascontiguousarray without dtype=, and let flatten_numba upcast the kept values into its float64 output. Avoids a full-image float64 temporary for non-float64 inputs (e.g. float32 microscopy data); bit-identical for float64. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/_detect.py | 16 +++------------- .../core/numba/measureobjectintensity.py | 2 +- src/cp_measure/primitives/_segment_numba.py | 3 ++- 3 files changed, 6 insertions(+), 15 deletions(-) diff --git a/src/cp_measure/_detect.py b/src/cp_measure/_detect.py index d9ca785..9dd8f9f 100644 --- a/src/cp_measure/_detect.py +++ b/src/cp_measure/_detect.py @@ -2,23 +2,13 @@ Detected ONCE at import via ``importlib.util.find_spec`` — availability is checked without importing the package or catching ImportErrors. Dispatch reads -these flags; the resolved backend path is then called directly and unguarded +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 -HAS_JAX: bool = importlib.util.find_spec("jax") is not None - - -def _detect_jax_gpu() -> bool: - if not HAS_JAX: - return False - import jax - - return jax.default_backend() != "cpu" - - -HAS_JAX_GPU: bool = _detect_jax_gpu() diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index 59926bc..47044f9 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -86,7 +86,7 @@ def get_intensity( values, seg0, xc, yc, zc = flatten_numba( numpy.ascontiguousarray(masks), - numpy.ascontiguousarray(masked_image, dtype=numpy.float64), + numpy.ascontiguousarray(masked_image), lut, ) has_objects = values.size > 0 diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py index 18453da..9bcdef4 100644 --- a/src/cp_measure/primitives/_segment_numba.py +++ b/src/cp_measure/primitives/_segment_numba.py @@ -20,7 +20,8 @@ def flatten_numba(masks, pixels, lut): 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 and ``pixels`` float64. + ``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 From 2ed44ceff746eee29fcd12a7412b7dca5f0b3f77 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 16:02:35 +0200 Subject: [PATCH 06/15] perf(numba): enumerate labels via find_objects, not numpy.unique MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit label_to_idx_lut used numpy.unique(masks) — a full-image sort — to find the present labels. scipy.ndimage.find_objects (scipy is already a core dep) returns the same ascending present-label set in one O(P) pass, giving a bit-identical LUT ~3-5x faster (12.4->3.5 ms at 1024^2; 21.9->4.4 ms on a 32x240x240 volume). Trick borrowed from Alan's pure-numpy speedup (afermg/cp_measure#55); unlike its percentile/MAD changes, this one preserves output exactly (verified identical). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/primitives/segment.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py index c671234..e5b5894 100644 --- a/src/cp_measure/primitives/segment.py +++ b/src/cp_measure/primitives/segment.py @@ -8,6 +8,7 @@ """ import numpy +import scipy.ndimage from numpy.typing import NDArray @@ -20,9 +21,17 @@ def label_to_idx_lut( 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. """ - unique = numpy.unique(masks) - labels = unique[unique > 0] + 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) From e277c5e2ab3237b21ed66f9b1735b74e2543a9b6 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Tue, 2 Jun 2026 16:29:02 +0200 Subject: [PATCH 07/15] perf(numba): compute edge inner-boundary in a kernel, not skimage (2D) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Profiling the sparse-large regime (1024^2, 64 obj, edge on) showed skimage.find_boundaries was ~37% of the call (~20-29 ms) — the morphology dominates, not the scan. A one-pass numba inner-boundary kernel (4-neighbour check, the cp_measure_fast approach) is bit-identical to find_boundaries( mode="inner") and 12-27x faster, verified exact across (H,W) and (1,H,W). Used for 2D planes (Z==1); true 3D keeps skimage (6-neighbourhood). Single-image 1024^2/64 edge-on drops ~47->32 ms, and per-image batch ~445->264 ms. No correctness change (exact boundary match; harness stays green). Co-Authored-By: Claude Opus 4.8 (1M context) --- .../core/numba/measureobjectintensity.py | 9 +++++-- src/cp_measure/primitives/_segment_numba.py | 26 +++++++++++++++++++ 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index 47044f9..d04e8c6 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -42,6 +42,7 @@ from cp_measure.primitives import label_to_idx_lut from cp_measure.primitives._segment_numba import ( flatten_numba, + inner_boundary, segment_moments, segment_quantiles, segment_resid_sumsq, @@ -140,8 +141,12 @@ def get_intensity( min_intensity_edge = numpy.zeros(nobjects) max_intensity_edge = numpy.zeros(nobjects) - outlines = skimage.segmentation.find_boundaries(masks, mode="inner") - emask = outlines > 0 + # 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]] diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py index 9bcdef4..ffb9148 100644 --- a/src/cp_measure/primitives/_segment_numba.py +++ b/src/cp_measure/primitives/_segment_numba.py @@ -129,6 +129,32 @@ def segment_stats(values, seg0, n): 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).""" From eefb4d442043c242c80be9e7fb3d4a91489b50eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Al=C3=A1n=20F=2E=20Mu=C3=B1oz?= Date: Tue, 2 Jun 2026 21:29:21 -0400 Subject: [PATCH 08/15] lock: add numba to lock --- uv.lock | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) 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" From 69c7e48d5e8e1374ff202c5ff51930470b5540aa Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 3 Jun 2026 15:45:39 +0200 Subject: [PATCH 09/15] refactor(numba): reword install hint; inline primitives import Address PR #54 review: - bulk._dispatch: reword the absent-numba RuntimeError from the imperative "install it via" to "you can install it via" (avoid issuing pip commands imperatively at the user). - primitives is an internal layer with no public API to curate; import label_to_idx_lut directly from primitives.segment (matching how the _segment_numba kernels are already imported) and drop the __init__ re-export. Documents the convention in the package docstring. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/bulk.py | 2 +- .../core/numba/measureobjectintensity.py | 2 +- src/cp_measure/primitives/__init__.py | 15 +++++++-------- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py index 1824431..8157ad6 100644 --- a/src/cp_measure/bulk.py +++ b/src/cp_measure/bulk.py @@ -77,7 +77,7 @@ def _dispatch(name: str) -> dict[str, Callable]: if not HAS_NUMBA: raise RuntimeError( "accelerator 'numba' selected but numba is not installed; " - "install it via `pip install cp_measure[numba]`" + "you can install it via `pip install cp_measure[numba]`" ) return _numba_registries()[name] if _ACCELERATOR == "fastest": diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py index d04e8c6..bd058c5 100644 --- a/src/cp_measure/core/numba/measureobjectintensity.py +++ b/src/cp_measure/core/numba/measureobjectintensity.py @@ -39,7 +39,7 @@ STD_INTENSITY_EDGE, UPPER_QUARTILE_INTENSITY, ) -from cp_measure.primitives import label_to_idx_lut +from cp_measure.primitives.segment import label_to_idx_lut from cp_measure.primitives._segment_numba import ( flatten_numba, inner_boundary, diff --git a/src/cp_measure/primitives/__init__.py b/src/cp_measure/primitives/__init__.py index 1f67af8..3619525 100644 --- a/src/cp_measure/primitives/__init__.py +++ b/src/cp_measure/primitives/__init__.py @@ -7,12 +7,11 @@ 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`` here builds the label→index lookup; the -flattening + reductions themselves live in ``_segment_numba``. -""" - -from cp_measure.primitives.segment import label_to_idx_lut +a rewrite. The numpy ``label_to_idx_lut`` (in ``segment``) builds the +label→index lookup; the flattening + reductions themselves live in +``_segment_numba``. -__all__ = [ - "label_to_idx_lut", -] +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. +""" From 4ed7d7912187bed4fe1a544d75cfc389a9ac88f2 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 3 Jun 2026 20:52:11 +0200 Subject: [PATCH 10/15] feat(primitives): canonical (B,Z,Y,X) batch-shape helper (to_bzyx) Shared foundational helper used by the numba intensity/granularity/zernike backends to normalise any input (2D/3D/4D/list) to the canonical batch-of-volumes form: single image = batch of 1, returning a dict for a lone image/volume and a list of dicts for a batch. Pure numpy, no numba. Extracted to its own PR so it can be reviewed first and unblock the feature backends (#56/#57/#58). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/primitives/shapes.py | 83 +++++++++++++++++++++++++++++ test/test_primitives_shapes.py | 80 +++++++++++++++++++++++++++ 2 files changed, 163 insertions(+) create mode 100644 src/cp_measure/primitives/shapes.py create mode 100644 test/test_primitives_shapes.py 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_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) From c261573ca7851a9bcdc24018fee19ea0a06d1b41 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 3 Jun 2026 17:13:15 +0200 Subject: [PATCH 11/15] feat(accelerator): numba granularity backend + (B,Z,Y,X) batch shape MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Second accelerator on the #49/#54 dispatch seam: set_accelerator("numba") now routes `granularity` (2D) to a numba backend alongside `intensity`. Kernels (core/numba/_granularity.py, single-threaded @njit, bit-exact vs skimage): - disk erosion/dilation via row-decomposed van-Herk/Gil-Werman sliding min/max (O(r·HW)) for the background opening - 5-tap disk(1) erosion/dilation for the spectrum loop - reconstruction by dilation via the Vincent (1993) FIFO-hybrid (1 fwd + 1 bwd raster, then a ring-buffer FIFO with a per-pixel in-queue flag — O(N), no overflow), replacing the per-iteration scipy/skimage morphology Wrapper (core/numba/measuregranularity.py): mirrors the numpy baseline but keeps resampling on scipy map_coordinates, swaps the morphology for the kernels, reads per-object means back via a sparse point-query of the reconstructed image (no full-res upsample) + dense bincount, and applies the cascaded reconstruction mask (rec_g <= rec_{g-1} <= pixels, exact by monotonicity). 3D / non-2D falls back to the numpy baseline. Batch shape: new primitives/shapes.py `to_bzyx` normalises any input to the canonical (B,Z,Y,X) form (2D/3D/4D array or ragged list), single image = batch of 1, returning a dict for a lone image/volume and a list of dicts for a batch. A 3D array stays one volume (never a batch), preserving existing semantics. Correctness: morphology bit-exact vs skimage (radii 1..10 + cascade chain); wrapper matches the numpy baseline within tol for fullres + all subsample modes, single/1-px/empty/label-gap objects, batch (4D + ragged), and 3D fallback. Full suite 177 passed. Benchmark (1080², 144 objects): fullres 9000->1369 ms (6.57x), default 1178->94 ms (12.47x). Follow-up: migrate the intensity backend onto to_bzyx (golden-test guarded). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/bulk.py | 13 +- src/cp_measure/core/numba/__init__.py | 7 +- src/cp_measure/core/numba/_granularity.py | 284 ++++++++++++++++++ .../core/numba/measuregranularity.py | 159 ++++++++++ src/cp_measure/primitives/shapes.py | 5 - test/test_backend_correctness.py | 21 ++ test/test_granularity_backend.py | 114 +++++++ test/test_granularity_kernels.py | 102 +++++++ 8 files changed, 694 insertions(+), 11 deletions(-) create mode 100644 src/cp_measure/core/numba/_granularity.py create mode 100644 src/cp_measure/core/numba/measuregranularity.py create mode 100644 test/test_granularity_backend.py create mode 100644 test/test_granularity_kernels.py diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py index 8157ad6..f824940 100644 --- a/src/cp_measure/bulk.py +++ b/src/cp_measure/bulk.py @@ -48,16 +48,23 @@ def _numba_registries() -> dict[str, dict[str, Callable]]: """Registries for the 'numba' accelerator. - Composes the numba implementations (currently ``intensity`` only) with the + 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_intensity as _numba_intensity + from cp_measure.core.numba import ( + get_granularity as _numba_granularity, + get_intensity as _numba_intensity, + ) return { - "core": {**_CORE, "intensity": _numba_intensity}, + "core": { + **_CORE, + "intensity": _numba_intensity, + "granularity": _numba_granularity, + }, "correlation": _CORRELATION, } diff --git a/src/cp_measure/core/numba/__init__.py b/src/cp_measure/core/numba/__init__.py index 671d723..65b9beb 100644 --- a/src/cp_measure/core/numba/__init__.py +++ b/src/cp_measure/core/numba/__init__.py @@ -4,11 +4,12 @@ or globally via ``cp_measure.set_accelerator("numba")``. Requires the optional ``numba`` extra; availability is gated by ``cp_measure._detect.HAS_NUMBA``. -This backend currently accelerates ``intensity`` only; the global "numba" -accelerator composes it with the numpy implementations of every other feature +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_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..80fd55d --- /dev/null +++ b/src/cp_measure/core/numba/_granularity.py @@ -0,0 +1,284 @@ +"""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 raster-until-convergence, 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 + + +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.empty((H, W), np.float64) + fill = -np.inf if is_max else np.inf + for i in range(H): + for j in range(W): + out[i, j] = fill + + 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 + for i in range(H): + pr = i + radius + dy + # forward prefix-min/max within blocks of size k + for x in range(L): + if x % k == 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 + # backward suffix-min/max within blocks of size k + for x in range(L - 1, -1, -1): + if x == L - 1 or x % k == 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 + # 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 reconstruction_by_dilation_2d(seed, mask): + """Morphological reconstruction by dilation (4-connectivity), seed under mask. + + Vincent (1993) hybrid: one forward (TL→BR) + one backward (BR→TL) raster + geodesic dilation under ``mask``, the 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)``. + + 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 + + # forward raster: propagate from up/left causal neighbours + 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 + + cap = H * W + 1 + queue = np.empty(cap, np.int64) + inq = np.zeros((H, W), np.uint8) + head = 0 + tail = 0 + + # 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] = i * W + 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 // W + j = code % W + 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] = (i - 1) * W + 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] = (i + 1) * W + 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] = i * W + (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] = i * W + (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..5fce5a7 --- /dev/null +++ b/src/cp_measure/core/numba/measuregranularity.py @@ -0,0 +1,159 @@ +"""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). Resampling stays on ``scipy.ndimage.map_coordinates``, as in every +backend. + +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 ( + 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)) + + 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) + if ero.max() == 0: + rec = ero + else: + rec = reconstruction_by_dilation_2d(ero, recon_mask) + recon_mask = rec # cascade: rec_g <= rec_{g-1} <= pixels (exact) + + if needs_resize: + rec_valid = scipy.ndimage.map_coordinates(rec, (sy, sx), order=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/primitives/shapes.py b/src/cp_measure/primitives/shapes.py index 0c94059..b5942f7 100644 --- a/src/cp_measure/primitives/shapes.py +++ b/src/cp_measure/primitives/shapes.py @@ -20,11 +20,6 @@ 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 diff --git a/test/test_backend_correctness.py b/test/test_backend_correctness.py index 55a8279..52ab275 100644 --- a/test/test_backend_correctness.py +++ b/test/test_backend_correctness.py @@ -61,6 +61,23 @@ def test_numba_intensity_matches_numpy(dim, 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") @@ -69,6 +86,9 @@ def test_set_accelerator_numba_composes_with_numpy(): 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" @@ -77,6 +97,7 @@ def test_set_accelerator_numba_composes_with_numpy(): 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): diff --git a/test/test_granularity_backend.py b/test/test_granularity_backend.py new file mode 100644 index 0000000..c1948f6 --- /dev/null +++ b/test/test_granularity_backend.py @@ -0,0 +1,114 @@ +"""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..1e9d148 --- /dev/null +++ b/test/test_granularity_kernels.py @@ -0,0 +1,102 @@ +"""Bit-exact validation of the numba granularity morphology kernels vs skimage.""" + +import numpy +import pytest +import skimage.morphology as M + +from cp_measure.core.numba._granularity import ( + _disk_halfwidths, + disk_dilation_2d, + disk_erosion_2d, + dilation_4conn_2d, + erosion_4conn_2d, + reconstruction_by_dilation_2d, +) + + +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 From 90b5131ea1f302057186cc90f0aee69d8c492281 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 3 Jun 2026 17:24:06 +0200 Subject: [PATCH 12/15] refactor(numba): tidy granularity kernels after review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Quality pass (no behaviour change; bit-exact, 107 tests pass): - _disk_reduce: np.full for the out init (was a hand-rolled nested loop); replace the `x % k` modulo in the per-pixel VHG prefix/suffix loops with a running counter (one hoisted modulo per row-band) - get_granularity: drop the `ero.max() == 0` early-out — it scanned the whole image every one of the 16 iterations, while reconstruction already returns zeros cheaply on an all-zero seed - shapes.to_bzyx: document that it is a pure batch normaliser (no mask-broadcast; per-element ndim dispatch stays in each backend) Benchmark (1080², 144 obj): fullres 6.57x->6.86x, default 12.47x->12.58x. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_granularity.py | 23 +++++++++++-------- .../core/numba/measuregranularity.py | 9 ++++---- src/cp_measure/primitives/shapes.py | 5 ++++ 3 files changed, 23 insertions(+), 14 deletions(-) diff --git a/src/cp_measure/core/numba/_granularity.py b/src/cp_measure/core/numba/_granularity.py index 80fd55d..fc34c57 100644 --- a/src/cp_measure/core/numba/_granularity.py +++ b/src/cp_measure/core/numba/_granularity.py @@ -44,22 +44,20 @@ def _disk_reduce(P, radius, hx, H, W, is_max): rows. ``is_max`` selects dilation (max) vs erosion (min). """ L = W + 2 * radius - out = np.empty((H, W), np.float64) - fill = -np.inf if is_max else np.inf - for i in range(H): - for j in range(W): - out[i, j] = fill + 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 + # forward prefix-min/max within blocks of size k (counter tracks x % k) + c = 0 for x in range(L): - if x % k == 0: + if c == 0: g[x] = P[pr, x] else: a = g[x - 1] @@ -68,9 +66,13 @@ def _disk_reduce(P, radius, hx, H, W, is_max): g[x] = a if a > b else b else: g[x] = a if a < b else b - # backward suffix-min/max within blocks of size k + 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 x % k == k - 1: + if x == L - 1 or c == k - 1: h[x] = P[pr, x] else: a = h[x + 1] @@ -79,6 +81,9 @@ def _disk_reduce(P, radius, hx, H, W, 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 diff --git a/src/cp_measure/core/numba/measuregranularity.py b/src/cp_measure/core/numba/measuregranularity.py index 5fce5a7..18d9701 100644 --- a/src/cp_measure/core/numba/measuregranularity.py +++ b/src/cp_measure/core/numba/measuregranularity.py @@ -141,11 +141,10 @@ def _label_mean(values_at_in_obj): recon_mask = pixels for granularity_id in range(1, ng + 1): ero = erosion_4conn_2d(ero) - if ero.max() == 0: - rec = ero - else: - rec = reconstruction_by_dilation_2d(ero, recon_mask) - recon_mask = rec # cascade: rec_g <= rec_{g-1} <= pixels (exact) + # 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 = scipy.ndimage.map_coordinates(rec, (sy, sx), order=1) diff --git a/src/cp_measure/primitives/shapes.py b/src/cp_measure/primitives/shapes.py index b5942f7..0c94059 100644 --- a/src/cp_measure/primitives/shapes.py +++ b/src/cp_measure/primitives/shapes.py @@ -20,6 +20,11 @@ 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 78788a4914425104f887285c6a3f562f5343b272 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 3 Jun 2026 17:28:46 +0200 Subject: [PATCH 13/15] style: ruff format the granularity backend + tests CI lint (nix fmt -> ruff-format) reformats lines over the 88-col default; apply ruff format to match. Also fix a stale module docstring in _granularity.py (reconstruction is the Vincent FIFO-hybrid, not raster-until-convergence). No behaviour change; 99 granularity tests pass. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_granularity.py | 20 ++++++---- .../core/numba/measuregranularity.py | 40 +++++++++++++++---- test/test_granularity_backend.py | 4 +- test/test_granularity_kernels.py | 20 +++++++--- 4 files changed, 63 insertions(+), 21 deletions(-) diff --git a/src/cp_measure/core/numba/_granularity.py b/src/cp_measure/core/numba/_granularity.py index fc34c57..25b47c0 100644 --- a/src/cp_measure/core/numba/_granularity.py +++ b/src/cp_measure/core/numba/_granularity.py @@ -11,7 +11,7 @@ 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 raster-until-convergence, matching + (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 @@ -51,7 +51,9 @@ def _disk_reduce(P, radius, hx, H, W, is_max): 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) + 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) @@ -106,7 +108,9 @@ def disk_erosion_2d(img: NDArray, radius: int) -> NDArray[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) + 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]: @@ -115,7 +119,9 @@ def disk_dilation_2d(img: NDArray, radius: int) -> NDArray[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) + return _disk_reduce( + P, radius, _disk_halfwidths(radius), a.shape[0], a.shape[1], True + ) @njit(cache=True) @@ -221,9 +227,9 @@ def reconstruction_by_dilation_2d(seed, mask): 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] - ) + 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] = i * W + j inq[i, j] = 1 diff --git a/src/cp_measure/core/numba/measuregranularity.py b/src/cp_measure/core/numba/measuregranularity.py index 18d9701..0f3c1a1 100644 --- a/src/cp_measure/core/numba/measuregranularity.py +++ b/src/cp_measure/core/numba/measuregranularity.py @@ -45,7 +45,12 @@ def get_granularity( masks_zyx, pixels_zyx, unwrap = to_bzyx(masks, pixels) results = [ _granularity_one( - m, p, subsample_size, image_sample_size, element_size, granular_spectrum_length + m, + p, + subsample_size, + image_sample_size, + element_size, + granular_spectrum_length, ) for m, p in zip(masks_zyx, pixels_zyx) ] @@ -58,7 +63,12 @@ def _granularity_one( """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 + 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 @@ -81,7 +91,10 @@ def _granularity_2d( # 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 + 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) @@ -89,7 +102,10 @@ def _granularity_2d( # 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 + 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 @@ -119,7 +135,9 @@ def _granularity_2d( 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) + 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: @@ -128,9 +146,13 @@ def _granularity_2d( sx = xx * (float(new_shape[1] - 1) / float(orig_shape[1] - 1)) def _label_mean(values_at_in_obj): - sums = numpy.bincount(labels_in, weights=values_at_in_obj, minlength=max_label + 1)[1:] + 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 + 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) @@ -152,7 +174,9 @@ def _label_mean(values_at_in_obj): rec_valid = rec.ravel()[flat_pos] new_mean = _label_mean(rec_valid) - results[f"Granularity_{granularity_id}"] = (current_mean - new_mean) * 100 / start_mean + results[f"Granularity_{granularity_id}"] = ( + (current_mean - new_mean) * 100 / start_mean + ) current_mean = new_mean return results diff --git a/test/test_granularity_backend.py b/test/test_granularity_backend.py index c1948f6..7c156f8 100644 --- a/test/test_granularity_backend.py +++ b/test/test_granularity_backend.py @@ -30,7 +30,9 @@ def _assert_same(a, b): 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)]) +@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( diff --git a/test/test_granularity_kernels.py b/test/test_granularity_kernels.py index 1e9d148..cb1c2e7 100644 --- a/test/test_granularity_kernels.py +++ b/test/test_granularity_kernels.py @@ -40,7 +40,9 @@ def _spike(H, W): @pytest.mark.parametrize("radius", [1, 2, 3, 5, 7, 10]) -@pytest.mark.parametrize("name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]) +@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)) @@ -48,20 +50,26 @@ def test_disk_erosion_bit_exact(radius, name): @pytest.mark.parametrize("radius", [1, 2, 3, 5, 7, 10]) -@pytest.mark.parametrize("name", ["random", "int_ties", "gradient", "border_hot", "single_spike"]) +@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"]) +@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"]) +@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))) @@ -84,7 +92,9 @@ def test_reconstruction_bit_exact(name): 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)) + got = reconstruction_by_dilation_2d( + seed.astype(numpy.float64), mask.astype(numpy.float64) + ) numpy.testing.assert_array_equal(got, expected) From 89dde43fd6dce9783fc0982c4e7f98a23e96bdf1 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Thu, 4 Jun 2026 07:03:52 +0200 Subject: [PATCH 14/15] perf(numba): numba bilinear gather for granularity resampling (2.88x) The default config (subsample<1) resamples the reconstructed image back to the object sample points via scipy.ndimage.map_coordinates ONCE PER granularity step (16x), re-deriving the floor/frac each time though the sample points never change -- ~67% of default-config runtime. Replace those 16 calls with _granularity.bilinear_gather (order-1, mode=constant cval=0), with floor/frac precomputed once and reused across steps. The three one-off grid resamples (subsample/background) stay on map_coordinates. Matches map_coordinates(order=1) to ~3e-16 per call (FP add-order only; verified vs scipy incl. boundary points in test_bilinear_gather_matches_map_coordinates); the full-output difference vs the old map_coordinates path is ~3e-14, far inside this backend's existing rtol=1e-6 contract (granularity is a within-tolerance lane, not bit-exact). Backend golden test (vs numpy, rtol=1e-6) stays green. Default config 1080^2/144 obj: 741 -> 257 ms (2.88x faster numba backend). Full suite 179 passed, lint clean. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_granularity.py | 29 +++++++++++++++++++ .../core/numba/measuregranularity.py | 20 +++++++++++-- test/test_granularity_kernels.py | 19 ++++++++++++ 3 files changed, 65 insertions(+), 3 deletions(-) diff --git a/src/cp_measure/core/numba/_granularity.py b/src/cp_measure/core/numba/_granularity.py index 25b47c0..e014dcb 100644 --- a/src/cp_measure/core/numba/_granularity.py +++ b/src/cp_measure/core/numba/_granularity.py @@ -24,6 +24,35 @@ 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)``. diff --git a/src/cp_measure/core/numba/measuregranularity.py b/src/cp_measure/core/numba/measuregranularity.py index 0f3c1a1..89dd4d4 100644 --- a/src/cp_measure/core/numba/measuregranularity.py +++ b/src/cp_measure/core/numba/measuregranularity.py @@ -5,8 +5,12 @@ 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). Resampling stays on ``scipy.ndimage.map_coordinates``, as in every -backend. +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 @@ -20,6 +24,7 @@ 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, @@ -144,6 +149,13 @@ def _granularity_2d( 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( @@ -169,7 +181,9 @@ def _label_mean(values_at_in_obj): recon_mask = rec # cascade: rec_g <= rec_{g-1} <= pixels (exact) if needs_resize: - rec_valid = scipy.ndimage.map_coordinates(rec, (sy, sx), order=1) + rec_valid = bilinear_gather( + rec, gy0, gx0, gfy, gfx, rec.shape[0], rec.shape[1] + ) else: rec_valid = rec.ravel()[flat_pos] diff --git a/test/test_granularity_kernels.py b/test/test_granularity_kernels.py index cb1c2e7..1aa15a7 100644 --- a/test/test_granularity_kernels.py +++ b/test/test_granularity_kernels.py @@ -2,10 +2,12 @@ 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, @@ -14,6 +16,23 @@ ) +@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 From c8b9414c661044646d95334efc43d52f3ebb6c06 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Thu, 4 Jun 2026 17:44:23 +0200 Subject: [PATCH 15/15] perf(numba): triple-raster + int32-packed FIFO for granularity reconstruction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Vincent/Robinson hybrid reconstruction: run three forward+backward raster geodesic passes before FIFO seeding (vs one pair) and pack FIFO coordinates as (i<<16)|j int32 decoded by shift/mask (vs i*W+j int64 with int-div/mod). The extra raster pairs cut the FIFO seed count ~50-95%, trading cheap cache-sequential scans for far fewer random-access FIFO updates; the packed coords drop the IDIV/MOD from the hot dequeue and halve the queue's cache footprint. Bit-identical to the previous kernel and to skimage.reconstruction (the result is independent of how many raster passes precede the FIFO) — the overflow-proof N+1 ring buffer + per-pixel in-queue dedup flag are kept. ~1.5x faster reconstruction; end-to-end ~3% (default subsample=0.25, reconstruction is ~20% post-gather) to ~8% (fullres subsample=1.0, reconstruction is ~81%). Non-seeding raster passes factored into _geodesic_raster_fwd/_bwd helpers. +long-geodesic regression test. Full suite 180 green. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_granularity.py | 94 +++++++++++++++-------- test/test_granularity_kernels.py | 20 +++++ 2 files changed, 84 insertions(+), 30 deletions(-) diff --git a/src/cp_measure/core/numba/_granularity.py b/src/cp_measure/core/numba/_granularity.py index e014dcb..27b3b22 100644 --- a/src/cp_measure/core/numba/_granularity.py +++ b/src/cp_measure/core/numba/_granularity.py @@ -200,21 +200,60 @@ def dilation_4conn_2d(img: NDArray) -> NDArray[np.float64]: 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) hybrid: one forward (TL→BR) + one backward (BR→TL) raster - geodesic dilation under ``mask``, the 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 + 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)``. - 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. + 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) @@ -224,27 +263,22 @@ def reconstruction_by_dilation_2d(seed, mask): m = mask[i, j] out[i, j] = s if s < m else m - # forward raster: propagate from up/left causal neighbours - 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 + # 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.int64) + 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 - # backward raster: propagate from down/right; seed the FIFO with pixels whose - # down/right neighbour can still grow from them + # 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] @@ -260,7 +294,7 @@ def reconstruction_by_dilation_2d(seed, mask): 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] = i * W + j + queue[tail] = np.int32((i << 16) | j) inq[i, j] = 1 tail += 1 if tail == cap: @@ -272,8 +306,8 @@ def reconstruction_by_dilation_2d(seed, mask): head += 1 if head == cap: head = 0 - i = code // W - j = code % W + i = code >> np.int32(16) + j = code & mask16 inq[i, j] = 0 v = out[i, j] # up, down, left, right @@ -284,7 +318,7 @@ def reconstruction_by_dilation_2d(seed, mask): out[i - 1, j] = nv if inq[i - 1, j] == 0: inq[i - 1, j] = 1 - queue[tail] = (i - 1) * W + j + queue[tail] = np.int32(((i - 1) << 16) | j) tail += 1 if tail == cap: tail = 0 @@ -295,7 +329,7 @@ def reconstruction_by_dilation_2d(seed, mask): out[i + 1, j] = nv if inq[i + 1, j] == 0: inq[i + 1, j] = 1 - queue[tail] = (i + 1) * W + j + queue[tail] = np.int32(((i + 1) << 16) | j) tail += 1 if tail == cap: tail = 0 @@ -306,7 +340,7 @@ def reconstruction_by_dilation_2d(seed, mask): out[i, j - 1] = nv if inq[i, j - 1] == 0: inq[i, j - 1] = 1 - queue[tail] = i * W + (j - 1) + queue[tail] = np.int32((i << 16) | (j - 1)) tail += 1 if tail == cap: tail = 0 @@ -317,7 +351,7 @@ def reconstruction_by_dilation_2d(seed, mask): out[i, j + 1] = nv if inq[i, j + 1] == 0: inq[i, j + 1] = 1 - queue[tail] = i * W + (j + 1) + queue[tail] = np.int32((i << 16) | (j + 1)) tail += 1 if tail == cap: tail = 0 diff --git a/test/test_granularity_kernels.py b/test/test_granularity_kernels.py index 1aa15a7..1b0150d 100644 --- a/test/test_granularity_kernels.py +++ b/test/test_granularity_kernels.py @@ -129,3 +129,23 @@ def test_reconstruction_cascaded_mask_chain(): 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)