Address inverse-FFT / edge-detection review findings

probeflow/processing/edge_detection.py

@@ -40,7 +40,11 @@ class EdgeDetectionResult:
     * ``display_image``       — array suitable for direct display (NaN-preserving).
     * ``edge_mask``           — boolean array suitable for the active-mask layer.
     * ``gradient_magnitude``  — continuous gradient response (Sobel/Scharr).
-    * ``gradient_orientation``— gradient direction in radians, ``arctan2(gy, gx)``.
+    * ``gradient_orientation``— gradient direction in radians, ``arctan2(gy, gx)``,
+      in **image coordinates** (the row axis y increases *downward*), so a
+      positive angle turns from +x toward the bottom of the image — not the
+      math/physical y-up convention.  NaN marks pixels with no valid gradient
+      (NaN-contaminated or outside the ROI).
     """
 
     method: str
@@ -237,13 +241,21 @@ def gradient_filter(
         When True, threshold the gradient *magnitude* (always, regardless of
         *output*) at the *threshold* percentile to produce ``edge_mask``.
     roi_mask:
-        Optional boolean mask; gradients outside it are zeroed and excluded
-        from any threshold mask.
+        Optional boolean mask; gradients outside it are excluded (magnitude → 0,
+        orientation → NaN) and dropped from any threshold mask.
     pixel_size_x_nm, pixel_size_y_nm:
         Physical pixel spacings.  When given, the column/row derivatives are
         scaled to physical units (∂z/∂x, ∂z/∂y) before forming the magnitude
         and orientation, so anisotropic pixels yield physically correct values.
         Falls back to *pixel_size_nm* (isotropic) and then to 1 px.
+
+    Notes
+    -----
+    Pixels whose 3×3 operator footprint touches a non-finite value are treated
+    as invalid (their gradient is a mean-fill artefact, not a real edge):
+    magnitude → 0, orientation → NaN, display → NaN, and they are excluded from
+    the threshold mask.  ``gradient_orientation`` is in image coordinates (row
+    axis points down); see :class:`EdgeDetectionResult`.
     """
     from skimage import filters as _skf
 
