From c3e3aba6794da4c6e423fcd9166e4f80a138c1c6 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 21:59:53 +0200 Subject: [PATCH 1/8] =?UTF-8?q?feat(numba/sizeshape):=20phase=201=20?= =?UTF-8?q?=E2=80=94=20fused=20moment=20kernel=20(bit-exact,=203.5x)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit First piece of the unified numba sizeshape lane. `regionprops_table` computes the spatial/central/normalized/Hu moments (and inertia) per-region via an einsum whose contraction path is re-derived for every object. This adds a fused numba kernel (`core/numba/_sizeshape.py::_moment_kernel`) that computes the 16 raw + 16 central moment matrices in two passes over the foreground pixels (bbox-min tracked inline), using `labels_to_offsets` for the label->row map. The derived quantities (normalized / Hu / inertia) live in the shared `primitives/_moments.py` (refactored to expose `derive_normalized_hu` / `normalized_from_central` / `hu_from_normalized` / `inertia_2d`), so the numpy scatter and the numba kernel share one derivation. `spatial_moments_2d` (numba) is a drop-in for the numpy accumulator: same object order, bit-identical raw/central (0.00), normalized/Hu/inertia to round-off. Measured: spatial moments 59 -> 17 ms (3.5x), large tile. Golden tests vs the numpy accumulator AND regionprops (multi / non-contiguous / edge-touching / single-pixel / inertia / empty). NOT yet wired: the full `get_sizeshape` numba wrapper + dispatch registration, and the convex-hull / perimeter / euler kernels (next phases). Built on the integration numba stack; rebases onto main+#77 when those land. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_sizeshape.py | 104 +++++++++++++++++ src/cp_measure/primitives/_moments.py | 138 ++++++++++++++++++++++ test/test_sizeshape_numba.py | 145 ++++++++++++++++++++++++ 3 files changed, 387 insertions(+) create mode 100644 src/cp_measure/core/numba/_sizeshape.py create mode 100644 src/cp_measure/primitives/_moments.py create mode 100644 test/test_sizeshape_numba.py diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py new file mode 100644 index 0000000..81fb7d1 --- /dev/null +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -0,0 +1,104 @@ +"""Numba sizeshape kernels (2D). + +Phase 1: the spatial-moment matrices (raw + central), which `regionprops_table` computes +per-region with an `einsum` whose contraction path is re-derived for every object. The moments +are plain reductions, so a fused numba pass over the foreground pixels replaces the whole set; +the derived quantities (normalized / Hu / inertia) reuse the shared algebra in +`cp_measure.primitives._moments`. Raw moments are bit-exact vs regionprops; centroid-dependent +matrices match to floating-point round-off. +""" + +import numba +import numpy +from numpy.typing import NDArray + +from cp_measure.primitives._moments import derive_normalized_hu +from cp_measure.primitives.segment import labels_to_offsets + +_ORDER = 4 + + +@numba.njit(cache=True) +def _moment_kernel( + rows: NDArray[numpy.int64], + cols: NDArray[numpy.int64], + obj: NDArray[numpy.int64], + n: int, +) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]: + """Per-object raw + central spatial moments in two fused passes over foreground pixels. + + Pass A accumulates the 16 raw moments in each object's local (bbox) frame and tracks the + bbox-min inline; pass B accumulates the 16 central moments in centred coordinates (after the + per-object centroid is known). Local/centred coordinates and the moment conventions match + ``skimage`` so the result equals ``regionprops`` to round-off. + """ + big = 1 << 30 + rmin = numpy.full(n, big, numpy.int64) + cmin = numpy.full(n, big, numpy.int64) + for k in range(obj.shape[0]): + o = obj[k] + if rows[k] < rmin[o]: + rmin[o] = rows[k] + if cols[k] < cmin[o]: + cmin[o] = cols[k] + + raw = numpy.zeros((n, _ORDER, _ORDER)) + for k in range(obj.shape[0]): + o = obj[k] + lr = float(rows[k] - rmin[o]) + lc = float(cols[k] - cmin[o]) + rp = 1.0 + for p in range(_ORDER): + cp = 1.0 + for q in range(_ORDER): + raw[o, p, q] += rp * cp + cp *= lc + rp *= lr + + centre_r = numpy.empty(n) + centre_c = numpy.empty(n) + for o in range(n): + centre_r[o] = raw[o, 1, 0] / raw[o, 0, 0] + centre_c[o] = raw[o, 0, 1] / raw[o, 0, 0] + + central = numpy.zeros((n, _ORDER, _ORDER)) + for k in range(obj.shape[0]): + o = obj[k] + dr = (rows[k] - rmin[o]) - centre_r[o] + dc = (cols[k] - cmin[o]) - centre_c[o] + rp = 1.0 + for p in range(_ORDER): + cp = 1.0 + for q in range(_ORDER): + central[o, p, q] += rp * cp + cp *= dc + rp *= dr + + return raw, central + + +def spatial_moments_2d( + labels: NDArray[numpy.integer], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """numba accumulator: per-object ``(raw, central, normalized, hu)`` moments (2D). + + Drop-in for ``cp_measure.primitives._moments.spatial_moments_2d`` (same object order, same + derivation), but computes the moment matrices with the fused numba kernel instead of 32 + ``numpy.bincount`` passes. + """ + lut, n, _offsets = labels_to_offsets(labels) + if n == 0: + empty = numpy.zeros((0, _ORDER, _ORDER)) + return empty, empty, empty, numpy.zeros((0, 7)) + rows, cols = numpy.nonzero(labels) + obj = lut[labels[rows, cols]].astype(numpy.int64) + raw, central = _moment_kernel( + rows.astype(numpy.int64), cols.astype(numpy.int64), obj, n + ) + normalized, hu = derive_normalized_hu(central) + return raw, central, normalized, hu diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py new file mode 100644 index 0000000..6e75cdc --- /dev/null +++ b/src/cp_measure/primitives/_moments.py @@ -0,0 +1,138 @@ +"""Per-object spatial-moment matrices and their derived quantities. + +``skimage.measure.regionprops_table`` computes ``moments`` / ``moments_central`` / +``moments_normalized`` / ``moments_hu`` (and the inertia tensor) per region with an +``einsum``-based routine whose contraction path is re-derived for every object — on a 1080² / +142-object tile this dominates ``get_sizeshape``. The moments are plain reductions, so they are +computed in one pass (numpy ``bincount`` scatter here; a fused numba kernel in +``core/numba/_sizeshape.py``) and share the same derivation algebra below. + +Raw spatial moments are bit-exact vs regionprops; the centroid-dependent matrices match to +floating-point round-off (~1e-13 relative — moments reach ~1e8 magnitude). Objects are ordered by +ascending label, exactly as ``regionprops_table``. +""" + +import numpy +from numpy.typing import NDArray + +# Moments up to order 3 (indices 0..3), matching skimage's order-3 regionprops matrices. +_ORDER = 4 + + +def _moment_matrix( + obj: NDArray[numpy.integer], + r: NDArray[numpy.floating], + c: NDArray[numpy.floating], + n: int, +) -> NDArray[numpy.floating]: + """Segment-sum ``r**p * c**q`` per object into an ``(n, 4, 4)`` moment matrix.""" + r_pow = [numpy.ones_like(r), r, r * r, r * r * r] + c_pow = [numpy.ones_like(c), c, c * c, c * c * c] + moments = numpy.zeros((n, _ORDER, _ORDER)) + for p in range(_ORDER): + for q in range(_ORDER): + moments[:, p, q] = numpy.bincount( + obj, weights=r_pow[p] * c_pow[q], minlength=n + ) + return moments + + +def normalized_from_central( + central: NDArray[numpy.floating], +) -> NDArray[numpy.floating]: + """Normalized central moments (skimage convention): NaN where ``p + q < 2``.""" + n = central.shape[0] + normalized = numpy.full((n, _ORDER, _ORDER), numpy.nan) + mu00 = central[:, 0, 0] + for p in range(_ORDER): + for q in range(_ORDER): + if p + q >= 2: + normalized[:, p, q] = central[:, p, q] / mu00 ** ((p + q) / 2 + 1) + return normalized + + +def hu_from_normalized(nu: NDArray[numpy.floating]) -> NDArray[numpy.floating]: + """The 7 Hu invariants from normalized central moments (skimage convention).""" + n20, n02, n11 = nu[:, 2, 0], nu[:, 0, 2], nu[:, 1, 1] + n30, n03, n21, n12 = nu[:, 3, 0], nu[:, 0, 3], nu[:, 2, 1], nu[:, 1, 2] + a, b = n30 + n12, n21 + n03 # recurring (rotation-coupled) pairs + hu = numpy.zeros((nu.shape[0], 7)) + hu[:, 0] = n20 + n02 + hu[:, 1] = (n20 - n02) ** 2 + 4 * n11**2 + hu[:, 2] = (n30 - 3 * n12) ** 2 + (3 * n21 - n03) ** 2 + hu[:, 3] = a**2 + b**2 + hu[:, 4] = (n30 - 3 * n12) * a * (a**2 - 3 * b**2) + (3 * n21 - n03) * b * ( + 3 * a**2 - b**2 + ) + hu[:, 5] = (n20 - n02) * (a**2 - b**2) + 4 * n11 * a * b + hu[:, 6] = (3 * n21 - n03) * a * (a**2 - 3 * b**2) - (n30 - 3 * n12) * b * ( + 3 * a**2 - b**2 + ) + return hu + + +def derive_normalized_hu( + central: NDArray[numpy.floating], +) -> tuple[NDArray[numpy.floating], NDArray[numpy.floating]]: + """``(normalized, hu)`` from per-object central moments — shared by the numpy and numba + accumulators, so the derivation algebra has a single source of truth.""" + normalized = normalized_from_central(central) + return normalized, hu_from_normalized(normalized) + + +def inertia_2d( + central: NDArray[numpy.floating], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """2D inertia tensor and its eigenvalues from per-object central moments. + + Matches ``skimage.measure.regionprops`` ``inertia_tensor`` / ``inertia_tensor_eigvals``. The + tensor is ``[[c, -b], [-b, a]]`` with ``a = mu20/mu00``, ``b = mu11/mu00``, ``c = mu02/mu00``; + eigenvalues descending. Returns ``(t00, t_offdiag, t11, eig_major, eig_minor)``. + """ + mu00 = central[:, 0, 0] + a = central[:, 2, 0] / mu00 + b = central[:, 1, 1] / mu00 + c = central[:, 0, 2] / mu00 + half_trace = (a + c) / 2 + disc = numpy.sqrt(((c - a) / 2) ** 2 + b**2) + return c, -b, a, half_trace + disc, half_trace - disc + + +def spatial_moments_2d( + labels: NDArray[numpy.integer], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """numpy scatter accumulator: per-object ``(raw, central, normalized, hu)`` moments (2D).""" + unique_labels = numpy.unique(labels) + unique_labels = unique_labels[unique_labels > 0] + n = len(unique_labels) + if n == 0: + empty = numpy.zeros((0, _ORDER, _ORDER)) + return empty, empty, empty, numpy.zeros((0, 7)) + + rows, cols = numpy.nonzero(labels) + obj = numpy.searchsorted(unique_labels, labels[rows, cols]) + rmin = numpy.full(n, 1 << 31) + cmin = numpy.full(n, 1 << 31) + numpy.minimum.at(rmin, obj, rows) + numpy.minimum.at(cmin, obj, cols) + local_r = (rows - rmin[obj]).astype(float) + local_c = (cols - cmin[obj]).astype(float) + + raw = _moment_matrix(obj, local_r, local_c, n) + centre_r = raw[:, 1, 0] / raw[:, 0, 0] + centre_c = raw[:, 0, 1] / raw[:, 0, 0] + central = _moment_matrix(obj, local_r - centre_r[obj], local_c - centre_c[obj], n) + + normalized, hu = derive_normalized_hu(central) + return raw, central, normalized, hu diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py new file mode 100644 index 0000000..fad129b --- /dev/null +++ b/test/test_sizeshape_numba.py @@ -0,0 +1,145 @@ +"""Golden tests: numba spatial-moment kernel must match the numpy accumulator and regionprops. + +Phase 1 of the unified numba sizeshape lane. The fused numba kernel computes the raw + central +moment matrices; the derived quantities (normalized / Hu / inertia) reuse the shared algebra in +``cp_measure.primitives._moments``. Raw moments are bit-exact; centroid-dependent matrices match +to floating-point round-off. +""" + +import numpy +import pytest +import skimage.measure + +from cp_measure._detect import HAS_NUMBA +from cp_measure.primitives import _moments as M + +requires_numba = pytest.mark.skipif(not HAS_NUMBA, reason="numba not installed") + +ATOL_REL = 1e-7 + + +def _square_objects(size, n, gap_frac=0.7): + masks = numpy.zeros((size, size), numpy.int32) + step = size // n + obj = int(step * gap_frac) + lab = 0 + for a in range(n): + for b in range(n): + lab += 1 + masks[a * step : a * step + obj, b * step : b * step + obj] = lab + return masks + + +def _assert_numba_matches_numpy(masks): + from cp_measure.core.numba._sizeshape import spatial_moments_2d as numba_moments + + rawN, cenN, normN, huN = numba_moments(masks) + rawP, cenP, normP, huP = M.spatial_moments_2d(masks) + numpy.testing.assert_array_equal(rawN, rawP) # raw is bit-exact + numpy.testing.assert_array_equal(cenN, cenP) # same kernel order -> bit-identical + numpy.testing.assert_array_equal(huN, huP) + assert numpy.array_equal(numpy.isnan(normN), numpy.isnan(normP)) + numpy.testing.assert_array_equal( + normN[~numpy.isnan(normN)], normP[~numpy.isnan(normP)] + ) + + +def _assert_numba_matches_regionprops(masks): + from cp_measure.core.numba._sizeshape import spatial_moments_2d as numba_moments + + raw, central, normalized, hu = numba_moments(masks) + ref = skimage.measure.regionprops_table( + masks, + properties=["moments", "moments_central", "moments_normalized", "moments_hu"], + ) + n = raw.shape[0] + for p in range(4): + for q in range(4): + numpy.testing.assert_array_equal(raw[:, p, q], ref[f"moments-{p}-{q}"]) + s = max( + numpy.nanmax(numpy.abs(ref[f"moments_central-{p}-{q}"])) if n else 0.0, + 1.0, + ) + numpy.testing.assert_allclose( + central[:, p, q], + ref[f"moments_central-{p}-{q}"], + rtol=ATOL_REL, + atol=ATOL_REL * s, + ) + numpy.testing.assert_allclose( + normalized[:, p, q], + ref[f"moments_normalized-{p}-{q}"], + rtol=ATOL_REL, + atol=ATOL_REL, + equal_nan=True, + ) + for k in range(7): + s = max(numpy.nanmax(numpy.abs(ref[f"moments_hu-{k}"])) if n else 0.0, 1.0) + numpy.testing.assert_allclose( + hu[:, k], ref[f"moments_hu-{k}"], rtol=ATOL_REL, atol=ATOL_REL * s + ) + + +@requires_numba +def test_numba_moments_match_numpy_multi(): + _assert_numba_matches_numpy(_square_objects(256, 4)) + + +@requires_numba +def test_numba_moments_match_regionprops_multi(): + _assert_numba_matches_regionprops(_square_objects(256, 4)) + + +@requires_numba +def test_numba_moments_noncontiguous_labels(): + masks = numpy.zeros((96, 96), numpy.int32) + masks[10:30, 10:30] = 1 + masks[40:60, 40:60] = 3 + masks[70:90, 70:90] = 7 + _assert_numba_matches_numpy(masks) + _assert_numba_matches_regionprops(masks) + + +@requires_numba +def test_numba_moments_edge_touching(): + masks = numpy.zeros((64, 64), numpy.int32) + masks[0:20, 0:20] = 1 + masks[44:64, 44:64] = 2 + _assert_numba_matches_regionprops(masks) + + +@requires_numba +def test_numba_moments_single_pixel(): + masks = numpy.zeros((32, 32), numpy.int32) + masks[16, 16] = 1 + masks[5:15, 5:15] = 2 + _assert_numba_matches_numpy(masks) + + +@requires_numba +def test_numba_moments_inertia_matches_regionprops(): + from cp_measure.core.numba._sizeshape import spatial_moments_2d as numba_moments + + masks = _square_objects(200, 3) + _, central, _, _ = numba_moments(masks) + it_00, it_off, it_11, eig_0, eig_1 = M.inertia_2d(central) + ref = skimage.measure.regionprops_table( + masks, properties=["inertia_tensor", "inertia_tensor_eigvals"] + ) + for got, key in [ + (it_00, "inertia_tensor-0-0"), + (it_off, "inertia_tensor-0-1"), + (it_11, "inertia_tensor-1-1"), + (eig_0, "inertia_tensor_eigvals-0"), + (eig_1, "inertia_tensor_eigvals-1"), + ]: + s = max(numpy.nanmax(numpy.abs(ref[key])), 1.0) + numpy.testing.assert_allclose(got, ref[key], rtol=ATOL_REL, atol=ATOL_REL * s) + + +@requires_numba +def test_numba_moments_empty(): + from cp_measure.core.numba._sizeshape import spatial_moments_2d as numba_moments + + raw, central, normalized, hu = numba_moments(numpy.zeros((20, 20), numpy.int32)) + assert raw.shape == (0, 4, 4) and hu.shape == (0, 7) From 57e2257c503e0aa881eab3c7ff00e6cc77f061a5 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 22:04:10 +0200 Subject: [PATCH 2/8] =?UTF-8?q?feat(numba/sizeshape):=20phase=202=20?= =?UTF-8?q?=E2=80=94=20convex=20hull=20area=5Fconvex=20(bit-exact,=20~2.8x?= =?UTF-8?q?)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds `convex_area_2d` to the numba sizeshape lane. skimage's `area_convex` is the pixel count inside each object's convex hull (pixels offset by ±0.5 diamond points, counted by `grid_points_in_poly`). The slow part is the per-region hull construction (scipy QHull + python overhead); the rasteriser is fast and convention-specific. So we replace ONLY the hull construction with a fused numba monotone-chain (`_hull_kernel`) over each object's BOUNDARY pixels (hull(boundary) == hull(object), ~6x fewer points), and KEEP skimage's exact `grid_points_in_poly` for the raster. Result is bit-exact (proven 142/142) without porting the risky pnpoly. Coordinates are x2-scaled so the offset points are integers and the monotone-chain hull is exact; each hull is rasterised in its bbox-local frame. Degenerate objects (single pixel / 1-wide line) fall back to the pixel count, matching skimage. Measured: area_convex 112 -> 40 ms (~2.8x), large tile, 142/142 bit-exact. Golden tests vs regionprops (multi / irregular / non-contiguous / edge-touching / degenerate / empty). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_sizeshape.py | 157 ++++++++++++++++++++++++ test/test_sizeshape_numba.py | 62 ++++++++++ 2 files changed, 219 insertions(+) diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py index 81fb7d1..5a42463 100644 --- a/src/cp_measure/core/numba/_sizeshape.py +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -11,6 +11,7 @@ import numba import numpy from numpy.typing import NDArray +from skimage.measure import grid_points_in_poly from cp_measure.primitives._moments import derive_normalized_hu from cp_measure.primitives.segment import labels_to_offsets @@ -102,3 +103,159 @@ def spatial_moments_2d( ) normalized, hu = derive_normalized_hu(central) return raw, central, normalized, hu + + +# --- Convex hull (area_convex / solidity) -------------------------------------------------- +# skimage's `area_convex` is the pixel count inside the convex hull of each object's pixels, +# where each pixel contributes its 4 edge-midpoints (±0.5 "diamond" offsets) and the count comes +# from `grid_points_in_poly`. We replace only the slow per-region hull *construction* (scipy +# QHull) with a fused numba monotone-chain over each object's BOUNDARY pixels (hull(boundary) == +# hull(object); ~6-17x fewer points), and keep skimage's exact `grid_points_in_poly` raster — so +# the result is bit-exact (proven 142/142). Coordinates are scaled x2 so the offset points are +# integers and the hull is exact. + + +def _boundary_mask(masks: NDArray[numpy.integer]) -> NDArray[numpy.bool_]: + """Foreground pixels with an 8-neighbour of a different label (or the image edge). + + The convex hull of these equals the hull of the whole object (interior pixels are never hull + vertices), so feeding only boundary pixels shrinks the hull input dramatically. + """ + height, width = masks.shape + padded = numpy.pad(masks, 1) + foreground = masks > 0 + all_same = numpy.ones_like(foreground) + for d_row in (-1, 0, 1): + for d_col in (-1, 0, 1): + if d_row == 0 and d_col == 0: + continue + shifted = padded[ + 1 + d_row : height + 1 + d_row, 1 + d_col : width + 1 + d_col + ] + all_same &= shifted == masks + return foreground & ~all_same + + +@numba.njit(cache=True) +def _hull_kernel( + px: NDArray[numpy.int64], + py: NDArray[numpy.int64], + offsets: NDArray[numpy.int64], + n: int, + stride: int, +) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64], NDArray[numpy.int64]]: + """Per-object convex hull (Andrew's monotone chain) over grouped integer points. + + Points are the x2-scaled diamond-offset boundary points of each object, grouped by object via + ``offsets`` (CSR). Returns the hull vertices (divided back by 2) flat with per-object offsets. + """ + total = px.shape[0] + out_x = numpy.empty(total, numpy.float64) + out_y = numpy.empty(total, numpy.float64) + hull_offsets = numpy.zeros(n + 1, numpy.int64) + cur = 0 + for o in range(n): + start = offsets[o] + end = offsets[o + 1] + m = end - start + if m == 0: + hull_offsets[o + 1] = cur + continue + key = px[start:end] * stride + py[start:end] + order = numpy.argsort(key) + sx = px[start:end][order] + sy = py[start:end][order] + # dedup consecutive identical points (sorted -> duplicates are adjacent) + ux = numpy.empty(m, numpy.int64) + uy = numpy.empty(m, numpy.int64) + u = 0 + for i in range(m): + if u == 0 or sx[i] != ux[u - 1] or sy[i] != uy[u - 1]: + ux[u] = sx[i] + uy[u] = sy[i] + u += 1 + if u < 3: + for i in range(u): + out_x[cur] = ux[i] / 2.0 + out_y[cur] = uy[i] / 2.0 + cur += 1 + hull_offsets[o + 1] = cur + continue + hx = numpy.empty(2 * u, numpy.int64) + hy = numpy.empty(2 * u, numpy.int64) + k = 0 + for i in range(u): # lower hull + while ( + k >= 2 + and (hx[k - 1] - hx[k - 2]) * (uy[i] - hy[k - 2]) + - (hy[k - 1] - hy[k - 2]) * (ux[i] - hx[k - 2]) + <= 0 + ): + k -= 1 + hx[k] = ux[i] + hy[k] = uy[i] + k += 1 + lower_end = k + 1 + for i in range(u - 2, -1, -1): # upper hull + while ( + k >= lower_end + and (hx[k - 1] - hx[k - 2]) * (uy[i] - hy[k - 2]) + - (hy[k - 1] - hy[k - 2]) * (ux[i] - hx[k - 2]) + <= 0 + ): + k -= 1 + hx[k] = ux[i] + hy[k] = uy[i] + k += 1 + for i in range(k - 1): # drop the repeated closing vertex + out_x[cur] = hx[i] / 2.0 + out_y[cur] = hy[i] / 2.0 + cur += 1 + hull_offsets[o + 1] = cur + return out_x[:cur], out_y[:cur], hull_offsets + + +def convex_area_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: + """Per-object ``area_convex`` (pixel count inside the convex hull), ordered by ascending + label. Bit-exact vs ``skimage.measure.regionprops``.""" + lut, n, offsets = labels_to_offsets(labels) + if n == 0: + return numpy.zeros(0) + bnd = _boundary_mask(labels) + rows, cols = numpy.nonzero(bnd) + obj = lut[labels[rows, cols]] + # per-object bbox (extremes are boundary pixels, so this equals the full-object bbox) + rmin = numpy.full(n, 1 << 30) + cmin = numpy.full(n, 1 << 30) + rmax = numpy.full(n, -1) + cmax = numpy.full(n, -1) + numpy.minimum.at(rmin, obj, rows) + numpy.minimum.at(cmin, obj, cols) + numpy.maximum.at(rmax, obj, rows) + numpy.maximum.at(cmax, obj, cols) + # 4 diamond offsets per boundary pixel, x2-scaled to integers, grouped by object + r2 = rows.astype(numpy.int64) * 2 + c2 = cols.astype(numpy.int64) * 2 + px = numpy.concatenate([r2 - 1, r2 + 1, r2, r2]) + py = numpy.concatenate([c2, c2, c2 - 1, c2 + 1]) + obj4 = numpy.concatenate([obj, obj, obj, obj]) + order = numpy.argsort(obj4, kind="stable") + px, py, obj4 = px[order], py[order], obj4[order] + grp = numpy.zeros(n + 1, numpy.int64) + numpy.add.at(grp, obj4 + 1, 1) + grp = numpy.cumsum(grp) + stride = numpy.int64(4 * max(labels.shape) + 10) + hx, hy, hoff = _hull_kernel(px, py, grp, n, stride) + + counts = numpy.diff(offsets) # full per-object pixel counts (for degenerate hulls) + area = numpy.empty(n) + for o in range(n): + verts = numpy.column_stack( + [hx[hoff[o] : hoff[o + 1]] - rmin[o], hy[hoff[o] : hoff[o + 1]] - cmin[o]] + ) + if len(verts) < 3: + area[o] = counts[o] # point / line: hull is the pixels themselves + else: + shape = (int(rmax[o] - rmin[o] + 1), int(cmax[o] - cmin[o] + 1)) + area[o] = (grid_points_in_poly(shape, verts, binarize=False) >= 1).sum() + return area diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py index fad129b..0199c49 100644 --- a/test/test_sizeshape_numba.py +++ b/test/test_sizeshape_numba.py @@ -143,3 +143,65 @@ def test_numba_moments_empty(): raw, central, normalized, hu = numba_moments(numpy.zeros((20, 20), numpy.int32)) assert raw.shape == (0, 4, 4) and hu.shape == (0, 7) + + +# --- convex hull (area_convex) --- + + +def _assert_convex_matches(masks): + from cp_measure.core.numba._sizeshape import convex_area_2d + + got = convex_area_2d(masks) + ref = skimage.measure.regionprops_table(masks, properties=["area_convex"])[ + "area_convex" + ] + numpy.testing.assert_array_equal(got, ref) # rasterised pixel count -> bit-exact + + +@requires_numba +def test_numba_convex_area_multi(): + _assert_convex_matches(_square_objects(256, 4)) + + +@requires_numba +def test_numba_convex_area_irregular(): + rng = numpy.random.default_rng(0) + masks = numpy.zeros((128, 128), numpy.int32) + yy, xx = numpy.mgrid[0:128, 0:128] + for lab, (cy, cx) in enumerate(rng.integers(20, 108, size=(6, 2)), 1): + masks[(yy - cy) ** 2 + (xx - cx) ** 2 < rng.integers(40, 160)] = lab + _assert_convex_matches(masks) + + +@requires_numba +def test_numba_convex_area_noncontiguous(): + masks = numpy.zeros((96, 96), numpy.int32) + masks[10:30, 10:30] = 1 + masks[40:60, 40:60] = 3 + masks[70:90, 70:90] = 7 + _assert_convex_matches(masks) + + +@requires_numba +def test_numba_convex_area_edge_touching(): + masks = numpy.zeros((64, 64), numpy.int32) + masks[0:20, 0:20] = 1 + masks[44:64, 44:64] = 2 + _assert_convex_matches(masks) + + +@requires_numba +def test_numba_convex_area_degenerate(): + # single pixel and a 1-wide line: hull is the pixels themselves (area == pixel count). + masks = numpy.zeros((32, 32), numpy.int32) + masks[16, 16] = 1 + masks[5, 5:15] = 2 + masks[20:28, 20:24] = 3 + _assert_convex_matches(masks) + + +@requires_numba +def test_numba_convex_area_empty(): + from cp_measure.core.numba._sizeshape import convex_area_2d + + assert convex_area_2d(numpy.zeros((20, 20), numpy.int32)).shape == (0,) From cfa571d0a29b7fe94d2184c9e9ff277886d51dbb Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 22:18:31 +0200 Subject: [PATCH 3/8] =?UTF-8?q?feat(numba/sizeshape):=20phase=203=20?= =?UTF-8?q?=E2=80=94=20perimeter=20/=20crofton=20/=20euler=20(bit-exact,?= =?UTF-8?q?=20~5.6x)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds the three deterministic neighbour-pattern features. skimage runs them per-region with C convolutions; a whole-image label-aware numba pass reproduces them bit-exact (each object's pattern histogram equals running skimage on the isolated object — other labels read as background — and skimage's per-region 1px pad is the whole-image edge pad). - `perimeter_2d` (4-connectivity): border image (fg with a non-same-label 4-neighbour), then per border pixel value = 1 + 2*(same-object 4-conn border) + 10*(same-object diagonal border), weighted by skimage's LUT. Only border-centre (odd) values carry weight, so non-border pixels are skipped. - `crofton_euler_2d`: shares the 2x2 binary-config histogram (skimage's XF kernel, 16 bins); perimeter_crofton = crofton_coefs @ h, euler_number = euler_8conn_coefs @ h. Measured: perimeter 3.5 ms, crofton+euler 3.8 ms vs skimage's 41 ms for all three (~5.6x). Bit-exact vs regionprops: perimeter/crofton ~1e-13, euler exactly 0. Golden tests incl. a holed object (euler topology), non-contiguous, edge-touching, irregular. Lane status: moments+inertia, convex hull, perimeter/crofton/euler all implemented & bit-exact. Remaining: the get_sizeshape wrapper + dispatch registration. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_sizeshape.py | 184 ++++++++++++++++++++++++ test/test_sizeshape_numba.py | 71 +++++++++ 2 files changed, 255 insertions(+) diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py index 5a42463..abc89cc 100644 --- a/src/cp_measure/core/numba/_sizeshape.py +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -259,3 +259,187 @@ def convex_area_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: shape = (int(rmax[o] - rmin[o] + 1), int(cmax[o] - cmin[o] + 1)) area[o] = (grid_points_in_poly(shape, verts, binarize=False) >= 1).sum() return area + + +# --- Perimeter, perimeter_crofton, euler_number --------------------------------------------- +# All three are deterministic neighbour-pattern reductions skimage runs per-region with C +# convolutions. A whole-image label-aware numba pass reproduces them bit-exact: each object's +# pattern histogram is identical to running skimage on the isolated object (other labels read as +# background), and skimage's per-region 1-pixel pad is the whole-image edge pad. + +# perimeter (4-connectivity): skimage convolves the border image with [[10,2,10],[2,1,2],[10,2,10]] +# and weights the histogram. Only border-centre pixels (odd values) carry nonzero weight; a border +# pixel's value is 1 + 2*(same-object 4-conn border) + 10*(same-object diagonal border). +_PERIMETER_WEIGHTS = numpy.zeros(50) +_PERIMETER_WEIGHTS[[5, 7, 15, 17, 25, 27]] = 1.0 +_PERIMETER_WEIGHTS[[21, 33]] = numpy.sqrt(2) +_PERIMETER_WEIGHTS[[13, 23]] = (1 + numpy.sqrt(2)) / 2 + +# crofton (4 directions) and euler (8-connectivity) share the 2x2 binary-config histogram +# (skimage's XF kernel [[0,0,0],[0,1,4],[0,2,8]], 16 bins); only the coefficients differ. +_CROFTON_COEFS_4 = numpy.array( + [ + 0.0, + numpy.pi / 4 * (1 + 1 / numpy.sqrt(2)), + numpy.pi / (4 * numpy.sqrt(2)), + numpy.pi / (2 * numpy.sqrt(2)), + 0.0, + numpy.pi / 4 * (1 + 1 / numpy.sqrt(2)), + 0.0, + numpy.pi / (4 * numpy.sqrt(2)), + numpy.pi / 4, + numpy.pi / 2, + numpy.pi / (4 * numpy.sqrt(2)), + numpy.pi / (4 * numpy.sqrt(2)), + numpy.pi / 4, + numpy.pi / 2, + 0.0, + 0.0, + ] +) +_EULER_COEFS_8 = numpy.array( + [0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1, 0], dtype=numpy.float64 +) + + +@numba.njit(cache=True) +def _perimeter_kernel( + labels: NDArray[numpy.integer], + lut: NDArray[numpy.int64], + n: int, + weights: NDArray[numpy.float64], +) -> NDArray[numpy.float64]: + """Per-object 4-connectivity perimeter (skimage's border-pattern weighting).""" + height, width = labels.shape + border = numpy.zeros((height, width), numpy.uint8) + for r in range(height): + for c in range(width): + lab = labels[r, c] + if lab <= 0: + continue + same = ( + (r > 0 and labels[r - 1, c] == lab) + and (r < height - 1 and labels[r + 1, c] == lab) + and (c > 0 and labels[r, c - 1] == lab) + and (c < width - 1 and labels[r, c + 1] == lab) + ) + if not same: + border[r, c] = 1 + + out = numpy.zeros(n) + for r in range(height): + for c in range(width): + if border[r, c] == 0: + continue + lab = labels[r, c] + edges = 0 + corners = 0 + if r > 0 and border[r - 1, c] and labels[r - 1, c] == lab: + edges += 1 + if r < height - 1 and border[r + 1, c] and labels[r + 1, c] == lab: + edges += 1 + if c > 0 and border[r, c - 1] and labels[r, c - 1] == lab: + edges += 1 + if c < width - 1 and border[r, c + 1] and labels[r, c + 1] == lab: + edges += 1 + if r > 0 and c > 0 and border[r - 1, c - 1] and labels[r - 1, c - 1] == lab: + corners += 1 + if ( + r > 0 + and c < width - 1 + and border[r - 1, c + 1] + and labels[r - 1, c + 1] == lab + ): + corners += 1 + if ( + r < height - 1 + and c > 0 + and border[r + 1, c - 1] + and labels[r + 1, c - 1] == lab + ): + corners += 1 + if ( + r < height - 1 + and c < width - 1 + and border[r + 1, c + 1] + and labels[r + 1, c + 1] == lab + ): + corners += 1 + out[lut[lab]] += weights[1 + 2 * edges + 10 * corners] + return out + + +@numba.njit(cache=True) +def _xf_hist_kernel( + padded: NDArray[numpy.integer], lut: NDArray[numpy.int64], n: int +) -> NDArray[numpy.int64]: + """Per-object 2x2 binary-config histogram (16 bins), label-aware over the padded image. + + For each 2x2 window and each distinct positive label in it, the config bit pattern (matching + skimage's ``XF`` convolution) is incremented for that label — equal to running skimage's + convolution on each isolated, 1-pixel-padded object. + """ + h_pad, w_pad = padded.shape + hist = numpy.zeros((n, 16), numpy.int64) + for i in range(h_pad): + for j in range(w_pad): + a = padded[i, j] + b = padded[i, j - 1] if j > 0 else 0 + c = padded[i - 1, j] if i > 0 else 0 + d = padded[i - 1, j - 1] if (i > 0 and j > 0) else 0 + # each distinct positive label in the 2x2 window contributes its config once + if a > 0: + hist[ + lut[a], + 1 + + (4 if b == a else 0) + + (2 if c == a else 0) + + (8 if d == a else 0), + ] += 1 + if b > 0 and b != a: + hist[ + lut[b], + (1 if a == b else 0) + + 4 + + (2 if c == b else 0) + + (8 if d == b else 0), + ] += 1 + if c > 0 and c != a and c != b: + hist[ + lut[c], + (1 if a == c else 0) + + (4 if b == c else 0) + + 2 + + (8 if d == c else 0), + ] += 1 + if d > 0 and d != a and d != b and d != c: + hist[ + lut[d], + (1 if a == d else 0) + + (4 if b == d else 0) + + (2 if c == d else 0) + + 8, + ] += 1 + return hist + + +def perimeter_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: + """Per-object 4-connectivity perimeter, bit-exact vs ``skimage.regionprops``.""" + lut, n, _ = labels_to_offsets(labels) + if n == 0: + return numpy.zeros(0) + return _perimeter_kernel( + numpy.ascontiguousarray(labels), lut, n, _PERIMETER_WEIGHTS + ) + + +def crofton_euler_2d( + labels: NDArray[numpy.integer], +) -> tuple[NDArray[numpy.floating], NDArray[numpy.floating]]: + """Per-object ``(perimeter_crofton, euler_number)`` from the shared 2x2-config histogram.""" + lut, n, _ = labels_to_offsets(labels) + if n == 0: + return numpy.zeros(0), numpy.zeros(0) + padded = numpy.pad(numpy.ascontiguousarray(labels), 1) + hist = _xf_hist_kernel(padded, lut, n) + return hist @ _CROFTON_COEFS_4, hist @ _EULER_COEFS_8 diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py index 0199c49..1f481f9 100644 --- a/test/test_sizeshape_numba.py +++ b/test/test_sizeshape_numba.py @@ -205,3 +205,74 @@ def test_numba_convex_area_empty(): from cp_measure.core.numba._sizeshape import convex_area_2d assert convex_area_2d(numpy.zeros((20, 20), numpy.int32)).shape == (0,) + + +# --- perimeter / perimeter_crofton / euler_number --- + + +def _ring(size=40): + """An object with a hole (euler_number = 0) plus a solid one (euler_number = 1).""" + m = numpy.zeros((size, size), numpy.int32) + m[5:25, 5:25] = 1 + m[10:20, 10:20] = 0 # punch a hole -> object 1 has euler 0 + m[28:36, 28:36] = 2 + return m + + +def _assert_per_crofton_euler(masks): + from cp_measure.core.numba._sizeshape import crofton_euler_2d, perimeter_2d + + ref = skimage.measure.regionprops_table( + masks, properties=["perimeter", "perimeter_crofton", "euler_number"] + ) + per = perimeter_2d(masks) + cro, eul = crofton_euler_2d(masks) + numpy.testing.assert_allclose(per, ref["perimeter"], rtol=1e-9, atol=1e-9) + numpy.testing.assert_allclose(cro, ref["perimeter_crofton"], rtol=1e-9, atol=1e-9) + numpy.testing.assert_array_equal(eul, ref["euler_number"]) # integer -> exact + + +@requires_numba +def test_numba_perimeter_euler_multi(): + _assert_per_crofton_euler(_square_objects(256, 4)) + + +@requires_numba +def test_numba_perimeter_euler_holed(): + _assert_per_crofton_euler(_ring()) + + +@requires_numba +def test_numba_perimeter_euler_noncontiguous(): + masks = numpy.zeros((96, 96), numpy.int32) + masks[10:30, 10:30] = 1 + masks[40:60, 40:60] = 3 + masks[70:90, 70:90] = 7 + _assert_per_crofton_euler(masks) + + +@requires_numba +def test_numba_perimeter_euler_edge_touching(): + masks = numpy.zeros((64, 64), numpy.int32) + masks[0:20, 0:20] = 1 + masks[44:64, 44:64] = 2 + _assert_per_crofton_euler(masks) + + +@requires_numba +def test_numba_perimeter_euler_irregular(): + rng = numpy.random.default_rng(1) + masks = numpy.zeros((128, 128), numpy.int32) + yy, xx = numpy.mgrid[0:128, 0:128] + for lab, (cy, cx) in enumerate(rng.integers(20, 108, size=(6, 2)), 1): + masks[(yy - cy) ** 2 + (xx - cx) ** 2 < rng.integers(40, 160)] = lab + _assert_per_crofton_euler(masks) + + +@requires_numba +def test_numba_perimeter_euler_empty(): + from cp_measure.core.numba._sizeshape import crofton_euler_2d, perimeter_2d + + assert perimeter_2d(numpy.zeros((20, 20), numpy.int32)).shape == (0,) + cro, eul = crofton_euler_2d(numpy.zeros((20, 20), numpy.int32)) + assert cro.shape == (0,) and eul.shape == (0,) From 79d733646a459a87b588834671d7249c860e07e0 Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 22:51:18 +0200 Subject: [PATCH 4/8] =?UTF-8?q?feat(numba/sizeshape):=20phase=204=20?= =?UTF-8?q?=E2=80=94=20get=5Fsizeshape=20wrapper=20+=20dispatch=20(option?= =?UTF-8?q?=20B,=20~3x)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Ties the kernels into the full `get_sizeshape` feature dict and registers the numba sizeshape backend (previously a NO-GO). Option B: axes / eccentricity / orientation are derived from the kernel central moments (`primitives/_moments.axes_eccentricity_orientation`, bit-exact incl. the symmetric ±pi/4 fallback), so the regionprops call is fully moment-free (verified 0 ms einsum) — only cheap props (area/bbox/centroid/extent/area_filled/image) and the scipy Euclidean EDT radius loop stay on host. Sources: moments/inertia + axes/ecc/orientation -> moment kernel; area_convex/solidity -> hull kernel; perimeter/crofton/euler -> pattern kernels; radii -> scipy EDT (kept). `to_bzyx` 2D-only; 3D volumes fall back to the numpy baseline. Also fixes `bulk._numba_registries`, which the integration merge had left as un-parseable garbage (interleaved docstrings / em-dashes / duplicate returns) — rewritten cleanly to compose all numba backends (intensity, granularity, zernike, radial_zernikes, radial_distribution, texture, feret, sizeshape, coloc+costes) and registers `sizeshape`. Measured: get_sizeshape 305 -> 100 ms (~3x), all 78 features match the numpy backend (raw moments / area_convex / euler exact; the rest <=4e-7). Tests: end-to-end vs numpy get_sizeshape (multi / non-contiguous / calculate_advanced+new_features variants), 3D fallback, and set_accelerator("numba") dispatch routing. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/bulk.py | 96 ++++++---------- src/cp_measure/core/numba/__init__.py | 2 + src/cp_measure/core/numba/_sizeshape.py | 141 +++++++++++++++++++++++- src/cp_measure/primitives/_moments.py | 27 +++++ test/test_sizeshape_numba.py | 86 +++++++++++++++ 5 files changed, 286 insertions(+), 66 deletions(-) diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py index 4587a5c..5f9ab4e 100644 --- a/src/cp_measure/bulk.py +++ b/src/cp_measure/bulk.py @@ -48,81 +48,49 @@ def _numba_registries() -> dict[str, dict[str, Callable]]: """Registries for the 'numba' accelerator. - Composes the numba implementations (``intensity``, ``granularity``) with the - Composes the numba implementations (``intensity``, ``texture``) with the - numpy implementations of every other feature — a single global "numba" - Composes the numba implementations (``intensity`` and the ``pearson`` / - ``manders_fold`` / ``rwc`` / ``costes`` / ``overlap`` colocalization features) - 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. - - Note: ``overlap`` is not in the numpy ``_CORRELATION`` registry, so the numba - correlation registry intentionally exposes one feature the numpy one does not - (the numba ``overlap`` backend exists and is cheap to surface). Adding - ``overlap`` to the numpy ``_CORRELATION`` for symmetry is a separate call. + Composes the numba implementations (intensity, granularity, zernike, radial_zernikes, + radial_distribution, texture, feret, sizeshape, and the colocalization features including + costes and overlap) with the numpy implementations of every other feature, so a single + global "numba" selection yields a full, accelerated-where-available feature set. This is + explicit per-function composition, NOT an error-driven fallback. ``overlap`` is numba-only + (not present in the numpy ``_CORRELATION`` registry). """ from cp_measure.core.numba import ( - get_granularity as _numba_granularity, - get_intensity as _numba_intensity, - Composes the numba implementations (``intensity``, ``zernike``, - ``radial_zernikes``, ``radial_distribution``) 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, - get_radial_distribution as _numba_radial_distribution, - get_radial_zernikes as _numba_radial_zernikes, - get_zernike as _numba_zernike, + get_correlation_costes, + get_correlation_manders_fold, + get_correlation_overlap, + get_correlation_pearson, + get_correlation_rwc, + get_feret, + get_granularity, + get_intensity, + get_radial_distribution, + get_radial_zernikes, + get_sizeshape, + get_texture, + get_zernike, ) return { "core": { **_CORE, - "intensity": _numba_intensity, - "granularity": _numba_granularity, - "zernike": _numba_zernike, - "radial_zernikes": _numba_radial_zernikes, - "radial_distribution": _numba_radial_distribution, + "intensity": get_intensity, + "granularity": get_granularity, + "zernike": get_zernike, + "radial_zernikes": get_radial_zernikes, + "radial_distribution": get_radial_distribution, + "texture": get_texture, + "feret": get_feret, + "sizeshape": get_sizeshape, }, - "correlation": _CORRELATION, - get_correlation_costes as _numba_costes, - get_correlation_manders_fold as _numba_manders_fold, - get_correlation_overlap as _numba_overlap, - get_correlation_pearson as _numba_pearson, - get_correlation_rwc as _numba_rwc, - Composes the numba implementations (``intensity``, ``feret``) 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_feret as _numba_feret, - get_intensity as _numba_intensity, - ) - - return { - "core": {**_CORE, "intensity": _numba_intensity}, "correlation": { **_CORRELATION, - "pearson": _numba_pearson, - "manders_fold": _numba_manders_fold, - "rwc": _numba_rwc, - "costes": _numba_costes, - "overlap": _numba_overlap, + "pearson": get_correlation_pearson, + "manders_fold": get_correlation_manders_fold, + "rwc": get_correlation_rwc, + "costes": get_correlation_costes, + "overlap": get_correlation_overlap, }, - get_intensity as _numba_intensity, - get_texture as _numba_texture, - ) - - return { - "core": {**_CORE, "intensity": _numba_intensity, "texture": _numba_texture}, - "core": {**_CORE, "intensity": _numba_intensity, "feret": _numba_feret}, - "correlation": _CORRELATION, } diff --git a/src/cp_measure/core/numba/__init__.py b/src/cp_measure/core/numba/__init__.py index e8c8249..48b545c 100644 --- a/src/cp_measure/core/numba/__init__.py +++ b/src/cp_measure/core/numba/__init__.py @@ -9,6 +9,7 @@ """ from cp_measure.core.numba._feret import get_feret +from cp_measure.core.numba._sizeshape import get_sizeshape from cp_measure.core.numba.measurecolocalization import ( get_correlation_costes, get_correlation_manders_fold, @@ -26,6 +27,7 @@ from cp_measure.core.numba.measuretexture import get_texture __all__ = [ + "get_sizeshape", "get_correlation_costes", "get_correlation_manders_fold", "get_correlation_overlap", diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py index abc89cc..9fea543 100644 --- a/src/cp_measure/core/numba/_sizeshape.py +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -8,13 +8,22 @@ matrices match to floating-point round-off. """ +import centrosome.cpmorphology import numba import numpy +import scipy.ndimage from numpy.typing import NDArray -from skimage.measure import grid_points_in_poly +from skimage.measure import grid_points_in_poly, regionprops_table -from cp_measure.primitives._moments import derive_normalized_hu +from cp_measure.core import measureobjectsizeshape as _ss +from cp_measure.primitives._moments import ( + axes_eccentricity_orientation, + derive_normalized_hu, + inertia_2d, +) from cp_measure.primitives.segment import labels_to_offsets +from cp_measure.primitives.shapes import to_bzyx +from cp_measure.utils import _ensure_np_scalar _ORDER = 4 @@ -443,3 +452,131 @@ def crofton_euler_2d( padded = numpy.pad(numpy.ascontiguousarray(labels), 1) hist = _xf_hist_kernel(padded, lut, n) return hist @ _CROFTON_COEFS_4, hist @ _EULER_COEFS_8 + + +# --- Full get_sizeshape wrapper ------------------------------------------------------------- +# Assembles the complete sizeshape feature dict, sourcing the einsum-heavy / QHull / per-region +# pieces from the numba kernels above and keeping only cheap, moment-free regionprops props +# (area / bbox / centroid / extent / area_filled / image) plus the scipy Euclidean EDT radius +# loop. With option B (axes/eccentricity/orientation derived from the central moments), +# regionprops computes no moments at all. + +# moment-free regionprops props (verified 0 ms einsum); `image` feeds the EDT radius loop. +_CHEAP_PROPS = [ + "image", + "area", + "area_bbox", + "equivalent_diameter_area", + "bbox", + "centroid", + "extent", +] + + +def _sizeshape_2d(labels, pixels, calculate_advanced, new_features): + props_list = _CHEAP_PROPS + (["area_filled"] if new_features else []) + props = regionprops_table(labels, pixels, properties=props_list) + area = props["area"] + + raw, central, normalized, hu = spatial_moments_2d(labels) + area_convex = convex_area_2d(labels) + perimeter = perimeter_2d(labels) + crofton, euler = crofton_euler_2d(labels) + axis_major, axis_minor, eccentricity, orientation = axes_eccentricity_orientation( + central + ) + + formfactor = 4.0 * numpy.pi * area / perimeter**2 + denom = numpy.maximum(4.0 * numpy.pi * area, 1.0) + compactness = perimeter**2 / denom + solidity = area / area_convex + + nobjects = len(props["image"]) + max_radius = numpy.zeros(nobjects) + mean_radius = numpy.zeros(nobjects) + median_radius = numpy.zeros(nobjects) + for index, mini_image in enumerate(props["image"]): + mini_image = numpy.pad(mini_image, 1) + distances = scipy.ndimage.distance_transform_edt(mini_image) + max_radius[index] = _ensure_np_scalar( + scipy.ndimage.maximum(distances, mini_image) + ) + mean_radius[index] = _ensure_np_scalar( + scipy.ndimage.mean(distances, mini_image) + ) + median_radius[index] = _ensure_np_scalar( + centrosome.cpmorphology.median_of_labels( + distances, mini_image.astype("int"), [1] + ) + ) + + results = { + _ss.F_AREA: area, + _ss.F_BBOX_AREA: props["area_bbox"], + _ss.F_CONVEX_AREA: area_convex, + _ss.F_EQUIVALENT_DIAMETER: props["equivalent_diameter_area"], + _ss.F_PERIMETER: perimeter, + _ss.F_MAJOR_AXIS_LENGTH: axis_major, + _ss.F_MINOR_AXIS_LENGTH: axis_minor, + _ss.F_ECCENTRICITY: eccentricity, + _ss.F_ORIENTATION: orientation * (180 / numpy.pi), + _ss.F_CENTER_X: props["centroid-1"], + _ss.F_CENTER_Y: props["centroid-0"], + _ss.F_MIN_X: props["bbox-1"], + _ss.F_MAX_X: props["bbox-3"], + _ss.F_MIN_Y: props["bbox-0"], + _ss.F_MAX_Y: props["bbox-2"], + _ss.F_FORM_FACTOR: formfactor, + _ss.F_EXTENT: props["extent"], + _ss.F_SOLIDITY: solidity, + _ss.F_COMPACTNESS: compactness, + _ss.F_EULER_NUMBER: euler, + _ss.F_MAXIMUM_RADIUS: max_radius, + _ss.F_MEAN_RADIUS: mean_radius, + _ss.F_MEDIAN_RADIUS: median_radius, + } + if new_features: + results |= {_ss.F_FILLED_AREA: props["area_filled"]} + + if calculate_advanced: + it_00, it_off, it_11, eig_0, eig_1 = inertia_2d(central) + for p in range(3): # spatial / central exposed for p in {0,1,2}, q in {0,1,2,3} + for q in range(4): + results[f"SpatialMoment_{p}_{q}"] = raw[:, p, q] + results[f"CentralMoment_{p}_{q}"] = central[:, p, q] + for p in range(4): # normalized full 4x4 + for q in range(4): + results[f"NormalizedMoment_{p}_{q}"] = normalized[:, p, q] + for k in range(7): + results[f"HuMoment_{k}"] = hu[:, k] + results |= { + _ss.F_INERTIA_TENSOR_0_0: it_00, + _ss.F_INERTIA_TENSOR_0_1: it_off, + _ss.F_INERTIA_TENSOR_1_0: it_off, + _ss.F_INERTIA_TENSOR_1_1: it_11, + _ss.F_INERTIA_TENSOR_EIGENVALUES_0: eig_0, + _ss.F_INERTIA_TENSOR_EIGENVALUES_1: eig_1, + } + + if new_features: + results |= {_ss.F_PERIMETER_CROFTON: crofton} + + return results + + +def get_sizeshape( + masks: NDArray[numpy.integer], + pixels: NDArray[numpy.floating] | None = None, + calculate_advanced: bool = True, + new_features: bool = True, + spacing=None, +) -> dict[str, NDArray[numpy.floating]]: + """numba sizeshape backend (2D). 3D volumes fall back to the numpy baseline.""" + masks_zyx, _pixels_zyx, unwrap = to_bzyx(masks, masks if pixels is None else pixels) + results = [ + _sizeshape_2d(m[0], None, calculate_advanced, new_features) + if m.shape[0] == 1 + else _ss.get_sizeshape(m, None, calculate_advanced, new_features, spacing) + for m in masks_zyx + ] + return unwrap(results) diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py index 6e75cdc..66cf6ed 100644 --- a/src/cp_measure/primitives/_moments.py +++ b/src/cp_measure/primitives/_moments.py @@ -104,6 +104,33 @@ def inertia_2d( return c, -b, a, half_trace + disc, half_trace - disc +def axes_eccentricity_orientation( + central: NDArray[numpy.floating], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """``(axis_major, axis_minor, eccentricity, orientation)`` from per-object central moments. + + Matches ``skimage.measure.regionprops`` to floating-point round-off, including the symmetric + fallback (``it00 == it11`` -> ±pi/4). Derived from the same inertia tensor / eigenvalues as + :func:`inertia_2d`, so requesting these no longer forces regionprops' moment einsum. + """ + it00, it_off, it11, eig_major, eig_minor = inertia_2d(central) + axis_major = 4 * numpy.sqrt(eig_major) + axis_minor = 4 * numpy.sqrt(eig_minor) + with numpy.errstate(invalid="ignore", divide="ignore"): + eccentricity = numpy.where(eig_major == 0, 0.0, numpy.sqrt(1 - eig_minor / eig_major)) + orientation = numpy.where( + it00 - it11 == 0, + numpy.where(it_off < 0, numpy.pi / 4, -numpy.pi / 4), + 0.5 * numpy.arctan2(-2 * it_off, it11 - it00), + ) + return axis_major, axis_minor, eccentricity, orientation + + def spatial_moments_2d( labels: NDArray[numpy.integer], ) -> tuple[ diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py index 1f481f9..df654b3 100644 --- a/test/test_sizeshape_numba.py +++ b/test/test_sizeshape_numba.py @@ -276,3 +276,89 @@ def test_numba_perimeter_euler_empty(): assert perimeter_2d(numpy.zeros((20, 20), numpy.int32)).shape == (0,) cro, eul = crofton_euler_2d(numpy.zeros((20, 20), numpy.int32)) assert cro.shape == (0,) and eul.shape == (0,) + + +# --- end-to-end get_sizeshape wrapper --- + + +def _assert_full_sizeshape_matches(masks, pixels, **kw): + import cp_measure.core.measureobjectsizeshape as numpy_ss + from cp_measure.core.numba._sizeshape import get_sizeshape as numba_ss + + got = numba_ss(masks, pixels, **kw) + ref = numpy_ss.get_sizeshape(masks, pixels, **kw) + assert set(got) == set(ref), set(got).symmetric_difference(ref) + for k in ref: + g, r = numpy.asarray(got[k], float), numpy.asarray(ref[k], float) + assert g.shape == r.shape, k + finite = numpy.abs(r[numpy.isfinite(r)]) # NormalizedMoment_0_* are always NaN + s = max(finite.max() if finite.size else 0.0, 1.0) + numpy.testing.assert_allclose( + g, r, rtol=1e-7, atol=1e-7 * s, equal_nan=True, err_msg=k + ) + + +@requires_numba +def test_get_sizeshape_matches_numpy_multi(): + masks = _square_objects(256, 4) + _assert_full_sizeshape_matches( + masks, _pixels := numpy.random.default_rng(0).random(masks.shape) + ) + + +@requires_numba +def test_get_sizeshape_matches_numpy_noncontiguous(): + masks = numpy.zeros((96, 96), numpy.int32) + masks[10:30, 10:30] = 1 + masks[40:60, 40:60] = 3 + masks[70:90, 70:90] = 7 + _assert_full_sizeshape_matches( + masks, numpy.random.default_rng(1).random(masks.shape) + ) + + +@requires_numba +def test_get_sizeshape_flag_variants(): + masks = _square_objects(160, 3) + pixels = numpy.random.default_rng(2).random(masks.shape) + _assert_full_sizeshape_matches(masks, pixels, calculate_advanced=False) + _assert_full_sizeshape_matches(masks, pixels, new_features=False) + _assert_full_sizeshape_matches( + masks, pixels, calculate_advanced=False, new_features=False + ) + + +@requires_numba +def test_get_sizeshape_3d_falls_back_to_numpy(): + import cp_measure.core.measureobjectsizeshape as numpy_ss + from cp_measure.core.numba._sizeshape import get_sizeshape as numba_ss + + rng = numpy.random.default_rng(3) + masks = numpy.zeros((12, 32, 32), numpy.int32) + zz, yy, xx = numpy.mgrid[0:12, 0:32, 0:32] + masks[(zz - 6) ** 2 + (yy - 16) ** 2 + (xx - 16) ** 2 < 60] = 1 + pixels = rng.random((12, 32, 32)) + got = numba_ss(masks, pixels) + ref = numpy_ss.get_sizeshape(masks, pixels) + assert set(got) == set(ref) + for k in ref: + numpy.testing.assert_allclose( + numpy.asarray(got[k], float), + numpy.asarray(ref[k], float), + rtol=1e-7, + atol=1e-6, + equal_nan=True, + ) + + +@requires_numba +def test_get_sizeshape_dispatch(): + import cp_measure + from cp_measure.bulk import _dispatch + from cp_measure.core.numba._sizeshape import get_sizeshape as numba_ss + + cp_measure.set_accelerator("numba") + try: + assert _dispatch("core")["sizeshape"] is numba_ss + finally: + cp_measure.set_accelerator(None) From a4974d064177a3b3141952068f8421dbc4758fde Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 23:20:42 +0200 Subject: [PATCH 5/8] perf(numba/sizeshape): share the labels_to_offsets pass across primitives MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The four sizeshape primitives (spatial_moments_2d, convex_area_2d, perimeter_2d, crofton_euler_2d) each recomputed labels_to_offsets over the full raster, so it ran 4x per get_sizeshape call. Add a _Prep NamedTuple wrapping the labels_to_offsets result (lut, n, offsets) — the only prep every primitive shares — plus a _foreground_prep() helper; the primitives now take an optional `prep` arg and _sizeshape_2d threads one prep into all four. Each primitive still computes its own prep when called standalone (prep=None), so the public functions and their tests are unchanged. The full foreground pixel list (rows/cols/obj nonzero) is needed only by the moment kernel, so it stays inside spatial_moments_2d rather than the shared prep — standalone perimeter_2d/crofton_euler_2d/convex_area_2d no longer pay for a full-raster nonzero they never read. convex_area_2d does its own nonzero over the boundary mask (sparse, not redundant). Removes 3 redundant labels_to_offsets passes: ~9 ms / ~5% on the 1080^2 / 132-object tile. 24 tests unchanged (bit-exactness, dispatch, edge cases); ruff clean. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/core/numba/_sizeshape.py | 70 +++++++++++++++++++------ 1 file changed, 54 insertions(+), 16 deletions(-) diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py index 9fea543..f100632 100644 --- a/src/cp_measure/core/numba/_sizeshape.py +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -8,6 +8,8 @@ matrices match to floating-point round-off. """ +from typing import NamedTuple + import centrosome.cpmorphology import numba import numpy @@ -28,6 +30,25 @@ _ORDER = 4 +class _Prep(NamedTuple): + """The ``labels_to_offsets`` result, computed once and threaded into every primitive. + + All four sizeshape primitives otherwise recompute ``labels_to_offsets`` independently (4x over + the same raster); the wrapper computes it once and passes it in. ``lut``/``n``/``offsets`` are + the only prep shared by every primitive — the full foreground pixel list (rows/cols/object + index) is needed by the moment kernel alone, so it is built there from ``lut``. + """ + + lut: NDArray[numpy.int64] + n: int + offsets: NDArray[numpy.int64] + + +def _foreground_prep(labels: NDArray[numpy.integer]) -> _Prep: + """``label->offsets`` (``lut``, object count, CSR ``offsets``), shared by every primitive.""" + return _Prep(*labels_to_offsets(labels)) + + @numba.njit(cache=True) def _moment_kernel( rows: NDArray[numpy.int64], @@ -89,6 +110,7 @@ def _moment_kernel( def spatial_moments_2d( labels: NDArray[numpy.integer], + prep: _Prep | None = None, ) -> tuple[ NDArray[numpy.floating], NDArray[numpy.floating], @@ -99,16 +121,18 @@ def spatial_moments_2d( Drop-in for ``cp_measure.primitives._moments.spatial_moments_2d`` (same object order, same derivation), but computes the moment matrices with the fused numba kernel instead of 32 - ``numpy.bincount`` passes. + ``numpy.bincount`` passes. ``prep`` supplies the shared ``labels_to_offsets`` result when + called from the wrapper; left ``None`` it is computed here (standalone use). """ - lut, n, _offsets = labels_to_offsets(labels) - if n == 0: + if prep is None: + prep = _foreground_prep(labels) + if prep.n == 0: empty = numpy.zeros((0, _ORDER, _ORDER)) return empty, empty, empty, numpy.zeros((0, 7)) rows, cols = numpy.nonzero(labels) - obj = lut[labels[rows, cols]].astype(numpy.int64) + obj = prep.lut[labels[rows, cols]].astype(numpy.int64) raw, central = _moment_kernel( - rows.astype(numpy.int64), cols.astype(numpy.int64), obj, n + rows.astype(numpy.int64), cols.astype(numpy.int64), obj, prep.n ) normalized, hu = derive_normalized_hu(central) return raw, central, normalized, hu @@ -224,10 +248,15 @@ def _hull_kernel( return out_x[:cur], out_y[:cur], hull_offsets -def convex_area_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: +def convex_area_2d( + labels: NDArray[numpy.integer], prep: _Prep | None = None +) -> NDArray[numpy.floating]: """Per-object ``area_convex`` (pixel count inside the convex hull), ordered by ascending - label. Bit-exact vs ``skimage.measure.regionprops``.""" - lut, n, offsets = labels_to_offsets(labels) + label. Bit-exact vs ``skimage.measure.regionprops``. Reuses ``prep``'s ``lut``/``n``/ + ``offsets`` when supplied; the boundary scan below is its own (smaller) pass.""" + if prep is None: + prep = _foreground_prep(labels) + lut, n, offsets = prep.lut, prep.n, prep.offsets if n == 0: return numpy.zeros(0) bnd = _boundary_mask(labels) @@ -432,9 +461,13 @@ def _xf_hist_kernel( return hist -def perimeter_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: +def perimeter_2d( + labels: NDArray[numpy.integer], prep: _Prep | None = None +) -> NDArray[numpy.floating]: """Per-object 4-connectivity perimeter, bit-exact vs ``skimage.regionprops``.""" - lut, n, _ = labels_to_offsets(labels) + if prep is None: + prep = _foreground_prep(labels) + lut, n = prep.lut, prep.n if n == 0: return numpy.zeros(0) return _perimeter_kernel( @@ -443,10 +476,12 @@ def perimeter_2d(labels: NDArray[numpy.integer]) -> NDArray[numpy.floating]: def crofton_euler_2d( - labels: NDArray[numpy.integer], + labels: NDArray[numpy.integer], prep: _Prep | None = None ) -> tuple[NDArray[numpy.floating], NDArray[numpy.floating]]: """Per-object ``(perimeter_crofton, euler_number)`` from the shared 2x2-config histogram.""" - lut, n, _ = labels_to_offsets(labels) + if prep is None: + prep = _foreground_prep(labels) + lut, n = prep.lut, prep.n if n == 0: return numpy.zeros(0), numpy.zeros(0) padded = numpy.pad(numpy.ascontiguousarray(labels), 1) @@ -478,10 +513,13 @@ def _sizeshape_2d(labels, pixels, calculate_advanced, new_features): props = regionprops_table(labels, pixels, properties=props_list) area = props["area"] - raw, central, normalized, hu = spatial_moments_2d(labels) - area_convex = convex_area_2d(labels) - perimeter = perimeter_2d(labels) - crofton, euler = crofton_euler_2d(labels) + prep = _foreground_prep( + labels + ) # one labels_to_offsets shared by all four primitives + raw, central, normalized, hu = spatial_moments_2d(labels, prep) + area_convex = convex_area_2d(labels, prep) + perimeter = perimeter_2d(labels, prep) + crofton, euler = crofton_euler_2d(labels, prep) axis_major, axis_minor, eccentricity, orientation = axes_eccentricity_orientation( central ) From 04bfedadef2d115b24c33600631c8134ea7944cf Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Sat, 6 Jun 2026 23:57:30 +0200 Subject: [PATCH 6/8] refactor(numba/sizeshape): share moment-feature assembly + clarify wrapper 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) --- src/cp_measure/core/numba/_sizeshape.py | 31 +++++++--------- src/cp_measure/primitives/_moments.py | 47 ++++++++++++++++++++++++- test/test_sizeshape_numba.py | 9 +++++ 3 files changed, 67 insertions(+), 20 deletions(-) diff --git a/src/cp_measure/core/numba/_sizeshape.py b/src/cp_measure/core/numba/_sizeshape.py index f100632..74d59d9 100644 --- a/src/cp_measure/core/numba/_sizeshape.py +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -22,6 +22,7 @@ axes_eccentricity_orientation, derive_normalized_hu, inertia_2d, + moment_feature_dict, ) from cp_measure.primitives.segment import labels_to_offsets from cp_measure.primitives.shapes import to_bzyx @@ -578,23 +579,10 @@ def _sizeshape_2d(labels, pixels, calculate_advanced, new_features): if calculate_advanced: it_00, it_off, it_11, eig_0, eig_1 = inertia_2d(central) - for p in range(3): # spatial / central exposed for p in {0,1,2}, q in {0,1,2,3} - for q in range(4): - results[f"SpatialMoment_{p}_{q}"] = raw[:, p, q] - results[f"CentralMoment_{p}_{q}"] = central[:, p, q] - for p in range(4): # normalized full 4x4 - for q in range(4): - results[f"NormalizedMoment_{p}_{q}"] = normalized[:, p, q] - for k in range(7): - results[f"HuMoment_{k}"] = hu[:, k] - results |= { - _ss.F_INERTIA_TENSOR_0_0: it_00, - _ss.F_INERTIA_TENSOR_0_1: it_off, - _ss.F_INERTIA_TENSOR_1_0: it_off, - _ss.F_INERTIA_TENSOR_1_1: it_11, - _ss.F_INERTIA_TENSOR_EIGENVALUES_0: eig_0, - _ss.F_INERTIA_TENSOR_EIGENVALUES_1: eig_1, - } + # off-diagonal is symmetric (it_0_1 == it_1_0); shared assembler owns the feature names + results |= moment_feature_dict( + raw, central, normalized, hu, (it_00, it_off, it_off, it_11, eig_0, eig_1) + ) if new_features: results |= {_ss.F_PERIMETER_CROFTON: crofton} @@ -609,8 +597,13 @@ def get_sizeshape( new_features: bool = True, spacing=None, ) -> dict[str, NDArray[numpy.floating]]: - """numba sizeshape backend (2D). 3D volumes fall back to the numpy baseline.""" - masks_zyx, _pixels_zyx, unwrap = to_bzyx(masks, masks if pixels is None else pixels) + """numba sizeshape backend (2D). 3D volumes fall back to the numpy baseline. + + ``pixels`` is accepted for dispatch-signature parity but unused: sizeshape is purely + geometric (no intensity-weighted features), exactly like the numpy backend. ``masks`` is + passed to ``to_bzyx`` only to normalise the batch/spatial axes. + """ + masks_zyx, _, unwrap = to_bzyx(masks, masks) results = [ _sizeshape_2d(m[0], None, calculate_advanced, new_features) if m.shape[0] == 1 diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py index 66cf6ed..1e3b6aa 100644 --- a/src/cp_measure/primitives/_moments.py +++ b/src/cp_measure/primitives/_moments.py @@ -122,7 +122,9 @@ def axes_eccentricity_orientation( axis_major = 4 * numpy.sqrt(eig_major) axis_minor = 4 * numpy.sqrt(eig_minor) with numpy.errstate(invalid="ignore", divide="ignore"): - eccentricity = numpy.where(eig_major == 0, 0.0, numpy.sqrt(1 - eig_minor / eig_major)) + eccentricity = numpy.where( + eig_major == 0, 0.0, numpy.sqrt(1 - eig_minor / eig_major) + ) orientation = numpy.where( it00 - it11 == 0, numpy.where(it_off < 0, numpy.pi / 4, -numpy.pi / 4), @@ -163,3 +165,46 @@ def spatial_moments_2d( normalized, hu = derive_normalized_hu(central) return raw, central, normalized, hu + + +def moment_feature_dict( + raw: NDArray[numpy.floating], + central: NDArray[numpy.floating], + normalized: NDArray[numpy.floating], + hu: NDArray[numpy.floating], + inertia: tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + ], +) -> dict[str, NDArray[numpy.floating]]: + """Assemble the ``calculate_advanced`` moment + inertia features of 2D ``get_sizeshape``. + + Single source of truth for these 53 feature names and the ``(p, q)`` orders exposed, shared by + the numpy and numba sizeshape backends (the key strings match the ``F_*`` constants in + ``core.measureobjectsizeshape``; the sizeshape golden test cross-checks the two). ``raw`` and + ``central`` are ``(n, 4, 4)`` with only ``p in {0, 1, 2}`` exposed; ``normalized`` is the full + ``(n, 4, 4)``; ``hu`` is ``(n, 7)``. ``inertia`` is the tuple + ``(it_0_0, it_0_1, it_1_0, it_1_1, eigenvalue_0, eigenvalue_1)``. + """ + it_0_0, it_0_1, it_1_0, it_1_1, eig_0, eig_1 = inertia + features: dict[str, NDArray[numpy.floating]] = {} + for p in range(3): # spatial / central expose p in {0,1,2}, q in {0,1,2,3} + for q in range(_ORDER): + features[f"SpatialMoment_{p}_{q}"] = raw[:, p, q] + features[f"CentralMoment_{p}_{q}"] = central[:, p, q] + for p in range(_ORDER): # normalized full 4x4 + for q in range(_ORDER): + features[f"NormalizedMoment_{p}_{q}"] = normalized[:, p, q] + for k in range(7): + features[f"HuMoment_{k}"] = hu[:, k] + features["InertiaTensor_0_0"] = it_0_0 + features["InertiaTensor_0_1"] = it_0_1 + features["InertiaTensor_1_0"] = it_1_0 + features["InertiaTensor_1_1"] = it_1_1 + features["InertiaTensorEigenvalues_0"] = eig_0 + features["InertiaTensorEigenvalues_1"] = eig_1 + return features diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py index df654b3..f5f3f3b 100644 --- a/test/test_sizeshape_numba.py +++ b/test/test_sizeshape_numba.py @@ -328,6 +328,15 @@ def test_get_sizeshape_flag_variants(): ) +@requires_numba +def test_get_sizeshape_empty(): + # all-background tile: exercises the full _sizeshape_2d dict assembly with zero objects + # (the individual kernels have empty tests, but the wrapper's assembly did not). + masks = numpy.zeros((24, 24), numpy.int32) + pixels = numpy.zeros((24, 24)) + _assert_full_sizeshape_matches(masks, pixels) + + @requires_numba def test_get_sizeshape_3d_falls_back_to_numpy(): import cp_measure.core.measureobjectsizeshape as numpy_ss From 3b77ad41f981a5da555fbaaa7b35eaaaea1d68ab Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 10 Jun 2026 01:07:19 +0200 Subject: [PATCH 7/8] fix(numba/sizeshape): clip inertia eigenvalues + group moment_feature_dict order MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Mirror the numpy-side fixes so the shared _moments.py converges: - inertia_2d clips eigenvalues to >=0 like skimage, so thin/oblique objects get axis_minor=0 / eccentricity=1 instead of NaN (axes_eccentricity_orientation reuses inertia_2d, so the fix propagates to the axis features too). - moment_feature_dict emits the moment keys in the grouped order of the PyPI 0.1.19 release (all Spatial, then Central, then Normalized, ...) instead of interleaving Spatial/Central by (p, q) — the numba get_sizeshape output column order now matches the numpy backend and the release. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cp_measure/primitives/_moments.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py index 1e3b6aa..28b44c5 100644 --- a/src/cp_measure/primitives/_moments.py +++ b/src/cp_measure/primitives/_moments.py @@ -101,7 +101,11 @@ def inertia_2d( c = central[:, 0, 2] / mu00 half_trace = (a + c) / 2 disc = numpy.sqrt(((c - a) / 2) ** 2 + b**2) - return c, -b, a, half_trace + disc, half_trace - disc + # Clip eigenvalues to >= 0 (skimage does the same): tiny-negative float error on degenerate / + # thin objects would otherwise give NaN axis lengths and eccentricity > 1. + eig_major = numpy.clip(half_trace + disc, 0.0, None) + eig_minor = numpy.clip(half_trace - disc, 0.0, None) + return c, -b, a, eig_major, eig_minor def axes_eccentricity_orientation( @@ -192,9 +196,13 @@ def moment_feature_dict( """ it_0_0, it_0_1, it_1_0, it_1_1, eig_0, eig_1 = inertia features: dict[str, NDArray[numpy.floating]] = {} + # Grouped order (all Spatial, then Central, then Normalized, ...) matching the CellProfiler / + # PyPI 0.1.19 release column order — NOT interleaved by (p, q). for p in range(3): # spatial / central expose p in {0,1,2}, q in {0,1,2,3} for q in range(_ORDER): features[f"SpatialMoment_{p}_{q}"] = raw[:, p, q] + for p in range(3): + for q in range(_ORDER): features[f"CentralMoment_{p}_{q}"] = central[:, p, q] for p in range(_ORDER): # normalized full 4x4 for q in range(_ORDER): From 289e9971a15697806f6c9b4ebfc17f96eb418eea Mon Sep 17 00:00:00 2001 From: Tim Treis Date: Wed, 10 Jun 2026 01:09:47 +0200 Subject: [PATCH 8/8] test(numba): fix backend-correctness syntax error + stale sizeshape assertion A missing ')' in test_set_accelerator_numba_composes_with_numpy made the whole test_backend_correctness.py file uncollectable, which masked a stale assertion: it expected sizeshape to stay on the numpy backend, but this branch adds a numba sizeshape backend, so under the numba accelerator it composes cp_measure.core.numba._sizeshape. Close the paren and assert the numba module. Co-Authored-By: Claude Opus 4.8 (1M context) --- test/test_backend_correctness.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_backend_correctness.py b/test/test_backend_correctness.py index b4e808e..34a8f16 100644 --- a/test/test_backend_correctness.py +++ b/test/test_backend_correctness.py @@ -113,6 +113,7 @@ def test_set_accelerator_numba_composes_with_numpy(): ) assert core["granularity"].__module__ == ( "cp_measure.core.numba.measuregranularity" + ) assert core["zernike"].__module__ == ( "cp_measure.core.numba.measureobjectsizeshape" ) @@ -130,8 +131,7 @@ def test_set_accelerator_numba_composes_with_numpy(): ), feature assert core["texture"].__module__ == "cp_measure.core.numba.measuretexture" assert core["feret"].__module__ == "cp_measure.core.numba._feret" - # Every other feature stays on the numpy backend. - assert core["sizeshape"].__module__ == "cp_measure.core.measureobjectsizeshape" + assert core["sizeshape"].__module__ == "cp_measure.core.numba._sizeshape" finally: cp_measure.set_accelerator(None)