From 8db5df0a4661796bc0faa6ce2dd8a68d7e2ca961 Mon Sep 17 00:00:00 2001 From: anon Date: Fri, 3 Jul 2026 19:20:07 +0200 Subject: [PATCH 1/2] perf(sizeshape): batched separable moments + inertia, replacing regionprops einsum MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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) --- .github/scripts/run_benchmark.sh | 5 + src/cp_measure/core/measureobjectsizeshape.py | 116 ++------ src/cp_measure/primitives/_moments.py | 216 ++++++++++++++ test/test_sizeshape_moments.py | 281 ++++++++++++++++++ 4 files changed, 535 insertions(+), 83 deletions(-) create mode 100644 src/cp_measure/primitives/_moments.py create mode 100644 test/test_sizeshape_moments.py diff --git a/.github/scripts/run_benchmark.sh b/.github/scripts/run_benchmark.sh index fb101ae..c5c6dbb 100755 --- a/.github/scripts/run_benchmark.sh +++ b/.github/scripts/run_benchmark.sh @@ -4,6 +4,11 @@ # Usage: run_benchmark.sh set -euo pipefail +# Pin BLAS/OpenMP to one thread so timings reflect algorithmic cost, not incidental parallelism +# (cp_measure keeps core functions single-threaded; the batch layer is what parallelises). +export OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 +export NUMEXPR_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 + OUT="${1:-bench-out}" COMMIT="${2:-}" HEAD_DIR="$(pwd)" diff --git a/src/cp_measure/core/measureobjectsizeshape.py b/src/cp_measure/core/measureobjectsizeshape.py index 1adf67d..286398e 100644 --- a/src/cp_measure/core/measureobjectsizeshape.py +++ b/src/cp_measure/core/measureobjectsizeshape.py @@ -7,6 +7,12 @@ import scipy.ndimage import skimage.measure from centrosome import zernike +from cp_measure.primitives._moments import ( + axes_eccentricity_orientation, + inertia_2d, + moment_feature_dict, + spatial_moments_2d, +) from cp_measure.utils import _zernike_scores, masks_to_ijv __doc__ = """\ @@ -588,39 +594,25 @@ def get_sizeshape( "centroid", "euler_number", "extent", - "axis_major_length", - "axis_minor_length", ] # Features not in CellProfiler 4 if new_features: desired_properties += ["area_filled"] - # 2d specific properties if masks.ndim == 2: - desired_properties += [ - "eccentricity", - "orientation", - "perimeter", - "solidity", - ] + # 2D requests NOTHING moment-related from regionprops: the spatial / central / normalized + # / Hu moments, the inertia tensor, AND the axis lengths / eccentricity / orientation are + # all derived from the `spatial_moments_2d` central moments below (option B), so + # regionprops never runs its per-region moment einsum. + desired_properties += ["perimeter", "solidity"] if new_features: desired_properties += ["perimeter_crofton"] - - if calculate_advanced: - if masks.ndim == 2: - desired_properties += [ - "inertia_tensor", - "inertia_tensor_eigvals", - "moments", - "moments_hu", - "moments_central", - "moments_normalized", - ] - - else: + else: + # 3D: axis lengths still come from regionprops (3D moments are out of scope here). + desired_properties += ["axis_major_length", "axis_minor_length"] + if calculate_advanced: desired_properties += ["solidity"] - # These advanced props were not in CP4 for 3D images if new_features: desired_properties += [ @@ -632,13 +624,22 @@ def get_sizeshape( ] labels = masks - nobjects = (numpy.unique(masks) > 0).sum() + nobjects = int(masks.max()) # contiguous 1..N (see cp_measure._sanitize) results: dict[str, NDArray[numpy.floating]] = {} if labels.ndim == 2: props = skimage.measure.regionprops_table( labels, pixels, properties=desired_properties ) + # Option B: every moment-derived 2D feature (spatial / central / normalized / Hu moments, + # the inertia tensor, AND axis lengths / eccentricity / orientation) comes from this one + # batched-moment pass, so regionprops above ran no moment einsum. + raw, central, normalized, hu = spatial_moments_2d(labels) + inertia = inertia_2d(central) + axis_major, axis_minor, eccentricity, orientation = ( + axes_eccentricity_orientation(inertia) + ) + formfactor = 4.0 * numpy.pi * props["area"] / props["perimeter"] ** 2 denom = [max(x, 1) for x in 4.0 * numpy.pi * props["area"]] compactness = props["perimeter"] ** 2 / denom @@ -662,10 +663,10 @@ def get_sizeshape( F_CONVEX_AREA: props["area_convex"], F_EQUIVALENT_DIAMETER: props["equivalent_diameter_area"], F_PERIMETER: props["perimeter"], - F_MAJOR_AXIS_LENGTH: props["axis_major_length"], - F_MINOR_AXIS_LENGTH: props["axis_minor_length"], - F_ECCENTRICITY: props["eccentricity"], - F_ORIENTATION: props["orientation"] * (180 / numpy.pi), + F_MAJOR_AXIS_LENGTH: axis_major, + F_MINOR_AXIS_LENGTH: axis_minor, + F_ECCENTRICITY: eccentricity, + F_ORIENTATION: orientation * (180 / numpy.pi), F_CENTER_X: props["centroid-1"], F_CENTER_Y: props["centroid-0"], F_MIN_X: props["bbox-1"], @@ -689,61 +690,10 @@ def get_sizeshape( } if calculate_advanced: - results |= { - F_SPATIAL_MOMENT_0_0: props["moments-0-0"], - F_SPATIAL_MOMENT_0_1: props["moments-0-1"], - F_SPATIAL_MOMENT_0_2: props["moments-0-2"], - F_SPATIAL_MOMENT_0_3: props["moments-0-3"], - F_SPATIAL_MOMENT_1_0: props["moments-1-0"], - F_SPATIAL_MOMENT_1_1: props["moments-1-1"], - F_SPATIAL_MOMENT_1_2: props["moments-1-2"], - F_SPATIAL_MOMENT_1_3: props["moments-1-3"], - F_SPATIAL_MOMENT_2_0: props["moments-2-0"], - F_SPATIAL_MOMENT_2_1: props["moments-2-1"], - F_SPATIAL_MOMENT_2_2: props["moments-2-2"], - F_SPATIAL_MOMENT_2_3: props["moments-2-3"], - F_CENTRAL_MOMENT_0_0: props["moments_central-0-0"], - F_CENTRAL_MOMENT_0_1: props["moments_central-0-1"], - F_CENTRAL_MOMENT_0_2: props["moments_central-0-2"], - F_CENTRAL_MOMENT_0_3: props["moments_central-0-3"], - F_CENTRAL_MOMENT_1_0: props["moments_central-1-0"], - F_CENTRAL_MOMENT_1_1: props["moments_central-1-1"], - F_CENTRAL_MOMENT_1_2: props["moments_central-1-2"], - F_CENTRAL_MOMENT_1_3: props["moments_central-1-3"], - F_CENTRAL_MOMENT_2_0: props["moments_central-2-0"], - F_CENTRAL_MOMENT_2_1: props["moments_central-2-1"], - F_CENTRAL_MOMENT_2_2: props["moments_central-2-2"], - F_CENTRAL_MOMENT_2_3: props["moments_central-2-3"], - F_NORMALIZED_MOMENT_0_0: props["moments_normalized-0-0"], - F_NORMALIZED_MOMENT_0_1: props["moments_normalized-0-1"], - F_NORMALIZED_MOMENT_0_2: props["moments_normalized-0-2"], - F_NORMALIZED_MOMENT_0_3: props["moments_normalized-0-3"], - F_NORMALIZED_MOMENT_1_0: props["moments_normalized-1-0"], - F_NORMALIZED_MOMENT_1_1: props["moments_normalized-1-1"], - F_NORMALIZED_MOMENT_1_2: props["moments_normalized-1-2"], - F_NORMALIZED_MOMENT_1_3: props["moments_normalized-1-3"], - F_NORMALIZED_MOMENT_2_0: props["moments_normalized-2-0"], - F_NORMALIZED_MOMENT_2_1: props["moments_normalized-2-1"], - F_NORMALIZED_MOMENT_2_2: props["moments_normalized-2-2"], - F_NORMALIZED_MOMENT_2_3: props["moments_normalized-2-3"], - F_NORMALIZED_MOMENT_3_0: props["moments_normalized-3-0"], - F_NORMALIZED_MOMENT_3_1: props["moments_normalized-3-1"], - F_NORMALIZED_MOMENT_3_2: props["moments_normalized-3-2"], - F_NORMALIZED_MOMENT_3_3: props["moments_normalized-3-3"], - F_HU_MOMENT_0: props["moments_hu-0"], - F_HU_MOMENT_1: props["moments_hu-1"], - F_HU_MOMENT_2: props["moments_hu-2"], - F_HU_MOMENT_3: props["moments_hu-3"], - F_HU_MOMENT_4: props["moments_hu-4"], - F_HU_MOMENT_5: props["moments_hu-5"], - F_HU_MOMENT_6: props["moments_hu-6"], - F_INERTIA_TENSOR_0_0: props["inertia_tensor-0-0"], - F_INERTIA_TENSOR_0_1: props["inertia_tensor-0-1"], - F_INERTIA_TENSOR_1_0: props["inertia_tensor-1-0"], - F_INERTIA_TENSOR_1_1: props["inertia_tensor-1-1"], - F_INERTIA_TENSOR_EIGENVALUES_0: props["inertia_tensor_eigvals-0"], - F_INERTIA_TENSOR_EIGENVALUES_1: props["inertia_tensor_eigvals-1"], - } + # Spatial / central / normalized / Hu moments and the inertia tensor all come from the + # batched `central` + the single `inertia_2d` computed above — regionprops ran no + # moment einsum. moment_feature_dict owns the 53 feature names / orders. + results |= moment_feature_dict(raw, central, normalized, hu, inertia) if new_features: results |= { diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py new file mode 100644 index 0000000..11b0338 --- /dev/null +++ b/src/cp_measure/primitives/_moments.py @@ -0,0 +1,216 @@ +"""Per-object spatial-moment matrices via batched separable matrix products. + +``skimage.measure.regionprops_table`` computes ``moments`` / ``moments_central`` / +``moments_normalized`` / ``moments_hu`` per region with an ``einsum``-based routine whose +contraction path is re-derived for every object. Spatial moments are separable — for one object +``M = R.T @ image @ C`` where ``R`` / ``C`` hold the row / column coordinate powers — so the whole +set over all objects is two batched ``numpy.matmul`` calls over the stacked per-object bounding +boxes: no per-object Python loop, no ``einsum`` path re-derivation, and no full-image scatter. The +matmuls are tall-and-thin (inner size 4), so BLAS runs them single-threaded — parallelism stays +the batch layer's job, not this kernel's. + +Raw moments are bit-exact to ``regionprops``; the centroid-dependent matrices +(``central`` / ``normalized`` / ``hu``) match to floating-point round-off. Objects are ordered by +ascending label, exactly as ``regionprops_table``. The per-object boxes are padded to the largest +box, so a mask mixing one very large object with many tiny ones trades memory for the batching. +""" + +import numpy +import scipy.ndimage +from numpy.typing import NDArray + +# Moments up to order 3 (indices 0..3), matching skimage's order-3 regionprops matrices. +_ORDER = 4 + + +def _powers(x: NDArray[numpy.floating]) -> NDArray[numpy.floating]: + """Coordinate powers ``0..3`` stacked on a trailing axis (via multiply; ``x**k`` is ~20x slower).""" + x2 = x * x + return numpy.stack([numpy.ones_like(x), x, x2, x2 * x], axis=-1) + + +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 spatial_moments_2d( + labels: NDArray[numpy.integer], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """Per-object ``(raw, central, normalized, hu)`` spatial moments for a 2D label image. + + ``raw`` / ``central`` / ``normalized`` are ``(n, 4, 4)`` and ``hu`` is ``(n, 7)``, ordered by + ascending label. Drop-in for the matching ``regionprops_table`` columns + (``moments-p-q`` etc.); ``normalized`` is NaN where ``p + q < 2`` (skimage convention). + """ + n = int(labels.max()) + if n == 0: + empty = numpy.zeros((0, _ORDER, _ORDER)) + return empty, empty, empty, numpy.zeros((0, 7)) + + # Contiguous 1..N (see cp_measure._sanitize): find_objects index i is label i+1. Stack each + # object's bounding-box mask into (n, H, W), padded to the largest box; moments are taken in + # each object's local (bounding-box) frame, exactly as regionprops. + bboxes = scipy.ndimage.find_objects(labels) + height = max(sl[0].stop - sl[0].start for sl in bboxes) + width = max(sl[1].stop - sl[1].start for sl in bboxes) + masks = numpy.zeros((n, height, width)) + for i, sl in enumerate(bboxes): + box = labels[sl] == i + 1 + masks[i, : box.shape[0], : box.shape[1]] = box + + rows = numpy.arange(height, dtype=float) + cols = numpy.arange(width, dtype=float) + # Separable raw moments R.T @ (mask @ C), batched over objects as one BLAS matmul each. + raw = _powers(rows).T @ (masks @ _powers(cols)) + centre_r = raw[:, 1, 0] / raw[:, 0, 0] + centre_c = raw[:, 0, 1] / raw[:, 0, 0] + # Central moments by direct centred summation with per-object shifted grids (binomial-from-raw + # loses ~1e-4 to cancellation). + col_moments = masks @ _powers(cols[None, :] - centre_c[:, None]) + central = ( + _powers(rows[None, :] - centre_r[:, None]).transpose(0, 2, 1) @ col_moments + ) + + 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 raw, central, 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`` to + floating-point round-off. The tensor is ``[[c, -b], [-b, a]]`` with ``a = mu20/mu00``, + ``b = mu11/mu00``, ``c = mu02/mu00`` (skimage's row/col convention); eigenvalues are returned + in descending order. Reuses the central moments from :func:`spatial_moments_2d`, so the + inertia features need no separate regionprops einsum. + + Returns ``(t00, t_offdiag, t11, eig_major, eig_minor)`` — ``t_offdiag`` is both + off-diagonal entries (the tensor is symmetric). + """ + 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) + # 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( + inertia: tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + ], +) -> tuple[ + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], + NDArray[numpy.floating], +]: + """``(axis_major, axis_minor, eccentricity, orientation)`` from the per-object inertia tuple. + + Takes the output of :func:`inertia_2d` (so the eigendecomposition is computed once and shared + with the inertia features) and matches ``skimage.measure.regionprops`` to floating-point + round-off, including the symmetric fallback (``it00 == it11`` -> ±pi/4). CellProfiler reports + orientation in degrees; callers apply the ``180/pi`` conversion. + """ + it00, it_off, it11, eig_major, eig_minor = inertia + 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 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], + ], +) -> 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. Keys are emitted + in the *grouped* order of the CellProfiler / PyPI release (all Spatial, then Central, then + Normalized, then Hu, then the inertia tensor) so the output column order is unchanged. ``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 :func:`inertia_2d`'s + ``(it_0_0, it_off_diag, it_1_1, eig_0, eig_1)`` — the off-diagonal fills both + ``InertiaTensor_0_1`` and ``_1_0`` (symmetric tensor). + """ + it_00, it_off, it_11, 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] + 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): + features[f"NormalizedMoment_{p}_{q}"] = normalized[:, p, q] + for k in range(7): + features[f"HuMoment_{k}"] = hu[:, k] + features["InertiaTensor_0_0"] = it_00 + features["InertiaTensor_0_1"] = it_off + features["InertiaTensor_1_0"] = it_off + features["InertiaTensor_1_1"] = it_11 + features["InertiaTensorEigenvalues_0"] = eig_0 + features["InertiaTensorEigenvalues_1"] = eig_1 + return features diff --git a/test/test_sizeshape_moments.py b/test/test_sizeshape_moments.py new file mode 100644 index 0000000..acd91cf --- /dev/null +++ b/test/test_sizeshape_moments.py @@ -0,0 +1,281 @@ +"""Equivalence tests for the scatter-based spatial moments in ``get_sizeshape``. + +These lock **numerical fidelity to skimage regionprops**: the scatter ``spatial_moments_2d`` / +``inertia_2d`` must reproduce ``regionprops_table``'s ``moments`` / ``moments_central`` / +``moments_normalized`` / ``moments_hu`` / ``inertia_tensor`` columns to round-off (raw moments +bit-exact; centroid-dependent matrices ~1e-13 relative), and the derived axis lengths / +eccentricity / orientation must match regionprops. The rest of ``get_sizeshape`` is unchanged. +""" + +import numpy +import skimage.measure + +from cp_measure.core.measureobjectsizeshape import get_sizeshape +from cp_measure.primitives._moments import inertia_2d, spatial_moments_2d + +ATOL_REL = 1e-7 # relative tolerance; moments span many orders of magnitude + + +def _generate_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_moments_match(masks): + raw, central, normalized, hu = spatial_moments_2d(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): + scale = max( + numpy.nanmax(numpy.abs(ref[f"moments-{p}-{q}"])) if n else 0.0, 1.0 + ) + numpy.testing.assert_allclose( + raw[:, p, q], + ref[f"moments-{p}-{q}"], + rtol=ATOL_REL, + atol=ATOL_REL * scale, + ) + cscale = 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 * cscale, + ) + # normalized: NaN where p+q<2 in both + assert numpy.array_equal( + numpy.isnan(normalized[:, p, q]), + numpy.isnan(ref[f"moments_normalized-{p}-{q}"]), + ), f"NaN pattern mismatch normalized {p}{q}" + 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): + kscale = 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 * kscale + ) + + # inertia tensor + eigenvalues derived from the same central moments + it_00, it_off, it_11, eig_0, eig_1 = inertia_2d(central) + iref = 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_off, "inertia_tensor-1-0"), + (it_11, "inertia_tensor-1-1"), + (eig_0, "inertia_tensor_eigvals-0"), + (eig_1, "inertia_tensor_eigvals-1"), + ]: + scale = max(numpy.nanmax(numpy.abs(iref[key])) if n else 0.0, 1.0) + numpy.testing.assert_allclose( + got, iref[key], rtol=ATOL_REL, atol=ATOL_REL * scale + ) + + +def test_raw_moments_bit_exact(): + # raw spatial moments are integer-coordinate sums -> exactly equal to regionprops. + masks = _generate_square_objects(256, 4) + raw, *_ = spatial_moments_2d(masks) + ref = skimage.measure.regionprops_table(masks, properties=["moments"]) + for p in range(4): + for q in range(4): + numpy.testing.assert_array_equal(raw[:, p, q], ref[f"moments-{p}-{q}"]) + + +def test_moments_match_multi_object(): + _assert_moments_match(_generate_square_objects(256, 4)) + + +def test_moments_match_empty(): + # No objects -> (0, ...) arrays; the shared assert covers it against regionprops. + _assert_moments_match(numpy.zeros((20, 20), numpy.int32)) + + +def test_moments_match_edge_touching(): + masks = numpy.zeros((64, 64), numpy.int32) + masks[0:20, 0:20] = 1 + masks[44:64, 44:64] = 2 + _assert_moments_match(masks) + + +def test_moments_single_pixel_object(): + # degenerate: central/normalized are 0/NaN; must match regionprops' handling. + masks = numpy.zeros((32, 32), numpy.int32) + masks[16, 16] = 1 + masks[5:15, 5:15] = 2 + _assert_moments_match(masks) + + +def test_get_sizeshape_wires_moment_features(): + # get_sizeshape exposes spatial_moments_2d under the public F_* names; the helper-vs- + # regionprops accuracy is covered by the tests above, so this only checks the wiring. + masks = _generate_square_objects(200, 3) + pixels = numpy.random.default_rng(0).random(masks.shape) + out = get_sizeshape(masks, pixels) + raw, central, normalized, hu = spatial_moments_2d(masks) + 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): + numpy.testing.assert_array_equal( + out[f"SpatialMoment_{p}_{q}"], raw[:, p, q] + ) + numpy.testing.assert_array_equal( + out[f"CentralMoment_{p}_{q}"], central[:, p, q] + ) + for p in range(4): # normalized exposed for the full 4x4 + for q in range(4): + numpy.testing.assert_array_equal( + out[f"NormalizedMoment_{p}_{q}"], normalized[:, p, q] + ) + for k in range(7): + numpy.testing.assert_array_equal(out[f"HuMoment_{k}"], hu[:, k]) + it_00, it_off, it_11, eig_0, eig_1 = inertia_2d(central) + numpy.testing.assert_array_equal(out["InertiaTensor_0_0"], it_00) + numpy.testing.assert_array_equal(out["InertiaTensor_0_1"], it_off) + numpy.testing.assert_array_equal(out["InertiaTensor_1_0"], it_off) + numpy.testing.assert_array_equal(out["InertiaTensor_1_1"], it_11) + numpy.testing.assert_array_equal(out["InertiaTensorEigenvalues_0"], eig_0) + numpy.testing.assert_array_equal(out["InertiaTensorEigenvalues_1"], eig_1) + + +def _assert_axes_match(masks): + # Option B: axis lengths / eccentricity / orientation are derived from the scatter central + # moments instead of regionprops, so get_sizeshape requests no moments at all. They must still + # match regionprops to round-off (orientation reported in degrees: radians * 180/pi). + out = get_sizeshape(masks, masks.astype(float)) + ref = skimage.measure.regionprops_table( + masks, + properties=[ + "axis_major_length", + "axis_minor_length", + "eccentricity", + "orientation", + ], + ) + numpy.testing.assert_allclose( + out["MajorAxisLength"], ref["axis_major_length"], rtol=1e-9, atol=1e-9 + ) + numpy.testing.assert_allclose( + out["MinorAxisLength"], ref["axis_minor_length"], rtol=1e-9, atol=1e-9 + ) + numpy.testing.assert_allclose( + out["Eccentricity"], ref["eccentricity"], rtol=1e-9, atol=1e-9 + ) + numpy.testing.assert_allclose( + out["Orientation"], ref["orientation"] * (180 / numpy.pi), rtol=1e-9, atol=1e-9 + ) + + +def test_axes_match_single_object(): + masks = numpy.zeros((64, 64), numpy.int32) + masks[18:50, 20:45] = 1 # rectangle -> nonzero orientation, distinct axis lengths + _assert_axes_match(masks) + + +# Exact output column order of the PyPI 0.1.19 release (2D, new_features + calculate_advanced on). +# moment_feature_dict must keep this grouped order (all Spatial, then Central, then Normalized, +# then Hu, then the inertia tensor); a reorder is a silent schema break. +_RELEASE_KEY_ORDER = [ + "Area", + "BoundingBoxArea", + "ConvexArea", + "EquivalentDiameter", + "Perimeter", + "MajorAxisLength", + "MinorAxisLength", + "Eccentricity", + "Orientation", + "Center_X", + "Center_Y", + "BoundingBoxMinimum_X", + "BoundingBoxMaximum_X", + "BoundingBoxMinimum_Y", + "BoundingBoxMaximum_Y", + "FormFactor", + "Extent", + "Solidity", + "Compactness", + "EulerNumber", + "MaximumRadius", + "MeanRadius", + "MedianRadius", + "FilledArea", + *[f"SpatialMoment_{p}_{q}" for p in range(3) for q in range(4)], + *[f"CentralMoment_{p}_{q}" for p in range(3) for q in range(4)], + *[f"NormalizedMoment_{p}_{q}" for p in range(4) for q in range(4)], + *[f"HuMoment_{k}" for k in range(7)], + "InertiaTensor_0_0", + "InertiaTensor_0_1", + "InertiaTensor_1_0", + "InertiaTensor_1_1", + "InertiaTensorEigenvalues_0", + "InertiaTensorEigenvalues_1", + "PerimeterCrofton", +] + + +def test_sizeshape_key_order_matches_release(): + masks = numpy.zeros((40, 40), numpy.int32) + masks[5:25, 5:25] = 1 + assert list(get_sizeshape(masks, masks.astype(float))) == _RELEASE_KEY_ORDER + + +def test_axes_clip_thin_objects_match_skimage(): + # 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. + rng = numpy.random.default_rng(1) + for _ in range(50): + masks = numpy.zeros((40, 40), numpy.int32) + r0, c0 = rng.integers(2, 18, 2) + length = int(rng.integers(6, 18)) + dr, dc = rng.integers(-2, 3, 2) + for t in range(length): + r, c = r0 + t * dr, c0 + t * dc + if 0 <= r < 40 and 0 <= c < 40: + masks[r, c] = 1 + if masks.max() == 0: + continue + out = get_sizeshape(masks, masks.astype(float)) + assert not numpy.isnan(out["MinorAxisLength"]).any() + assert not numpy.isnan(out["Eccentricity"]).any() + ref = skimage.measure.regionprops_table( + masks, properties=["axis_minor_length", "eccentricity"] + ) + # atol=1e-6: a 1px line's minor axis is float noise (~4e-8) in skimage vs exactly 0 + # from our clip; both mean "zero width", so compare only above that floor. + numpy.testing.assert_allclose( + out["MinorAxisLength"], ref["axis_minor_length"], rtol=1e-7, atol=1e-6 + ) + numpy.testing.assert_allclose( + out["Eccentricity"], ref["eccentricity"], rtol=1e-7, atol=1e-6 + ) + + +def test_axes_match_single_pixel_and_line(): + # degenerate objects: single pixel (axes 0) and a 1-wide line (minor axis 0). + masks = numpy.zeros((40, 40), numpy.int32) + masks[20, 20] = 1 + masks[5:15, 8] = 2 + _assert_axes_match(masks) From 11e0251ead0b68fa09cc4841cbd96acd77b24589 Mon Sep 17 00:00:00 2001 From: anon Date: Tue, 7 Jul 2026 13:42:09 +0200 Subject: [PATCH 2/2] =?UTF-8?q?review(sizeshape):=20resolve=20#77=20commen?= =?UTF-8?q?ts=20=E2=80=94=20unified=20scaled=20tolerance=20+=20set-based?= =?UTF-8?q?=20schema=20test?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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) --- src/cp_measure/primitives/_moments.py | 2 +- test/test_sizeshape_moments.py | 111 +++++++++++++++----------- 2 files changed, 66 insertions(+), 47 deletions(-) diff --git a/src/cp_measure/primitives/_moments.py b/src/cp_measure/primitives/_moments.py index 11b0338..920fd6e 100644 --- a/src/cp_measure/primitives/_moments.py +++ b/src/cp_measure/primitives/_moments.py @@ -63,7 +63,7 @@ def spatial_moments_2d( ascending label. Drop-in for the matching ``regionprops_table`` columns (``moments-p-q`` etc.); ``normalized`` is NaN where ``p + q < 2`` (skimage convention). """ - n = int(labels.max()) + n = labels.max() # numpy int: fine for == 0 and as an array dimension if n == 0: empty = numpy.zeros((0, _ORDER, _ORDER)) return empty, empty, empty, numpy.zeros((0, 7)) diff --git a/test/test_sizeshape_moments.py b/test/test_sizeshape_moments.py index acd91cf..27f62d6 100644 --- a/test/test_sizeshape_moments.py +++ b/test/test_sizeshape_moments.py @@ -2,9 +2,9 @@ These lock **numerical fidelity to skimage regionprops**: the scatter ``spatial_moments_2d`` / ``inertia_2d`` must reproduce ``regionprops_table``'s ``moments`` / ``moments_central`` / -``moments_normalized`` / ``moments_hu`` / ``inertia_tensor`` columns to round-off (raw moments -bit-exact; centroid-dependent matrices ~1e-13 relative), and the derived axis lengths / -eccentricity / orientation must match regionprops. The rest of ``get_sizeshape`` is unchanged. +``moments_normalized`` / ``moments_hu`` / ``inertia_tensor`` columns, and the derived axis lengths / +eccentricity / orientation must match regionprops. Everything agrees to ~1 ULP once compared +against the right scale (see ``RTOL`` below). The rest of ``get_sizeshape`` is unchanged. """ import numpy @@ -13,7 +13,11 @@ from cp_measure.core.measureobjectsizeshape import get_sizeshape from cp_measure.primitives._moments import inertia_2d, spatial_moments_2d -ATOL_REL = 1e-7 # relative tolerance; moments span many orders of magnitude +# Relative equivalence tolerance. The batched matmul reorders sums vs regionprops' einsum, so +# every quantity matches to ~1 ULP *relative to its natural magnitude*; 1e-12 leaves a ~1000x +# margin. Central moments cancel toward 0, so they are compared against the raw-moment magnitude +# of the same order (the size of the summed terms) rather than their own near-zero value. +RTOL = 1e-12 def _generate_square_objects(size, n, gap_frac=0.7): @@ -28,35 +32,36 @@ def _generate_square_objects(size, n, gap_frac=0.7): return masks +def _generate_disk(size, cy, cx, r, label=1): + yy, xx = numpy.mgrid[0:size, 0:size] + masks = numpy.zeros((size, size), numpy.int32) + masks[(yy - cy) ** 2 + (xx - cx) ** 2 < r * r] = label + return masks + + def _assert_moments_match(masks): raw, central, normalized, hu = spatial_moments_2d(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): - scale = max( - numpy.nanmax(numpy.abs(ref[f"moments-{p}-{q}"])) if n else 0.0, 1.0 - ) + # Scale the absolute floor by the raw-moment magnitude: raw moments set the size of the + # terms summed for the central moments, so it bounds the cancellation round-off of the + # (structurally near-zero) central moments too. `initial=1.0` floors it and covers the + # empty-mask case. + scale = numpy.nanmax(numpy.abs(ref[f"moments-{p}-{q}"]), initial=1.0) numpy.testing.assert_allclose( - raw[:, p, q], - ref[f"moments-{p}-{q}"], - rtol=ATOL_REL, - atol=ATOL_REL * scale, - ) - cscale = max( - numpy.nanmax(numpy.abs(ref[f"moments_central-{p}-{q}"])) if n else 0.0, - 1.0, + raw[:, p, q], ref[f"moments-{p}-{q}"], rtol=RTOL, atol=RTOL * scale ) numpy.testing.assert_allclose( central[:, p, q], ref[f"moments_central-{p}-{q}"], - rtol=ATOL_REL, - atol=ATOL_REL * cscale, + rtol=RTOL, + atol=RTOL * scale, ) - # normalized: NaN where p+q<2 in both + # normalized moments are O(1); NaN where p+q<2 in both assert numpy.array_equal( numpy.isnan(normalized[:, p, q]), numpy.isnan(ref[f"moments_normalized-{p}-{q}"]), @@ -64,14 +69,13 @@ def _assert_moments_match(masks): numpy.testing.assert_allclose( normalized[:, p, q], ref[f"moments_normalized-{p}-{q}"], - rtol=ATOL_REL, - atol=ATOL_REL, + rtol=RTOL, + atol=RTOL, equal_nan=True, ) - for k in range(7): - kscale = max(numpy.nanmax(numpy.abs(ref[f"moments_hu-{k}"])) if n else 0.0, 1.0) + for k in range(7): # Hu invariants are O(1) numpy.testing.assert_allclose( - hu[:, k], ref[f"moments_hu-{k}"], rtol=ATOL_REL, atol=ATOL_REL * kscale + hu[:, k], ref[f"moments_hu-{k}"], rtol=RTOL, atol=RTOL ) # inertia tensor + eigenvalues derived from the same central moments @@ -87,10 +91,8 @@ def _assert_moments_match(masks): (eig_0, "inertia_tensor_eigvals-0"), (eig_1, "inertia_tensor_eigvals-1"), ]: - scale = max(numpy.nanmax(numpy.abs(iref[key])) if n else 0.0, 1.0) - numpy.testing.assert_allclose( - got, iref[key], rtol=ATOL_REL, atol=ATOL_REL * scale - ) + scale = numpy.nanmax(numpy.abs(iref[key]), initial=1.0) + numpy.testing.assert_allclose(got, iref[key], rtol=RTOL, atol=RTOL * scale) def test_raw_moments_bit_exact(): @@ -112,6 +114,14 @@ def test_moments_match_empty(): _assert_moments_match(numpy.zeros((20, 20), numpy.int32)) +def test_moments_match_disk(): + # Non-trivial fixture: a disk has an off-grid centroid, so its central moments are genuine + # near-total cancellations (many are structurally ~0). Axis-aligned rectangles hide this — + # their odd central moments are bit-exactly 0 — so this is the case that actually exercises + # the raw-magnitude-scaled tolerance in _assert_moments_match. + _assert_moments_match(_generate_disk(240, 119.3, 121.7, 80)) + + def test_moments_match_edge_touching(): masks = numpy.zeros((64, 64), numpy.int32) masks[0:20, 0:20] = 1 @@ -173,16 +183,16 @@ def _assert_axes_match(masks): ], ) numpy.testing.assert_allclose( - out["MajorAxisLength"], ref["axis_major_length"], rtol=1e-9, atol=1e-9 + out["MajorAxisLength"], ref["axis_major_length"], rtol=RTOL, atol=RTOL ) numpy.testing.assert_allclose( - out["MinorAxisLength"], ref["axis_minor_length"], rtol=1e-9, atol=1e-9 + out["MinorAxisLength"], ref["axis_minor_length"], rtol=RTOL, atol=RTOL ) numpy.testing.assert_allclose( - out["Eccentricity"], ref["eccentricity"], rtol=1e-9, atol=1e-9 + out["Eccentricity"], ref["eccentricity"], rtol=RTOL, atol=RTOL ) numpy.testing.assert_allclose( - out["Orientation"], ref["orientation"] * (180 / numpy.pi), rtol=1e-9, atol=1e-9 + out["Orientation"], ref["orientation"] * (180 / numpy.pi), rtol=RTOL, atol=RTOL ) @@ -192,10 +202,11 @@ def test_axes_match_single_object(): _assert_axes_match(masks) -# Exact output column order of the PyPI 0.1.19 release (2D, new_features + calculate_advanced on). -# moment_feature_dict must keep this grouped order (all Spatial, then Central, then Normalized, -# then Hu, then the inertia tensor); a reorder is a silent schema break. -_RELEASE_KEY_ORDER = [ +# Output features of the PyPI 0.1.19 release (2D, new_features + calculate_advanced on). Order is +# not part of the schema (downstream selects features by name), but the feature *set* is: a missing +# feature and an accidental extra (typo/dup) are both silent schema breaks. A genuine new feature is +# a deliberate one-line addition here. +_RELEASE_KEYS = [ "Area", "BoundingBoxArea", "ConvexArea", @@ -234,17 +245,24 @@ def test_axes_match_single_object(): ] -def test_sizeshape_key_order_matches_release(): +def test_sizeshape_feature_set_matches_release(): + # The emitted feature set must equal the release set (order ignored): catches both a dropped + # feature and an accidental extra. A genuine new feature means adding it to _RELEASE_KEYS. masks = numpy.zeros((40, 40), numpy.int32) masks[5:25, 5:25] = 1 - assert list(get_sizeshape(masks, masks.astype(float))) == _RELEASE_KEY_ORDER + keys = set(get_sizeshape(masks, masks.astype(float))) + expected = set(_RELEASE_KEYS) + assert keys == expected, { + "missing": sorted(expected - keys), + "unexpected": sorted(keys - expected), + } def test_axes_clip_thin_objects_match_skimage(): - # 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. + # A thin / 1px-wide object has almost no width, so the smaller eigenvalue of its inertia tensor + # is ~0; float round-off can push that value slightly below zero. inertia_2d clips it to 0 (as + # skimage does), 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. rng = numpy.random.default_rng(1) for _ in range(50): masks = numpy.zeros((40, 40), numpy.int32) @@ -263,13 +281,14 @@ def test_axes_clip_thin_objects_match_skimage(): ref = skimage.measure.regionprops_table( masks, properties=["axis_minor_length", "eccentricity"] ) - # atol=1e-6: a 1px line's minor axis is float noise (~4e-8) in skimage vs exactly 0 - # from our clip; both mean "zero width", so compare only above that floor. + # atol=1e-6 (looser than the global RTOL): a 1px line's minor axis is float noise (~1e-7) + # in skimage vs exactly 0 from our clip; both mean "zero width", so compare only above + # that floor. numpy.testing.assert_allclose( - out["MinorAxisLength"], ref["axis_minor_length"], rtol=1e-7, atol=1e-6 + out["MinorAxisLength"], ref["axis_minor_length"], rtol=RTOL, atol=1e-6 ) numpy.testing.assert_allclose( - out["Eccentricity"], ref["eccentricity"], rtol=1e-7, atol=1e-6 + out["Eccentricity"], ref["eccentricity"], rtol=RTOL, atol=1e-6 )