@@ -274,6 +286,28 @@ def gradient_filter(
     magnitude = np.hypot(gx, gy)
     orientation = np.arctan2(gy, gx)
 
+    # Mean-filling non-finite pixels (see ``_prepare``) creates a step at the
+    # hole rim, so every pixel whose 3×3 operator footprint touches a NaN has a
+    # spurious gradient.  Treat that 1-px band as invalid, not just the NaN
+    # pixels themselves.
+    if nan_mask.any():
+        from scipy.ndimage import binary_dilation
+        invalid = binary_dilation(nan_mask, structure=np.ones((3, 3), dtype=bool))
+    else:
+        invalid = nan_mask  # all-False, correct shape
+
+    # Pixels to keep: finite-and-clean, restricted to the ROI when given.
+    keep = ~invalid if roi is None else (~invalid & roi)
+
+    # Null the gradient outside ``keep`` *before* normalisation/thresholding so
+    # the spurious rim never sets the normalisation peak or the percentile cut.
+    # Magnitude/components use 0 (= no gradient); orientation uses NaN because
+    # 0.0 is a valid direction (+x) and would be mistaken for genuine data.
+    gx = np.where(keep, gx, 0.0)
+    gy = np.where(keep, gy, 0.0)
+    magnitude = np.where(keep, magnitude, 0.0)
+    orientation = np.where(keep, orientation, np.nan)
+
     if output == "magnitude":
         chosen = magnitude
     elif output == "x":
@@ -283,37 +317,24 @@ def gradient_filter(
     else:  # orientation
         chosen = orientation
 
-    if roi is not None:
-        chosen = np.where(roi, chosen, 0.0)
-        magnitude = np.where(roi, magnitude, 0.0)
-        # Keep the returned orientation field ROI-bounded too, so downstream
-        # consumers never see out-of-ROI gradient directions.
-        orientation = np.where(roi, orientation, 0.0)
-
     if normalize and output != "orientation":
         peak = float(np.nanmax(np.abs(chosen))) if chosen.size else 0.0
         if peak > 0:
             chosen = chosen / peak
 
     display = chosen.astype(np.float64, copy=True)
-    if nan_mask.any():
-        display[nan_mask] = np.nan
+    display[invalid] = np.nan
 
     edge_mask = None
     if threshold_to_mask:
-        ref = magnitude.copy()
-        if roi is not None:
-            ref = ref[roi]
+        ref = magnitude[keep]
         finite_ref = ref[np.isfinite(ref)] if ref.size else ref
         if finite_ref.size:
             cut = float(np.percentile(finite_ref, float(threshold)))
             # Require a strictly positive gradient: when the percentile cut is
             # 0 (flat or sparse-step images) ``>= cut`` would mark the entire
             # zero-gradient background as an edge.
-            edge_mask = (magnitude >= cut) & (magnitude > 0.0)
-            edge_mask &= ~nan_mask
-            if roi is not None:
-                edge_mask &= roi
+            edge_mask = (magnitude >= cut) & (magnitude > 0.0) & keep
         else:
             edge_mask = np.zeros(a.shape, dtype=bool)
 

probeflow/processing/inverse_fft.py

@@ -25,7 +25,7 @@
 from __future__ import annotations
 
 import math
-from dataclasses import dataclass
+from dataclasses import dataclass, replace
 
 import numpy as np
 
@@ -80,7 +80,7 @@ class FourierStrokes:
     offsets from DC; ``radius`` is the shared brush radius in pixels.  The union
     of the stamped discs forms the region.
     """
-    stamps: tuple
+    stamps: tuple[tuple[float, float], ...]
     radius: float
 
 
@@ -180,7 +180,6 @@ def fourier_region_mask(
                     mask, _strokes_contribution(xx, yy, cx, cy, region, sign, soft_px))
             continue
         # Ellipse / rect: build the primary, then mirror it through DC.
-        from dataclasses import replace
         variants = [region]
         if conjugate:
             variants.append(replace(region, dx=-region.dx, dy=-region.dy))

tests/test_edge_detection.py

@@ -179,15 +179,37 @@ def test_x_and_y_outputs_differ_for_diagonal(self):
         # y-slope is twice the x-slope, so |gy| median should exceed |gx|.
         assert np.median(np.abs(gy)) > np.median(np.abs(gx))
 
-    def test_roi_bounds_orientation_field(self):
-        # ROI restriction must also bound the returned orientation field, not
-        # just display/magnitude/edge_mask.
+    def test_roi_bounds_orientation_field_with_nan_sentinel(self):
+        # ROI restriction bounds the returned orientation field, and uses NaN
+        # (not 0.0, a valid +x direction) outside the ROI.
         img = _step_image(edge_col=32)
         roi = np.zeros_like(img, dtype=bool)
         roi[:, :16] = True
         res = gradient_filter(img, output="magnitude", roi_mask=roi)
         assert res.gradient_orientation is not None
-        assert not np.any(res.gradient_orientation[:, 16:])
+        assert np.all(np.isnan(res.gradient_orientation[:, 16:]))
+        assert not np.any(res.gradient_magnitude[:, 16:])  # magnitude → 0 outside ROI
+
+    def test_nan_border_does_not_create_spurious_edges(self):
+        # Flat background with a distant bright block, so the global mean used
+        # to fill the NaN hole differs from the local (zero) surroundings and
+        # would step the rim. The 1-px contaminated band must be excluded.
+        img = np.zeros((40, 40), dtype=np.float64)
+        img[0:5, 0:5] = 100.0          # distant feature → nonzero global mean
+        img[18:22, 18:22] = np.nan     # hole in the flat (zero) region
+        res = gradient_filter(img, output="magnitude", normalize=False,
+                              threshold_to_mask=True, threshold=90)
+        # No edge pixels in the contaminated band around the hole.
+        assert not res.edge_mask[16:24, 16:24].any()
+        # Magnitude is nulled (not spuriously large) on the rim.
+        assert np.nanmax(res.gradient_magnitude[16:24, 16:24]) < 1e-9
+
+    def test_clean_image_unaffected_by_nan_guard(self):
+        # With no NaNs, behaviour is unchanged (full keep region).
+        res = gradient_filter(_step_image(edge_col=32), output="magnitude",
+                              threshold_to_mask=True, threshold=90)
+        assert res.edge_mask.any()
+        assert np.isfinite(res.gradient_orientation).all()
 
     def test_invalid_operator_and_output_raise(self):
         with pytest.raises(ValueError):