perf(sizeshape): scatter-based spatial moments + inertia, replacing regionprops einsum - #77
Conversation
1f77175 to
010232e
Compare
010232e to
303cd48
Compare
…apper Address /code-review notes on the sizeshape lane: - Extract `moment_feature_dict` into `primitives/_moments.py` as the single source of truth for the 53 `calculate_advanced` moment/inertia feature names and the (p,q) orders exposed. The numba `_sizeshape_2d` now calls it instead of building those keys with inline f-string loops, removing the duplication of the numpy assembly and the f-string/constant drift risk. The numpy `get_sizeshape` adopts the same helper when #77's moment rewrite rebases onto this lane (the sizeshape golden test cross-checks the keys vs the numpy F_* constants meanwhile). - Clarify the wrapper: `pixels` is accepted only for dispatch-signature parity and is unused (sizeshape is purely geometric, like the numpy backend); pass `masks` to `to_bzyx` for axis normalisation and drop the computed-then-discarded `_pixels_zyx`. - Add an end-to-end empty-mask test exercising the full dict assembly with zero objects (the kernels had empty tests; the wrapper assembly did not). 25 tests pass (incl. new empty case); ruff clean. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
6d961eb to
61b1a41
Compare
Benchmark —
|
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 1.03× | 1.05× | 1.01× |
| 512 | 0.99× | 0.98× | 0.94× |
| 1024 | 1.01× | 1.03× | 0.99× |
manders_fold
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 1.03× | 1.02× | 1.01× |
| 512 | 0.98× | 0.97× | 1.00× |
| 1024 | 0.97× | 0.94× | 1.00× |
pearson
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 1.03× | 0.99× | 0.99× |
| 512 | 0.99× | 0.98× | 0.99× |
| 1024 | 0.93× | 0.97× | 1.00× |
radial_zernikes
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 1.02× | 1.01× | 1.06× |
| 512 | 0.98× | 0.95× | 1.00× |
| 1024 | 1.02× | 1.03× | 1.04× |
sizeshape
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 1.25× | 1.39× | 1.48× |
| 512 | 1.08× | 1.24× | 1.42× |
| 1024 | 1.01× | 1.08× | 1.25× |
texture
| size \ objects | 16 | 64 | 256 |
|---|---|---|---|
| 256 | 0.96× | 0.98× | 0.99× |
| 512 | 0.98× | 0.97× | 0.95× |
| 1024 | 0.99× | 1.01× | 0.98× |
61b1a41 to
6d20661
Compare
…nprops einsum Replace regionprops' per-region moment einsum for the moments / central / normalized / Hu / inertia-tensor columns with batched separable matrix products (primitives/_moments.py): each object's moments are R.T @ mask @ C over its bounding box, computed for all objects at once as two BLAS matmuls — no per-object einsum path re-derivation, no full-image scatter. Axis lengths / eccentricity / orientation are derived from the same central moments (option B), so get_sizeshape requests no regionprops moments at all. Raw moments are bit-exact to regionprops; the rest match to round-off. Wins across image sizes and object counts, and the matmuls are tall-and-thin so BLAS runs them single-threaded. Assumes the contiguous 1..N label contract (see cp_measure._sanitize): find_objects index i is object i+1; per-object boxes are padded to the largest box. Also pins BLAS/OpenMP threads in the benchmark script so timings measure algorithmic cost, not incidental parallelism. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
6d20661 to
8db5df0
Compare
afermg
left a comment
There was a problem hiding this comment.
I requested some minor changes, mostly justifying the presence of global tolerance level and adjusting one to not assume a defined feature order. Once those are sorted out it should be an easy approval.
| # Thin / oblique objects have a (near-)singular inertia tensor; float error can drive the minor | ||
| # eigenvalue slightly negative. inertia_2d clips to 0 like skimage, so axis lengths and | ||
| # eccentricity match regionprops and are never NaN. Pre-clip, ~4% of these gave NaN axis_minor | ||
| # / eccentricity > 1 — this loop guards that regression. |
There was a problem hiding this comment.
Hmm, perhaps it is a bit excessive to test this, but I think it makes sense since we are ensuring that the original skimage clipping behaviour is maintained.
…set-based schema test - One relative tolerance (RTOL=1e-12) across all moment/axis asserts. Scale the absolute floor by the raw-moment magnitude of the same order, so the structurally-near-zero central moments stay robust on curved/asymmetric objects (their cancellation round-off is bounded by the summed terms, not their own ~0 value). Add a disk fixture that actually exercises this — plain axis-aligned rectangles hide it (their odd central moments are bit-exact 0). - Schema test asserts the feature *set* (order-independent), catching both a dropped feature and an accidental extra key. - Drop the unnecessary int() cast on labels.max(); reword the thin-object "near-singular" comment in plain terms. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
afermg
left a comment
There was a problem hiding this comment.
Propagating the global tolerances seems like the right choice. I'm still not fully convinced about the new test (test_moments_match_disk), but I think it's simple enough that it doesn't greatly increase complexity.
|
The test was to actually test a situation in which scaling can land "inbetween pixels" 👌 |
Replace regionprops' per-region einsum for the
moments/moments_central/moments_normalized/moments_hu/inertia_tensorcolumns with a single label-scatter pass (primitives/_moments.py), and derive axis lengths / eccentricity / orientation from the scatter central moments (option B) soget_sizeshaperequests no regionprops moments at all (~1.6×). Matches regionprops to round-off (raw moments bit-exact); the eigenvalue clip keeps thin objects NaN-free like skimage.Assumes the contiguous 1..N contract: segment index is
label - 1and per-object origins come from onescipy.ndimage.find_objectspass — no lookup table (this rebase drops thelabel_to_idx_lut/primitives/segment.pydependency that #89 removed).