feat(sanitize): central non-contiguous label-ID sanitation

README.md

@@ -189,7 +189,7 @@ measurecolocalization.get_correlation_overlap
 
 ### Important notes
 
-- **Contiguous labels**: The input labels must be sequential (e.g., `[1,2,3]`, not `[1,3,4]`). You can use `skimage.segmentation.relabel_sequential` to ensure compliance.
+- **Labels**: Any positive integer labels work — non-contiguous IDs (e.g. `[1, 3, 4]`) are relabelled to `1..N` internally without modifying your array, and results are reported against your original IDs. `featurize` and the bulk `get_*` registries sanitize by default (`sanitize=False` to opt out); raw measurement functions assume contiguous `1..N`, so wrap them with `cp_measure._sanitize.sanitize` if you call them directly with gapped IDs.
 - **Fidelity**: If you need to match CellProfiler measurements 1:1, you must convert your image arrays to float values between 0 and 1. For instance, if you have an array of data type uint16, you must divide them all by 65535. This is important for radial distribution measurements. For the four intensity quantile features (`LowerQuartileIntensity`, `MedianIntensity`, `UpperQuartileIntensity`, `MADIntensity`) you additionally need `legacy=True` — see below.
 - **Speed**: The Granularity measurement is particularly slow (~80% of the compute time). Skip this one it if speed is of utmost importance.
 - **Legacy percentile convention**: `get_intensity` (numpy and numba backends), `get_core_measurements`, `get_core_measurements_3d`, and `make_featurizer_config` accept `legacy: bool = False`. The default uses `numpy.percentile` 'linear' (`(n-1)*q`) quartiles and the textbook `median(|x - median(x)|)` MAD. Pass `legacy=True` to reproduce the original cp_measure / CellProfiler behavior: `n*q` quartiles and the `(1/ndim)`-quantile MAD (which returns the 33rd percentile in 3D rather than the median). All other intensity features are identical either way.

src/cp_measure/_sanitize.py

@@ -0,0 +1,63 @@
+"""Relabel non-contiguous mask label IDs (gaps or arbitrary values) to contiguous
+``1..N`` without mutating the caller's array, applied at the entry points. See
+:func:`sanitize`.
+"""
+
+import functools
+import inspect
+from typing import Callable
+
+import numpy
+from numpy.typing import NDArray
+from skimage.segmentation import relabel_sequential
+
+# Argument names that hold the label image across the ``get_*`` functions.
+_MASK_PARAMS = ("masks", "labels", "mask")
+
+
+def sanitize_masks(masks: NDArray) -> tuple[NDArray, NDArray[numpy.int64]]:
+    """Relabel positive labels to contiguous ``1..N``.
+
+    Returns ``(clean, ids)`` where ``clean`` has labels ``1..N`` (the input is
+    returned unchanged when already contiguous) and ``ids[i]`` is the original
+    label of rank ``i + 1``. Never mutates the input. Raises ``ValueError`` for
+    non-integer (non-boolean) or negative labels.
+    """
+    if masks.dtype != bool and not numpy.issubdtype(masks.dtype, numpy.integer):
+        raise ValueError(f"labels must be an integer array, got dtype {masks.dtype!r}")
+    if masks.min(initial=0) < 0:
+        raise ValueError("labels must be non-negative")
+    mx = int(masks.max(initial=0))
+    if mx == 0:
+        return masks, numpy.empty(0, dtype=numpy.int64)
+    # bincount is faster than unique but needs a bounded range; fall back for huge labels.
+    if mx <= masks.size:
+        ids = numpy.flatnonzero(numpy.bincount(masks.ravel(), minlength=mx + 1))
+    else:
+        ids = numpy.unique(masks)
+    ids = ids[ids > 0].astype(numpy.int64)
+    if ids.size == mx:
+        # already 1..N: no copy (cast only a bool mask up to an integer dtype)
+        return (masks if masks.dtype != bool else masks.astype(numpy.intp)), ids
+    clean, _forward, _inverse = relabel_sequential(masks)
+    return clean, ids
+
+
+def sanitize(func: Callable) -> Callable:
+    """Wrap a ``get_*`` function to relabel its label argument (named in
+    :data:`_MASK_PARAMS`) to ``1..N`` before the call; functions with no such
+    argument are returned unchanged. Use this only to call a raw measurement
+    function directly with gapped IDs — the entry points already sanitize.
+    """
+    sig = inspect.signature(func)
+    param = next((name for name in _MASK_PARAMS if name in sig.parameters), None)
+    if param is None:
+        return func
+
+    @functools.wraps(func)
+    def wrapper(*args, **kwargs):
+        bound = sig.bind(*args, **kwargs)
+        bound.arguments[param], _ids = sanitize_masks(bound.arguments[param])
+        return func(*bound.args, **bound.kwargs)
+
+    return wrapper

