feat(accelerator): numba texture (Haralick) backend#64
Draft
timtreis wants to merge 13 commits into
Draft
Conversation
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
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>
Replaces mahotas.features.haralick (the ~99% cost per the step-0 profile) with a fused per-object numba kernel: build the symmetric GLCM for each direction (bit-exact to mahotas.cooccurence), apply ignore_zeros, compute the 13 Haralick features. One kernel covers 2D (4 directions) and 3D (13 directions) over (dz,dy,dx) offsets = distance*delta. GLCM sized to crop.max()+1 (NOT 256) since feature 9 (px_minus_y.var()) is length-sensitive. Edge cases mirror mahotas exactly: empty GLCM after ignore_zeros -> whole-object NaN; sx==0 -> Correlation 1; max(HX,HY)==0; InfoMeas2 max(0,.) clamp; entropy 0*log2(1)=0. error_model="numpy". skimage img_as_ubyte/regionprops + gray_levels rescale stay host-side (~1%, exact). to_bzyx-shaped; texture is already per-object (regionprops crops) so it matches the numpy baseline directly (no Issue-#22 analogue). 4.8x (1938->402ms, 1080^2/144 obj). Tests: kernel vs mahotas (2D scale 1/3, 3D, empty->NaN, constant->Corr 1), backend golden vs numpy (2D/3D/batch/non-default scale+gray/empty), dispatch wiring. Full suite 100 passed, lint clean. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Two per-direction O(fm1^2 ~ 256^2) loops were the bottleneck. Both removed, bit-exact: - InfoMeas cross-entropies: for a SYMMETRIC GLCM, HXY1 == HXY2 == 2*HX (both collapse to the marginals; verified vs mahotas to ~1e-15). Replaces the fm1^2 double-loop with O(fm1). - T (GLCM total after ignore_zeros): accumulated during the build (2 per pair), then T = total - 2*row0 + cmat[0,0] in O(fm1) — no separate fm1^2 sum. - The per-cell feature pass now skips the background row/col 0 structurally (range(1, fm1)) instead of zeroing then re-scanning. 1984->118ms (16.8x, was 4.8x). Exactness unchanged: kernel vs mahotas + backend golden vs numpy all green; full suite 100 passed, lint clean. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Collaborator
Author
|
Follow-up perf (1f15132): collapsed the two per-direction O(fm1²≈256²) loops — 4.8× → 16.8× (1984→118 ms), bit-exact.
Kernel-vs-mahotas + backend golden all green; full suite 100. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stacked on #59 (
to_bzyx), sibling to #56/#57/#58/#60. Independent of the other lanes.Replaces
mahotas.features.haralick(the ~99% cost per the step-0 profile,tasks/profile_texture.py) with a fused per-object numba kernel (core/numba/_texture.py::haralick_object):mahotas.features.texture.cooccurence(integer histogram, algorithm-independent),ignore_zeros, compute the 13 Haralick features.One kernel covers 2D (4 directions) and 3D (13 directions) over
(dz,dy,dx)offsets =distance*delta. GLCM sized tocrop.max()+1(not 256) since feature 9 (px_minus_y.var()) is length-sensitive. Edge cases mirror mahotas exactly (empty GLCM → whole-object NaN;sx==0→ Correlation 1;max(HX,HY)==0; InfoMeas2max(0,·); entropy0·log2(1)=0).error_model="numpy".skimageimg_as_ubyte/regionprops+gray_levelsrescale stay host-side (~1%, exact).to_bzyx-shaped. No Issue-#22 analogue — texture is already per-object (regionprops crops), so it matches the numpy baseline directly.16.8× (1938→402 ms, 1080²/144 obj) — on the slowest feature, so ~1.5 s saved, bit-exact.
Tests: kernel vs mahotas (2D scale 1/3, 3D, empty→NaN, constant→Corr 1), backend golden vs numpy (2D/3D/batch/non-default scale+gray/empty), dispatch wiring. Full suite 100 passed, lint clean. Stack: #59 → this.