diff --git a/README.md b/README.md
index 470c25d..b4212f7 100644
--- a/README.md
+++ b/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.
diff --git a/src/cp_measure/_sanitize.py b/src/cp_measure/_sanitize.py
new file mode 100644
index 0000000..7706eb7
--- /dev/null
+++ b/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
diff --git a/src/cp_measure/bulk.py b/src/cp_measure/bulk.py
index 1456b83..c59266b 100644
--- a/src/cp_measure/bulk.py
+++ b/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]:
diff --git a/src/cp_measure/core/measurecolocalization.py b/src/cp_measure/core/measurecolocalization.py
index 9c97a5f..2bd435e 100644
--- a/src/cp_measure/core/measurecolocalization.py
+++ b/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] = []
diff --git a/src/cp_measure/core/measuregranularity.py b/src/cp_measure/core/measuregranularity.py
index d0cffee..1b26e26 100644
--- a/src/cp_measure/core/measuregranularity.py
+++ b/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
diff --git a/src/cp_measure/core/measureobjectintensity.py b/src/cp_measure/core/measureobjectintensity.py
index 3ca99dd..9913cd1 100644
--- a/src/cp_measure/core/measureobjectintensity.py
+++ b/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
diff --git a/src/cp_measure/core/measureobjectintensitydistribution.py b/src/cp_measure/core/measureobjectintensitydistribution.py
index f142d2a..4c781c4 100644
--- a/src/cp_measure/core/measureobjectintensitydistribution.py
+++ b/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)
diff --git a/src/cp_measure/core/measureobjectsizeshape.py b/src/cp_measure/core/measureobjectsizeshape.py
index 0e32b17..7179604 100644
--- a/src/cp_measure/core/measureobjectsizeshape.py
+++ b/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)
diff --git a/src/cp_measure/core/measuretexture.py b/src/cp_measure/core/measuretexture.py
index 659fe58..d54b37c 100644
--- a/src/cp_measure/core/measuretexture.py
+++ b/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"]
diff --git a/src/cp_measure/core/numba/measureobjectintensity.py b/src/cp_measure/core/numba/measureobjectintensity.py
index 57cf841..11fd08c 100644
--- a/src/cp_measure/core/numba/measureobjectintensity.py
+++ b/src/cp_measure/core/numba/measureobjectintensity.py
@@ -39,7 +39,6 @@
     STD_INTENSITY_EDGE,
     UPPER_QUARTILE_INTENSITY,
 )
-from cp_measure.primitives.segment import label_to_idx_lut
 from cp_measure.primitives._segment_numba import (
     flatten_numba,
     inner_boundary,
@@ -58,6 +57,8 @@ def get_intensity(
 ) -> dict[str, NDArray[numpy.floating]]:
     """masks is a labeled array where 0 are background.
 
+    Labels must be contiguous ``1..N`` (see :func:`cp_measure._sanitize.sanitize`).
+
     ``legacy`` mirrors the numpy backend: False (default) uses ``numpy.percentile``
     'linear' quartiles + textbook median MAD; True reproduces the original
     CellProfiler ``n*q`` quartiles + ``(1/ndim)``-quantile MAD.
@@ -72,7 +73,7 @@ def get_intensity(
     elif pixels.ndim == 3 and masks.ndim == 2:  # 3D image, 2D mask
         masks = masks.reshape(1, *masks.shape)
 
-    lut, nobjects = label_to_idx_lut(masks)
+    nobjects = int(masks.max())
 
     integrated_intensity = numpy.zeros(nobjects)
     mean_intensity = numpy.zeros(nobjects)
@@ -94,7 +95,6 @@ def get_intensity(
     values, seg0, xc, yc, zc = flatten_numba(
         numpy.ascontiguousarray(masks),
         numpy.ascontiguousarray(masked_image),
-        lut,
     )
     has_objects = values.size > 0
 
@@ -154,7 +154,7 @@ def get_intensity(
         else:
             emask = skimage.segmentation.find_boundaries(masks, mode="inner") > 0
         e_values = masked_image[emask].astype(numpy.float64)
-        e_seg0 = lut[masks[emask]]
+        e_seg0 = masks[emask] - 1
 
         if e_values.size > 0:
             ecount, esum, emin, emax = segment_stats(e_values, e_seg0, nobjects)
diff --git a/src/cp_measure/featurizer.py b/src/cp_measure/featurizer.py
index a1766a6..433770d 100644
--- a/src/cp_measure/featurizer.py
+++ b/src/cp_measure/featurizer.py
@@ -20,6 +20,8 @@
 
 import numpy as np
 
+from cp_measure._sanitize import sanitize_masks
+
 # Feature groups that only support 2D spatial data.
 _2D_ONLY = {"radial_distribution", "radial_zernikes", "zernike", "feret"}
 
@@ -184,9 +186,9 @@ def featurize(
     masks : numpy.ndarray
         Integer-labeled masks with shape ``(M, H, W)`` or
         ``(M, Z, H, W)``.  Must have the same ``ndim`` as *image*.
-        Background is 0; labels must be contiguous integers ``1..N``
-        (standard cp_measure convention, see
-        ``skimage.segmentation.relabel_sequential``).
+        Background is 0; any positive integer labels (non-contiguous IDs are
+        relabelled internally, array untouched, original IDs reported in the
+        rows; see :mod:`cp_measure._sanitize`).
     config : dict, optional
         Configuration dictionary produced by :func:`make_featurizer_config`.
         If ``None``, all features are enabled with default parameters.
@@ -217,12 +219,13 @@ def featurize(
 
     is_3d = image.ndim == 4
     legacy = config.get("legacy", False)
+    # Sanitize each mask once in the loop below, so fetch raw (unsanitized) funcs.
     core_funcs = (
-        get_core_measurements_3d(legacy=legacy)
+        get_core_measurements_3d(legacy=legacy, sanitize=False)
         if is_3d
-        else get_core_measurements(legacy=legacy)
+        else get_core_measurements(legacy=legacy, sanitize=False)
     )
-    corr_funcs = get_correlation_measurements()
+    corr_funcs = get_correlation_measurements(sanitize=False)
 
     if is_3d:
         _warn_2d_only_in_3d(config)
@@ -238,21 +241,20 @@ def featurize(
     columns: list[str] | None = None
 
     for mask_idx, object_name in enumerate(objects):
-        mask = masks[mask_idx]
-        # Assumes contiguous labels 1..max (standard cp_measure contract).
-        n_labels = int(mask.max())
-        if n_labels == 0:
+        # Relabel arbitrary IDs to 1..N once up front; keep originals for rows.
+        clean, ids = sanitize_masks(masks[mask_idx])
+        if ids.size == 0:
             continue
 
         results: dict[str, np.ndarray] = {}
 
         for func, params in shape_feats:
-            results.update(func(mask, dummy_pixels, **params))
+            results.update(func(clean, dummy_pixels, **params))
 
         for ch_idx, ch_name in enumerate(channels):
             pixels = image[ch_idx]
             for func, params in channel_feats:
-                for key, values in func(mask, pixels, **params).items():
+                for key, values in func(clean, pixels, **params).items():
                     results[f"{key}__{ch_name}"] = values
 
         n_ch = len(channels)
@@ -262,7 +264,7 @@ def featurize(
                 for key, values in func(
                     pixels_1=image[ch_i],
                     pixels_2=image[ch_j],
-                    masks=mask,
+                    masks=clean,
                     **params,
                 ).items():
                     results[f"{key}__{channels[ch_i]}__{channels[ch_j]}"] = values
@@ -284,9 +286,7 @@ def featurize(
         block = np.column_stack([results[c] for c in columns])
         all_blocks.append(block)
 
-        all_rows.extend(
-            (image_id, object_name, label) for label in range(1, n_labels + 1)
-        )
+        all_rows.extend((image_id, object_name, label) for label in ids)
 
     if not all_blocks:
         raise ValueError("all masks have no labels (all zeros)")
diff --git a/src/cp_measure/primitives/__init__.py b/src/cp_measure/primitives/__init__.py
index 3619525..79b534e 100644
--- a/src/cp_measure/primitives/__init__.py
+++ b/src/cp_measure/primitives/__init__.py
@@ -7,11 +7,11 @@
 segment index + per-axis coordinates) which the segment kernels reduce. All
 spatial structure (2D vs 3D) and any future batch/image axis are encoded in that
 flat segment index, so a single set of segment kernels covers every case without
-a rewrite. The numpy ``label_to_idx_lut`` (in ``segment``) builds the
-label→index lookup; the flattening + reductions themselves live in
-``_segment_numba``.
+a rewrite. Labels are the contiguous ``1..N`` cp_measure guarantees (see
+:mod:`cp_measure._sanitize`), so the segment index is ``label - 1``; the
+flattening + reductions live in ``_segment_numba``.
 
 This is an internal layer with no curated public API: import its building
-blocks directly from the concrete modules (``primitives.segment``,
-``primitives._segment_numba``) rather than re-exporting them here.
+blocks directly from the concrete modules (``primitives._segment_numba``)
+rather than re-exporting them here.
 """
diff --git a/src/cp_measure/primitives/_segment_numba.py b/src/cp_measure/primitives/_segment_numba.py
index a786d5f..b44249e 100644
--- a/src/cp_measure/primitives/_segment_numba.py
+++ b/src/cp_measure/primitives/_segment_numba.py
@@ -1,9 +1,9 @@
 """Numba segment kernels (single-threaded, cached).
 
 These are the numba implementation of the segment-reduce / segment-quantile
-primitives. They loop over the flat ``(values, seg0, coords)`` arrays produced
-by :mod:`cp_measure.primitives.segment` — no image shape, no batch axis — so one
-kernel set covers 2D, 3D, and (future) batched inputs unchanged.
+primitives. They loop over the flat ``(values, seg0, coords)`` arrays that
+:func:`flatten_numba` builds — no image shape, no batch axis — so one kernel set
+covers 2D, 3D, and (future) batched inputs unchanged.
 
 All kernels are ``@njit(cache=True)`` and serial: no ``prange``/``nogil``.
 Parallelism is the job of the (future) batch layer over images, not the kernel.
@@ -14,12 +14,14 @@
 
 
 @njit(cache=True)
-def flatten_numba(masks, pixels, lut):
+def flatten_numba(masks, pixels):
     """Flatten a labeled (Z, Y, X) image to ``(values, seg0, xc, yc, zc)``.
 
     Two grid scans (count, then fill) replace the numpy ``(masks>0)&isfinite``
     mask + ``numpy.nonzero`` + fancy-index gathers; coordinates are the loop
     indices. Background and non-finite pixels are dropped, in C (raster) order.
+    Labels are the contiguous ``1..N`` cp_measure guarantees (see
+    :mod:`cp_measure._sanitize`), so the segment index is ``label - 1``.
     ``masks`` and ``pixels`` must be C-contiguous; ``pixels`` may be any float
     dtype (kept values are upcast into the float64 ``values`` output).
     """
@@ -46,7 +48,7 @@ def flatten_numba(masks, pixels, lut):
                 if not np.isfinite(v):
                     continue
                 values[i] = v
-                seg0[i] = lut[L]
+                seg0[i] = L - 1
                 xc[i] = x
                 yc[i] = y
                 zc[i] = z
diff --git a/src/cp_measure/primitives/segment.py b/src/cp_measure/primitives/segment.py
deleted file mode 100644
index e5b5894..0000000
--- a/src/cp_measure/primitives/segment.py
+++ /dev/null
@@ -1,39 +0,0 @@
-"""Host-side segment helpers (numpy, backend-agnostic).
-
-A labeled image is reduced to a flat *segment* form — ``values[M]`` intensities,
-``seg0[M]`` 0-based segment indices, and ``xc/yc/zc[M]`` per-pixel coordinates —
-which the segment kernels then reduce. The flattening itself lives in the numba
-layer (:func:`cp_measure.primitives._segment_numba.flatten_numba`); this module
-holds only the numpy label→index lookup that feeds it.
-"""
-
-import numpy
-import scipy.ndimage
-from numpy.typing import NDArray
-
-
-def label_to_idx_lut(
-    masks: NDArray[numpy.integer],
-) -> tuple[NDArray[numpy.int64], int]:
-    """Build a ``label -> 0..n-1`` lookup over the sorted positive labels.
-
-    Returns ``(lut, n)`` where ``lut[label]`` is the segment index (and ``-1``
-    for absent labels / background) and ``n`` is the object count. Output arrays
-    are indexed by this segment rank, which equals ``label - 1`` when labels are
-    the contiguous ``1..n`` that cp_measure expects.
-
-    Present labels are enumerated with ``scipy.ndimage.find_objects`` (one pass)
-    rather than ``numpy.unique`` (a full-image sort): same ascending label set,
-    identical LUT, ~3-5x faster on large/3D masks.
-    """
-    if not numpy.issubdtype(masks.dtype, numpy.integer):
-        masks = masks.astype(numpy.intp, copy=False)
-    bboxes = scipy.ndimage.find_objects(masks)
-    labels = numpy.array(
-        [i + 1 for i, sl in enumerate(bboxes) if sl is not None], dtype=numpy.int64
-    )
-    n = int(labels.size)
-    max_label = int(labels[-1]) if n else 0
-    lut = numpy.full(max_label + 1, -1, dtype=numpy.int64)
-    lut[labels] = numpy.arange(n, dtype=numpy.int64)
-    return lut, n
diff --git a/test/test_sanitize.py b/test/test_sanitize.py
new file mode 100644
index 0000000..1c266bd
--- /dev/null
+++ b/test/test_sanitize.py
@@ -0,0 +1,100 @@
+"""Tests for central non-contiguous label sanitation (cp_measure._sanitize)."""
+
+import numpy as np
+import pytest
+
+from cp_measure._sanitize import sanitize, sanitize_masks
+from cp_measure.bulk import get_core_measurements
+from cp_measure.featurizer import featurize, make_featurizer_config
+
+
+def _three_objects(labels):
+    """3 squares in a 64x64 frame, labelled with the given (l1, l2, l3)."""
+    m = np.zeros((64, 64), dtype=np.int32)
+    m[5:19, 5:19] = labels[0]
+    m[5:19, 30:44] = labels[1]
+    m[35:49, 5:19] = labels[2]
+    return m
+
+
+_M3D = np.zeros((4, 16, 16), dtype=np.int32)
+_M3D[0, :4, :4] = 3
+_M3D[3, 10:14, 10:14] = 9
+
+
+@pytest.mark.parametrize(
+    "mask, clean, ids",
+    [
+        (_three_objects((1, 17, 5)), [1, 2, 3], [1, 5, 17]),  # gapped, non-monotonic
+        (_M3D, [1, 2], [3, 9]),  # 3D
+        (np.zeros((8, 8), np.int32), [], []),  # all background
+        (np.array([[False, True], [True, False]]), [1], [1]),  # bool -> one object
+    ],
+)
+def test_sanitize_relabels(mask, clean, ids):
+    out, out_ids = sanitize_masks(mask)
+    assert sorted(set(np.unique(out)) - {0}) == clean
+    assert out_ids.tolist() == ids
+
+
+def test_no_mutation_and_copy_policy():
+    contig, gapped = _three_objects((1, 2, 3)), _three_objects((1, 17, 5))
+    assert sanitize_masks(contig)[0] is contig  # already 1..N: no copy
+    before = gapped.copy()
+    assert sanitize_masks(gapped)[0] is not gapped  # relabelled: fresh copy
+    assert np.array_equal(gapped, before)  # input untouched
+
+
+@pytest.mark.parametrize("bad", [np.array([[0, -1, 2]]), np.array([[1.0, 2.0]])])
+def test_invalid_input_raises(bad):
+    with pytest.raises(ValueError):
+        sanitize_masks(bad)
+
+
+def test_entry_point_sanitizes_by_default():
+    # get_zernike assumes 1..N; the bulk entry point sanitizes by default, so
+    # gapped IDs yield the same values as their contiguous relabelling.
+    px = np.random.default_rng(0).random((64, 64))
+    zernike = get_core_measurements()["zernike"]
+    gapped = zernike(_three_objects((1, 17, 5)), px)
+    contig = zernike(_three_objects((1, 3, 2)), px)
+    for key in gapped:
+        np.testing.assert_allclose(gapped[key], contig[key], equal_nan=True)
+
+
+def test_raw_function_does_not_sanitize():
+    # sanitize=False returns the bare implementation: it assumes 1..N and breaks
+    # on gapped IDs (here zernike indexes past the relabelled range).
+    px = np.random.default_rng(0).random((64, 64))
+    raw = get_core_measurements(sanitize=False)["zernike"]
+    raw(_three_objects((1, 3, 2)), px)  # contiguous: fine
+    with pytest.raises(IndexError):
+        raw(_three_objects((1, 17, 5)), px)
+
+
+def test_sanitize_helper_wraps_raw_function():
+    # The user-facing escape hatch: wrap a raw function to handle gapped IDs.
+    px = np.random.default_rng(0).random((64, 64))
+    raw = get_core_measurements(sanitize=False)["zernike"]
+    wrapped = sanitize(raw)
+    gapped = wrapped(_three_objects((1, 17, 5)), px)
+    contig = raw(_three_objects((1, 3, 2)), px)
+    for key in gapped:
+        np.testing.assert_allclose(gapped[key], contig[key], equal_nan=True)
+
+
+def test_featurizer_uses_original_ids_and_sanitizes_once(monkeypatch):
+    import cp_measure.featurizer as fz
+
+    calls = []
+    real = fz.sanitize_masks
+    monkeypatch.setattr(fz, "sanitize_masks", lambda m: (calls.append(1), real(m))[1])
+    img = np.random.default_rng(0).random((2, 64, 64))
+    config = make_featurizer_config(["DNA", "ER"])
+
+    data_g, _, rows_g = featurize(img, _three_objects((1, 17, 5))[np.newaxis], config)
+    assert len(calls) == 1  # cost paid once per object-type mask
+    data_c, _, _ = featurize(img, _three_objects((1, 3, 2))[np.newaxis], config)
+    # original IDs in the rows, and identical geometry -> identical values
+    assert rows_g == [(None, "object", 1), (None, "object", 5), (None, "object", 17)]
+    np.testing.assert_allclose(data_g, data_c, equal_nan=True)