diff --git a/probeflow/processing/edge_detection.py b/probeflow/processing/edge_detection.py index 8008bf3..1d3685c 100644 --- a/probeflow/processing/edge_detection.py +++ b/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) diff --git a/probeflow/processing/inverse_fft.py b/probeflow/processing/inverse_fft.py index 9914f57..c0a84cf 100644 --- a/probeflow/processing/inverse_fft.py +++ b/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)) diff --git a/tests/test_edge_detection.py b/tests/test_edge_detection.py index 0bd8489..b148608 100644 --- a/tests/test_edge_detection.py +++ b/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):