src/cp_measure/bulk.py

@@ -6,6 +6,7 @@
 from functools import partial
 from typing import Callable
 
+from cp_measure._sanitize import sanitize as _sanitize
 from cp_measure.core import (
     measurecolocalization,
     measuregranularity,
@@ -105,25 +106,42 @@ def _apply_legacy(core: dict[str, Callable], legacy: bool) -> dict[str, Callable
     }
 
 
-def get_core_measurements(legacy: bool = False) -> dict[str, Callable]:
+def _maybe_sanitize(funcs: dict[str, Callable], sanitize: bool) -> dict[str, Callable]:
+    """Wrap each function to relabel arbitrary IDs to ``1..N``; no-op when False."""
+    if not sanitize:
+        return funcs
+    return {name: _sanitize(fn) for name, fn in funcs.items()}
+
+
+def get_core_measurements(
+    legacy: bool = False, sanitize: bool = True
+) -> dict[str, Callable]:
     """Core per-object measurement functions.
 
     ``legacy`` (default False) selects the original CellProfiler ``n*q`` percentile
     convention for the intensity quartile/MAD features instead of the default
     ``numpy.percentile`` ``(n-1)*q``; see
     :func:`cp_measure.core.measureobjectintensity.get_intensity`.
+
+    ``sanitize`` (default True) wraps each function so non-contiguous label IDs
+    (gaps or arbitrary values) are relabelled to ``1..N`` before the call. Set
+    False only when the caller has already guaranteed contiguous labels.
     """
-    return _apply_legacy(_dispatch("core"), legacy)
+    core = _apply_legacy(_dispatch("core"), legacy)
+    return _maybe_sanitize(core, sanitize)
 
 
-def get_core_measurements_3d(legacy: bool = False) -> dict[str, Callable]:
-    """Return only measurements that support 3D input (see ``legacy`` above)."""
+def get_core_measurements_3d(
+    legacy: bool = False, sanitize: bool = True
+) -> dict[str, Callable]:
+    """Return only measurements that support 3D input (see ``legacy``/``sanitize`` above)."""
     core = {k: _dispatch("core")[k] for k in _3D_FEATURES}
-    return _apply_legacy(core, legacy)
+    return _maybe_sanitize(_apply_legacy(core, legacy), sanitize)
 
 
-def get_correlation_measurements() -> dict[str, Callable]:
-    return _dispatch("correlation")
+def get_correlation_measurements(sanitize: bool = True) -> dict[str, Callable]:
+    """Correlation/colocalization measurement functions (see ``sanitize`` above)."""
+    return _maybe_sanitize(_dispatch("correlation"), sanitize)
 
 
 def get_multimask_measurements() -> dict[str, Callable]:

src/cp_measure/core/measurecolocalization.py

@@ -443,7 +443,7 @@ def linear_costes(
                 i -= i_step
         except ValueError:
             break
-    return thr_fi_c, thr_si_c
+    return float(thr_fi_c), float(thr_si_c)
 
 
 def bisection_costes(
@@ -517,7 +517,7 @@ def bisection_costes(
     thr_fi_c = (valid - 1) / scale_max
     thr_si_c = (a * thr_fi_c) + b
 
-    return thr_fi_c, thr_si_c
+    return float(thr_fi_c), float(thr_si_c)
 
 
 # MODIFIED: This reproduces the behaviour of the block at
@@ -549,6 +549,10 @@ def get_correlation_pearson(
     pixels_2: NDArray[numpy.floating],
     masks: NDArray[numpy.integer],
 ) -> dict[str, list[float]]:
+    """Per-object Pearson correlation and slope between two channels.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     corrs: list[float] = []
     slopes: list[float] = []
     for fi, si in _iter_label_pixels(pixels_1, pixels_2, masks):
@@ -568,6 +572,10 @@ def get_correlation_manders_fold(
     masks: NDArray[numpy.integer],
     thr: int = 15,
 ) -> dict[str, list[float]]:
+    """Per-object Manders fold coefficients between two channels.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     m1_list: list[float] = []
     m2_list: list[float] = []
     for fi, si in _iter_label_pixels(pixels_1, pixels_2, masks):
@@ -587,6 +595,10 @@ def get_correlation_rwc(
     masks: NDArray[numpy.integer],
     thr: int = 15,
 ) -> dict[str, list[float]]:
+    """Per-object rank-weighted colocalization coefficients between two channels.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     r1: list[float] = []
     r2: list[float] = []
     for fi, si in _iter_label_pixels(pixels_1, pixels_2, masks):
@@ -607,6 +619,10 @@ def get_correlation_costes(
     fast_costes: str = M_FASTER,
     thr: int = 15,
 ) -> dict[str, list[float]]:
+    """Per-object Costes colocalization coefficients between two channels.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     scale = infer_scale(pixels_1)
     c1: list[float] = []
     c2: list[float] = []
@@ -627,6 +643,10 @@ def get_correlation_overlap(
     masks: NDArray[numpy.integer],
     thr: int = 15,
 ) -> dict[str, list[float]]:
+    """Per-object overlap and k1/k2 colocalization coefficients between two channels.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     overlap_list: list[float] = []
     k1_list: list[float] = []
     k2_list: list[float] = []

src/cp_measure/core/measuregranularity.py

@@ -66,7 +66,10 @@ def get_granularity(
     element_size: int = 10,
     granular_spectrum_length: int = 16,
 ) -> dict[str, NDArray[numpy.floating]]:
-    """
+    """Per-object granularity spectrum features.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+
     1. (Outcommented) Subsample image
     2.  Remove background pixels using a greyscale tophat filter
     3.  Calculate granular spectrum (size distribution) for all masks
@@ -246,8 +249,6 @@ def get_granularity(
     # Per-object stats use original-scale labels so the cell boundaries are exact.
     # start_mean uses the raw (non-background-subtracted) original pixels, matching
     # CellProfiler's ObjectRecord initialisation which also uses im_pixel_data directly.
-    unique_labels = numpy.unique(orig_mask)
-    unique_labels = unique_labels[unique_labels > 0]
     range_ = numpy.arange(1, numpy.max(orig_mask) + 1)
 
     current_mean = fix(scipy.ndimage.mean(orig_pixels, orig_mask, range_))
@@ -278,7 +279,7 @@ def get_granularity(
 
         # Calculate the means for the objects
         gss = numpy.zeros((0,))
-        if unique_labels.any():
+        if range_.size:
             new_mean = fix(scipy.ndimage.mean(rec_orig, orig_mask, range_))
             gss = (current_mean - new_mean) * 100 / start_mean
             current_mean = new_mean  # update running mean for next iteration

src/cp_measure/core/measureobjectintensity.py

@@ -144,6 +144,8 @@ def get_intensity(
 ) -> dict[str, NDArray[numpy.floating]]:
     """Per-object intensity features.
 
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+
     Walks each object on its ``scipy.ndimage.find_objects`` bounding box rather
     than the full image; for each label the per-pixel reductions, quantiles, MAD,
     max position, and centroid run on the bbox crop only. 2D inputs are promoted

src/cp_measure/core/measureobjectintensitydistribution.py

@@ -99,6 +99,8 @@ def get_radial_distribution(
     """
     Radial features (2D only)
 
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+
     zernike_degree : int
         Maximum zernike moment.
 
@@ -138,11 +140,10 @@ def get_radial_distribution(
         return {}
 
     if labels.dtype == bool:
-        labels = labels.astype(numpy.integer)
+        labels = labels.astype(numpy.int64)
 
-    unique_labels = numpy.unique(labels)
-    unique_labels = unique_labels[unique_labels > 0]
-    nobjects = len(unique_labels)
+    nobjects = int(labels.max())
+    unique_labels = numpy.arange(1, nobjects + 1)
     d_to_edge = centrosome.cpmorphology.distance_to_edge(labels)
 
     # Find the point in each object farthest away from the edge.
@@ -282,7 +283,7 @@ def get_radial_distribution(
 
         mask = pixel_count == 0
 
-        radial_means: numpy.ma.MaskedArray = numpy.ma.masked_array(
+        radial_means: numpy.ma.MaskedArray = numpy.ma.MaskedArray(
             radial_values / pixel_count, mask
         )
 
@@ -310,13 +311,15 @@ def get_radial_zernikes(
     pixels: NDArray[numpy.floating],
     zernike_degree: int = 9,
 ) -> dict[str, NDArray[numpy.floating]]:
-    # Radial Zernike features (2D only)
+    """Per-object radial Zernike features (2D only).
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     if labels.ndim == 3:
         return {}
     zernike_indexes = centrosome.zernike.get_zernike_indexes(zernike_degree + 1)
 
-    unique_labels = numpy.unique(labels)  # Will be used later for scipy.ndimage.sum
-    unique_labels = unique_labels[unique_labels > 0]
+    unique_labels = numpy.arange(1, int(labels.max()) + 1)
     # MODIFIED: Delegate index generation to the minimum_enclosing_circle
     # MODIFIED: We assume non-overlapping labels for now
     # TODO Support label overlap (i.e., format in ijv)

src/cp_measure/core/measureobjectsizeshape.py

@@ -573,7 +573,10 @@ def get_sizeshape(
     new_features: bool = True,
     spacing: Optional[tuple[float, ...]] = None,
 ) -> dict[str, NDArray[numpy.floating]]:
-    """Compute the measurements for multiple object masks."""
+    """Compute the size/shape measurements for multiple object masks.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     # Properties available for both 2d and 3d
     desired_properties = [
         "image",
@@ -1008,9 +1011,10 @@ def get_zernike(
     pixels: NDArray[numpy.floating] | None,
     zernike_numbers: int = 9,
 ) -> dict[str, NDArray[numpy.floating]]:
-    #
-    # Zernike features (2D only)
-    #
+    """Per-object Zernike shape features (2D only).
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     if masks.ndim == 3:
         return {}
     unique_indices = numpy.unique(masks)
@@ -1030,7 +1034,10 @@ def get_zernike(
 def get_feret(
     masks: NDArray[numpy.integer], pixels: NDArray[numpy.floating] | None
 ) -> dict[str, NDArray[numpy.floating]]:
-    # Feret diameter (2D only)
+    """Per-object Feret diameter features (2D only).
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+    """
     if masks.ndim == 3:
         return {}
     ijv = masks_to_ijv(masks)

src/cp_measure/core/measuretexture.py

@@ -160,7 +160,10 @@ def get_texture(
     scale: int = 3,
     gray_levels: int = 256,
 ) -> dict[str, NDArray[numpy.floating]]:
-    """
+    """Per-object Haralick texture features.
+
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+
     Parameters
     ----------
     gray_levels : int, optional (default is 256)
@@ -191,9 +194,6 @@ def get_texture(
     CellProfiler 3 versions it was fixed at 256.  The minimum number of levels is
     2, the maximum is 256.
     """
-    unique_labels = numpy.unique(masks)
-    unique_labels = unique_labels[unique_labels > 0]
-
     # Modified to use the number of dimensions in pixels to determine the number of directions
     n_directions = 13 if pixels.ndim > 2 else 4
 
@@ -210,7 +210,7 @@ def get_texture(
         ).astype(numpy.uint8)
     props = skimage.measure.regionprops(masks, pixels)
 
-    features = numpy.empty((n_directions, 13, len(unique_labels)))
+    features = numpy.empty((n_directions, 13, len(props)))
 
     for index, prop in enumerate(props):
         label_data = prop["intensity_image"]