diff --git a/openhcs/processing/backends/cellprofiler/thresholding.py b/openhcs/processing/backends/cellprofiler/thresholding.py index 7cd50a34d1..295dbba3a5 100644 --- a/openhcs/processing/backends/cellprofiler/thresholding.py +++ b/openhcs/processing/backends/cellprofiler/thresholding.py @@ -2716,6 +2716,14 @@ def _li_threshold_float32_numpy(values: np.ndarray) -> float: def _li_tolerance_numpy(values: np.ndarray) -> float: if values.size < 2: return CELLPROFILER_LI_TOLERANCE + if values.size > 1024: + stride = (values.size + 1023) // 1024 + observed = np.unique(np.asarray(values.flat[::stride], dtype=np.float64)) + differences = np.diff(observed) + # This pair also occurs in the full array, whose minimum positive gap + # cannot be larger. A gap at most twice the floor proves the tolerance. + if np.any((differences > 0) & (differences <= 2.0 * CELLPROFILER_LI_TOLERANCE)): + return CELLPROFILER_LI_TOLERANCE unique_values = np.unique(np.asarray(values, dtype=np.float64).ravel()) if unique_values.size < 2: return CELLPROFILER_LI_TOLERANCE diff --git a/tests/unit/test_cellprofiler_threshold_hotpath.py b/tests/unit/test_cellprofiler_threshold_hotpath.py index 3ce92c8876..8e3423c8de 100644 --- a/tests/unit/test_cellprofiler_threshold_hotpath.py +++ b/tests/unit/test_cellprofiler_threshold_hotpath.py @@ -12,6 +12,8 @@ CellProfilerThresholdMethod, CellProfilerThresholdScope, NumbaNumpyThresholdDiagnosticsBackendStrategy, + CELLPROFILER_LI_TOLERANCE, + _li_tolerance_numpy, cellprofiler_get_adaptive_threshold, cellprofiler_get_global_threshold, threshold, @@ -31,6 +33,54 @@ UINT16_SCALE = int(np.iinfo(np.uint16).max) +@pytest.mark.parametrize("layout", ("contiguous", "strided", "transposed")) +@pytest.mark.parametrize("witness", (True, False)) +def test_li_tolerance_preserves_full_array_gap_for_sparse_and_dense_values( + layout: str, witness: bool +) -> None: + values = np.arange(4096, dtype=np.float64).reshape(64, 64) + if witness: + values.flat[4] = CELLPROFILER_LI_TOLERANCE + else: + # The close pair misses the sparse observation; the full owner is needed. + values.flat[1] = CELLPROFILER_LI_TOLERANCE + if layout == "strided": + backing = np.empty((64, 128), dtype=np.float64) + backing[:, ::2] = values + values = backing[:, ::2] + elif layout == "transposed": + values = values.T + differences = np.diff(np.unique(values.astype(np.float64).ravel())) + expected = max( + float(np.min(differences[differences > 0]) / 2), CELLPROFILER_LI_TOLERANCE + ) + assert _li_tolerance_numpy(values) == expected + + +@pytest.mark.parametrize( + "values", + ( + np.array([], dtype=np.float32), + np.ones(4096, dtype=np.float32), + np.tile([0.0, 1.0], 4096), + np.tile([np.nan, np.inf, -np.inf, 0.0, 1.0], 1024), + np.tile([np.nan, np.inf, -np.inf, 0.0, CELLPROFILER_LI_TOLERANCE], 1024), + ), +) +def test_li_tolerance_preserves_degenerate_and_nonfinite_inputs( + values: np.ndarray, +) -> None: + with np.errstate(invalid="ignore"): + differences = np.diff(np.unique(values.astype(np.float64).ravel())) + positive = differences[differences > 0] + expected = ( + max(float(np.min(positive) / 2), CELLPROFILER_LI_TOLERANCE) + if positive.size + else CELLPROFILER_LI_TOLERANCE + ) + assert _li_tolerance_numpy(values) == expected + + def _quantized_stack() -> tuple[np.ndarray, np.ndarray]: raw = ((np.arange(60 * 16 * 16, dtype=np.uint32) * 37 + 11) % 4096).astype( np.uint16