Merge pull request #440 from jhlegarreta/ref/avoid-divide-by-zero

REF: Avoid divide by zero warnings when denominator is zero

Author: jhlegarreta

src/nifreeze/data/filtering.py

@@ -35,6 +35,31 @@
 BVAL_ATOL = 100.0
 """b-value tolerance value."""
 
+CLIPPING_EMPTY_SELECTION_ERROR_MSG = """\
+Empty percentile selection after applying {constraints}finite filtering (p_min={p_min}, p_max={p_max}).
+"""
+"""Clipping empty selection error message."""
+
+CLIPPING_NONFINITE_THRESHOLDS_ERROR_MSG = """\
+Percentile thresholds are not finite (a_min={a_min}, a_max={a_max}, p_min={p_min}, p_max={p_max}, nonnegative={nonnegative}).
+"""
+"""Clipping non-finite thresholds error message."""
+
+CLIPPING_INVALID_THRESHOLDS_ERROR_MSG = """\
+Invalid percentile thresholds (a_max <= a_min) (a_min={a_min}, a_max={a_max}, p_min={p_min}, p_max={p_max}, nonnegative={nonnegative}).
+"""
+"""Clipping invalid percentile thresholds error message."""
+
+CLIPPING_DEGENERATE_RANGE_ERROR_MSG = """\
+Degenerate dynamic range after clipping/shift (den={den}, a_min={a_min}, a_max={a_max}, \
+nonnegative={nonnegative}, data_finite={data_finite}, unique~={unique}).
+"""
+"""Clipping degenerate dynamic range error message."""
+
+
+class ClippingValueError(RuntimeError):
+    """Raised when clipping cannot compute a valid clipping/scaling."""
+
 
 def advanced_clip(
     data: np.ndarray,
@@ -94,17 +119,55 @@ def advanced_clip(
     # Calculate stats on denoised version to avoid outlier bias
     denoised = median_filter(data, footprint=ball(3))
 
-    a_min = np.percentile(
-        np.asarray([denoised[denoised >= 0] if nonnegative else denoised]), p_min
-    )
-    a_max = np.percentile(
-        np.asarray([denoised[denoised >= 0] if nonnegative else denoised]), p_max
-    )
+    # Select values for percentile computation.
+    # If nonnegative=True, we must have at least one finite >=0 value.
+    sel = denoised[denoised >= 0] if nonnegative else denoised
+    sel = sel[np.isfinite(sel)]
+
+    if sel.size == 0:
+        constraints = "nonnegative constraint and " if nonnegative else ""
+        raise ClippingValueError(
+            CLIPPING_EMPTY_SELECTION_ERROR_MSG.format(
+                constraints=constraints, p_min=p_min, p_max=p_max
+            )
+        )
+
+    a_min = float(np.percentile(sel, p_min))
+    a_max = float(np.percentile(sel, p_max))
+
+    if not np.isfinite(a_min) or not np.isfinite(a_max):
+        raise ClippingValueError(
+            CLIPPING_NONFINITE_THRESHOLDS_ERROR_MSG.format(
+                a_min=a_min, a_max=a_max, p_min=p_min, p_max=p_max, nonnegative=nonnegative
+            )
+        )
+
+    if a_max <= a_min:
+        raise ClippingValueError(
+            CLIPPING_INVALID_THRESHOLDS_ERROR_MSG.format(
+                a_min=a_min, a_max=a_max, p_min=p_min, p_max=p_max, nonnegative=nonnegative
+            )
+        )
 
     # Clip and scale data
     np.clip(data, a_min=a_min, a_max=a_max, out=data)
     data -= data.min()
-    data /= data.max()
+    den = float(data.max())  # max-min because min is now 0
+    if not np.isfinite(den) or den == 0.0:
+        # Degenerate dynamic range after clipping
+        unique = len(np.unique(data)) if data.size < 1_000_000 else "too big"
+        raise ClippingValueError(
+            CLIPPING_DEGENERATE_RANGE_ERROR_MSG.format(
+                den=den,
+                a_min=a_min,
+                a_max=a_max,
+                nonnegative=nonnegative,
+                data_finite=bool(np.isfinite(data).all()),
+                unique=unique,
+            )
+        )
+
+    data /= den
 
     if invert:
         np.subtract(1.0, data, out=data)

test/test_filtering.py

@@ -30,13 +30,95 @@
 
 from nifreeze.data.filtering import (
     BVAL_ATOL,
+    CLIPPING_DEGENERATE_RANGE_ERROR_MSG,
+    CLIPPING_EMPTY_SELECTION_ERROR_MSG,
+    CLIPPING_INVALID_THRESHOLDS_ERROR_MSG,
+    CLIPPING_NONFINITE_THRESHOLDS_ERROR_MSG,
+    ClippingValueError,
     advanced_clip,
     dwi_select_shells,
     grand_mean_normalization,
     robust_minmax_normalization,
 )
 
 
+def test_advanced_clip_empty_selection():
+    data = -np.ones((5, 5, 5), dtype=np.float32)
+
+    with pytest.raises(ClippingValueError) as exc:
+        advanced_clip(data, inplace=False, nonnegative=True)
+
+    expected = CLIPPING_EMPTY_SELECTION_ERROR_MSG.format(
+        constraints="nonnegative constraint and ", p_min=35.0, p_max=99.98
+    )
+    assert str(exc.value) == expected
+
+
+def test_advanced_clip_nonfinite_thresholds(monkeypatch):
+    data = np.ones((5, 5, 5), dtype=np.float32)
+
+    # Force non-finite thresholds regardless of input data
+    def _percentile_returns_nan(*args, **kwargs):
+        return np.nan
+
+    from nifreeze.data import filtering
+
+    monkeypatch.setattr(filtering.np, "percentile", _percentile_returns_nan)
+
+    with pytest.raises(ClippingValueError) as exc:
+        advanced_clip(data, inplace=False, nonnegative=True, p_min=35.0, p_max=99.98)
+
+    expected = CLIPPING_NONFINITE_THRESHOLDS_ERROR_MSG.format(
+        a_min=float("nan"), a_max=float("nan"), p_min=35.0, p_max=99.98, nonnegative=True
+    )
+    assert str(exc.value) == expected
+
+
+def test_advanced_clip_invalid_thresholds():
+    data = np.arange(27, dtype=np.float32).reshape((3, 3, 3))
+    p = 50.0
+
+    with pytest.raises(ClippingValueError) as exc:
+        advanced_clip(data, inplace=False, nonnegative=True, p_min=p, p_max=p)
+
+    # With p_min == p_max, percentiles are equal (deterministically for this data)
+    sel = data[np.isfinite(data)]
+    a = float(np.percentile(sel, p))
+    expected = CLIPPING_INVALID_THRESHOLDS_ERROR_MSG.format(
+        a_min=a, a_max=a, p_min=p, p_max=p, nonnegative=True
+    )
+    assert str(exc.value) == expected
+
+
+def test_advanced_clip_degenerate_range(monkeypatch):
+    data = np.arange(27, dtype=np.float32).reshape((3, 3, 3))
+
+    # Ensure thresholds are valid (a_max > a_min) and deterministic.
+    # advanced_clip calls np.percentile twice (for a_min and a_max).
+    percentiles = iter([0.0, 1.0])
+
+    def _percentile(_arr, _q):
+        return next(percentiles)
+
+    monkeypatch.setattr(np, "percentile", _percentile)
+
+    # Force clipping to collapse all values to a constant after thresholds are deemed valid
+    def _clip(_data, a_min=None, a_max=None, out=None):
+        out.fill(0.5)
+        return out
+
+    monkeypatch.setattr(np, "clip", _clip)
+
+    with pytest.raises(ClippingValueError) as exc:
+        advanced_clip(data, inplace=False, nonnegative=True, p_min=35.0, p_max=99.98)
+
+    # After our fake clip: data is constant 0.5; subtract min: all zeros; den == 0.0
+    expected = CLIPPING_DEGENERATE_RANGE_ERROR_MSG.format(
+        den=0.0, a_min=0.0, a_max=1.0, nonnegative=True, data_finite=True, unique=1
+    )
+    assert str(exc.value) == expected
+
+
 @pytest.mark.random_uniform_ndim_data((32, 32, 32), 0.0, 2.0)
 @pytest.mark.parametrize(
     "p_min, p_max, nonnegative, dtype, invert, inplace",