Skip to content

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

Draft
timtreis wants to merge 13 commits into
feat/bzyx-shapefrom
feat/intensity-bzyx
Draft

refactor(numba): migrate intensity backend onto (B,Z,Y,X)#57
timtreis wants to merge 13 commits into
feat/bzyx-shapefrom
feat/intensity-bzyx

Conversation

@timtreis

@timtreis timtreis commented Jun 3, 2026

Copy link
Copy Markdown
Collaborator

Stacked on #56 (needs to_bzyx; retarget to main once #56 merges).

Migrates the numba intensity backend onto the (B,Z,Y,X) batch shape from #56: get_intensity becomes a thin to_bzyx wrapper, per-volume work in _intensity_volume. Single → dict, 4D/list → list of dicts. Behaviour-preserving (numba-vs-numpy golden test green); single-image overhead ~0.26%. 180 passed.

timtreis and others added 10 commits June 2, 2026 04:47
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
- 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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
- _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) <noreply@anthropic.com>
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 (#55); unlike its
percentile/MAD changes, this one preserves output exactly (verified identical).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
feat(accelerator): numba intensity backend
@timtreis timtreis force-pushed the feat/accelerator-numba-granularity branch from 8d84546 to 13f1ee7 Compare June 3, 2026 18:55
@timtreis timtreis force-pushed the feat/intensity-bzyx branch from 5d08a15 to 6666ace Compare June 3, 2026 18:58
@timtreis timtreis changed the base branch from feat/accelerator-numba-granularity to feat/bzyx-shape June 3, 2026 18:58
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) <noreply@anthropic.com>
timtreis and others added 2 commits June 3, 2026 21:08
get_intensity becomes a thin batching wrapper over primitives.shapes.to_bzyx;
the per-volume computation moves to _intensity_volume. Single image/volume in =
one dict out; a 4D (B,Z,Y,X) array or a list of images = a list of dicts. Single
= batch-of-1 (one code path), matching the granularity backend.

Behaviour preserved for all real inputs (the numba-vs-numpy golden test stays
green for 2D and 3D): the 2D-mask-on-3D-stack broadcast passes through to_bzyx
unchanged, and the MAD fraction is derived as `2 if Z == 1 else 3`, which equals
the old `pixels.ndim` for every case except a degenerate single-slice volume
(Z == 1 treated as 2D -> 1/2 instead of 1/3; documented).

Tests: +single-returns-dict and +batch-matches-per-image (4D and list). Full
suite 180 passed; ruff clean.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
After the to_bzyx migration removed the reshape block, `masked_image` was just an
unconditional alias of `pixels`; use `pixels` directly. No behaviour change
(/simplify finding). Ruff clean, intensity tests pass.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
@timtreis timtreis force-pushed the feat/intensity-bzyx branch from 6666ace to c0f077e Compare June 3, 2026 19:09
@timtreis timtreis added the numba label Jun 9, 2026
timtreis added a commit that referenced this pull request Jun 9, 2026
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) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants