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 new file mode 100644 index 0000000..74d59d9 --- /dev/null +++ b/src/cp_measure/core/numba/_sizeshape.py @@ -0,0 +1,613 @@ +"""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. +""" + +from typing import NamedTuple + +import centrosome.cpmorphology +import numba +import numpy +import scipy.ndimage +from numpy.typing import NDArray +from skimage.measure import grid_points_in_poly, regionprops_table + +from cp_measure.core import measureobjectsizeshape as _ss +from cp_measure.primitives._moments import ( + 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 +from cp_measure.utils import _ensure_np_scalar + +_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], + 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], + prep: _Prep | None = None, +) -> 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. ``prep`` supplies the shared ``labels_to_offsets`` result when + called from the wrapper; left ``None`` it is computed here (standalone use). + """ + 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 = prep.lut[labels[rows, cols]].astype(numpy.int64) + raw, central = _moment_kernel( + rows.astype(numpy.int64), cols.astype(numpy.int64), obj, prep.n + ) + 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], 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``. 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) + 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 + + +# --- 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], prep: _Prep | None = None +) -> NDArray[numpy.floating]: + """Per-object 4-connectivity perimeter, bit-exact vs ``skimage.regionprops``.""" + if prep is None: + prep = _foreground_prep(labels) + lut, n = prep.lut, prep.n + 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], prep: _Prep | None = None +) -> tuple[NDArray[numpy.floating], NDArray[numpy.floating]]: + """Per-object ``(perimeter_crofton, euler_number)`` from the shared 2x2-config histogram.""" + 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) + 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"] + + 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 + ) + + 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) + # 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} + + 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. + + ``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 + 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 new file mode 100644 index 0000000..28b44c5 --- /dev/null +++ b/src/cp_measure/primitives/_moments.py @@ -0,0 +1,218 @@ +"""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) + # 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( + 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[ + 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 + + +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]] = {} + # 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): + 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_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) diff --git a/test/test_sizeshape_numba.py b/test/test_sizeshape_numba.py new file mode 100644 index 0000000..f5f3f3b --- /dev/null +++ b/test/test_sizeshape_numba.py @@ -0,0 +1,373 @@ +"""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) + + +# --- 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,) + + +# --- 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,) + + +# --- 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_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 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)