Prove the shared Li tolerance floor without sorting all pixels

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

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