fix(radial): per-object measurement, label-independent (closes #22)

src/cp_measure/core/measureobjectintensitydistribution.py

@@ -95,6 +95,7 @@ def get_radial_distribution(
     scaled: bool = True,
     bin_count: int = 4,
     maximum_radius: int = 100,
+    legacy: bool = False,
 ) -> dict[str, NDArray[numpy.floating]]:
     """
     Radial features (2D only)
@@ -132,6 +133,17 @@ def get_radial_distribution(
         creates the number of bins that you specify and creates equally spaced bin
         boundaries up to the maximum radius. Parts of the object that are beyond this
         radius will be counted in an overflow bin. The radius is measured in pixels.
+
+    legacy : bool
+        When False (default) each object is measured on its own cropped mask, so its
+        result is independent of every other label in the field (the Issue #22 fix).
+        This deliberately diverges from CellProfiler on multi-object images.
+
+        When True, reproduce the whole-image result of CellProfiler / centrosome,
+        where an object's radial geometry can be perturbed by neighbouring labels
+        (Issue #22, rooted upstream in scipy). Use this for CellProfiler-faithful
+        values; otherwise the default gives the corrected, per-object-independent
+        measurement.
     """
 
     if labels.ndim == 3:
@@ -140,6 +152,50 @@ def get_radial_distribution(
     if labels.dtype == bool:
         labels = labels.astype(numpy.integer)
 
+    if legacy:
+        return _radial_distribution_image(
+            labels, pixels, scaled, bin_count, maximum_radius
+        )
+
+    # Issue #22 fix: measure each object on its own cropped + 1px-padded mask, so the
+    # imported geometry (distance_to_edge / propagate) and the binning see ONLY this
+    # object — its result equals the baseline run on that object in isolation.
+    slices = scipy.ndimage.find_objects(labels)
+    present = [(label, sl) for label, sl in enumerate(slices, 1) if sl is not None]
+    if not present:  # no objects -> empty arrays with the right keys
+        return _radial_distribution_image(
+            labels, pixels, scaled, bin_count, maximum_radius
+        )
+    per_object = [
+        _radial_distribution_image(
+            numpy.pad(labels[sl] == label, 1).astype(int),
+            numpy.pad(pixels[sl], 1),
+            scaled,
+            bin_count,
+            maximum_radius,
+        )
+        for label, sl in present
+    ]
+    return {
+        key: numpy.concatenate([obj[key] for obj in per_object])
+        for key in per_object[0]
+    }
+
+
+def _radial_distribution_image(
+    labels: NDArray[numpy.integer],
+    pixels: NDArray[numpy.floating],
+    scaled: bool = True,
+    bin_count: int = 4,
+    maximum_radius: int = 100,
+) -> dict[str, NDArray[numpy.floating]]:
+    """Whole-image radial distribution over every label (the original algorithm).
+
+    All objects are measured together, so an object's geometry can be influenced by
+    other labels in the field. :func:`get_radial_distribution` calls this once per
+    isolated object crop (the Issue #22 fix) or once on the whole image
+    (``legacy=True``).
+    """
     unique_labels = numpy.unique(labels)
     unique_labels = unique_labels[unique_labels > 0]
     nobjects = len(unique_labels)

test/test_radial_distribution.py

@@ -0,0 +1,57 @@
+"""Issue #22: get_radial_distribution must give per-object results independent of
+other labels. The default (legacy=False) measures each object on its own crop;
+legacy=True keeps the original whole-image behaviour."""
+
+import numpy as np
+
+from cp_measure.core.measureobjectintensitydistribution import get_radial_distribution
+
+
+def _three_touching():
+    """Three abutting rectangles — a layout where the whole-image geometry leaks."""
+    labels = np.zeros((48, 48), dtype=np.int32)
+    labels[5:40, 5:18] = 1
+    labels[5:40, 18:31] = 2
+    labels[5:40, 31:44] = 3
+    pixels = np.random.default_rng(1).random((48, 48))
+    return labels, pixels
+
+
+def test_radial_per_object_independent_of_neighbours():
+    """#22: an object's features are identical with or without other labels present.
+
+    Both sides use the per-object default, so this holds bit-exactly even for
+    symmetric objects whose centre is a plateau (the crop is the same array)."""
+    labels, pixels = _three_touching()
+    new = get_radial_distribution(labels, pixels)
+    for lbl in (1, 2, 3):
+        alone = np.where(labels == lbl, 1, 0).astype(np.int32)
+        iso = get_radial_distribution(alone, pixels)
+        for key in new:
+            np.testing.assert_allclose(
+                new[key][lbl - 1], iso[key][0], equal_nan=True, err_msg=key
+            )
+
+
+def test_radial_legacy_leaks_where_new_does_not():
+    """The whole-image (legacy) path is perturbed by neighbours on this layout while
+    the per-object default is not — they differ, proving #22 is real and fixed."""
+    labels, pixels = _three_touching()
+    new = get_radial_distribution(labels, pixels)
+    legacy = get_radial_distribution(labels, pixels, legacy=True)
+    assert any(not np.allclose(new[k], legacy[k], equal_nan=True) for k in new)
+
+
+def test_radial_unique_centre_object_new_equals_legacy():
+    """For an unambiguous (odd-sized) centre, the per-object crop reproduces the
+    whole-image result exactly — confirming the crop introduces no difference of its
+    own; the only systematic divergence is the centre tie-break on symmetric objects."""
+    labels = np.zeros((41, 41), dtype=np.int32)
+    labels[8:29, 8:29] = 1  # 21x21 odd square -> unique centre pixel
+    pixels = np.random.default_rng(2).random((41, 41))
+    new = get_radial_distribution(labels, pixels)
+    legacy = get_radial_distribution(labels, pixels, legacy=True)
+    for k in new:
+        np.testing.assert_allclose(
+            new[k], legacy[k], rtol=1e-6, atol=1e-8, equal_nan=True, err_msg=k
+        )