FFT physics review: window-envelope compensation, odd-size soft-border centring

probeflow/processing/filters.py

@@ -159,7 +159,18 @@ def fourier_filter(
 
     F_filtered = F * mask
     F_filtered = np.fft.ifftshift(F_filtered)
-    result = np.fft.ifft2(F_filtered).real + mean_val
+    filtered = np.fft.ifft2(F_filtered).real
+    # The data window must be divided back out after the inverse transform,
+    # or its envelope stays baked into the filtered image: with the default
+    # Hann the output was vignetted to the mean at the borders (edge std
+    # 0.03 vs centre 0.77 on unit-variance noise — 2026-06-12 FFT review).
+    # fft_soft_border has always compensated its taper; same here, with the
+    # gain capped at 20× so near-zero edge weights cannot amplify residual
+    # spectral leakage into spikes (the interior, where the window is ≫0.05,
+    # is compensated exactly).
+    if window_canonical != "none":
+        filtered = filtered / np.maximum(win2d, 0.05)
+    result = filtered + mean_val
     result[nan_mask] = np.nan
     return result
 
@@ -375,9 +386,12 @@ def _tukey(n: int, frac: float) -> np.ndarray:
     tapered = centered * win2d
 
     F = np.fft.fftshift(np.fft.fft2(tapered))
-    cy, cx = Ny / 2.0, Nx / 2.0
-    yr = (np.arange(Ny) - cy) / max(cy, 1e-9)
-    xr = (np.arange(Nx) - cx) / max(cx, 1e-9)
+    # Radial mask in cycles/pixel scaled so cutoff=1.0 is Nyquist — the same
+    # convention (and the same DC bin) as fourier_filter. The previous
+    # index-arithmetic version ((arange - N/2.0) / (N/2.0)) is identical for
+    # even sizes but sat half a pixel off the fftshift DC bin for odd ones.
+    xr = np.fft.fftshift(np.fft.fftfreq(Nx)) / 0.5
+    yr = np.fft.fftshift(np.fft.fftfreq(Ny)) / 0.5
     Xr, Yr = np.meshgrid(xr, yr)
     R = np.sqrt(Xr ** 2 + Yr ** 2)
     if mode == "low_pass":

tests/test_fft_physics.py

@@ -380,3 +380,83 @@ def test_fourier_filter_radial_cutoff_isotropic_on_rectangular_image(self):
         assert max(rms_x, rms_y) / max(min(rms_x, rms_y), 1e-12) < 5.0, (
             f"Isotropic cutoff broken: rms_x={rms_x}, rms_y={rms_y}"
         )
+
+
+class TestWindowEnvelopeCompensation:
+    """2026-06-12 FFT review: fourier_filter windowed the data before the
+    transform but never divided the window back out, so the filtered image
+    carried the window envelope — with the GUI-default Hann, the output was
+    vignetted toward the mean at the borders (edge std 0.03 vs centre 0.77
+    on unit-variance noise). The pre-existing near-identity tests used
+    cutoff exactly 1.0/0.0, which early-return before any windowing and
+    could not see it."""
+
+    def test_near_identity_low_pass_is_spatially_flat(self):
+        rng = np.random.default_rng(0)
+        img = rng.normal(size=(64, 64)) + 5.0
+        out = fourier_filter(img, mode="low_pass", cutoff=0.99, window="hanning")
+        centre = float(np.std(out[24:40, 24:40]))
+        edge = float(np.std(out[:8, :8]))
+        assert edge > 0.5 * centre, (
+            f"window envelope baked into output: edge std {edge:.3f} vs "
+            f"centre {centre:.3f}"
+        )
+
+    def test_high_pass_response_uniform_across_image(self):
+        """A uniform-amplitude high-frequency stripe must keep roughly the
+        same amplitude at the border as in the middle after a high-pass."""
+        yy, xx = np.mgrid[:96, :96]
+        img = np.sin(2 * np.pi * 0.3 * xx)
+        out = fourier_filter(img, mode="high_pass", cutoff=0.2, window="hanning")
+        centre_amp = float(np.std(out[40:56, 40:56]))
+        edge_amp = float(np.std(out[4:12, 40:56]))
+        assert edge_amp > 0.5 * centre_amp
+
+    def test_window_none_behaviour_unchanged(self):
+        """No window → no compensation: pin equality with the raw masked
+        transform so the compensation cannot leak into the 'none' path."""
+        rng = np.random.default_rng(1)
+        img = rng.normal(size=(32, 32))
+        out = fourier_filter(img, mode="low_pass", cutoff=0.5, window="none")
+        mean = img.mean()
+        F = np.fft.fftshift(np.fft.fft2(img - mean))
+        qx = np.fft.fftshift(np.fft.fftfreq(32))
+        Qx, Qy = np.meshgrid(qx, qx)
+        mask = (np.sqrt(Qx**2 + Qy**2) / 0.5 <= 0.5).astype(float)
+        expected = np.fft.ifft2(np.fft.ifftshift(F * mask)).real + mean
+        np.testing.assert_allclose(out, expected, atol=1e-12)
+
+
+class TestSoftBorderOddSizes:
+    def test_dc_preserved_on_odd_sized_image(self):
+        """The radial mask must sit on the fftshift DC bin for odd sizes —
+        the old index-arithmetic centre was half a pixel off, so a tiny
+        low-pass could miss DC entirely and shift the image mean."""
+        rng = np.random.default_rng(2)
+        img = rng.normal(size=(63, 65)) + 7.0
+        out = fft_soft_border(img, mode="low_pass", cutoff=0.02,
+                              border_frac=0.1)
+        assert abs(float(np.nanmean(out)) - 7.0) < 0.05
+
+    def test_even_size_results_unchanged_by_mask_rewrite(self):
+        """For even sizes the fftfreq mask is identical to the old index
+        arithmetic — pin a checksum so the rewrite is provably a no-op."""
+        rng = np.random.default_rng(3)
+        img = rng.normal(size=(64, 64))
+        out = fft_soft_border(img, mode="low_pass", cutoff=0.4, border_frac=0.12)
+        # Old-mask reference computed inline.
+        a = img.astype(np.float64).copy()
+        mean = a.mean()
+        n = 64
+        edge = max(1, int(round(0.12 * n)))
+        ramp = 0.5 * (1.0 - np.cos(np.linspace(0.0, np.pi, edge)))
+        w = np.ones(n); w[:edge] = ramp; w[-edge:] = ramp[::-1]
+        win2d = np.outer(w, w)
+        F = np.fft.fftshift(np.fft.fft2((a - mean) * win2d))
+        cyx = n / 2.0
+        r1 = (np.arange(n) - cyx) / cyx
+        Xr, Yr = np.meshgrid(r1, r1)
+        mask = (np.sqrt(Xr**2 + Yr**2) <= 0.4).astype(float)
+        ref = np.fft.ifft2(np.fft.ifftshift(F * mask)).real
+        ref = ref / np.where(win2d > 1e-6, win2d, 1.0) + mean
+        np.testing.assert_allclose(out, ref, atol=1e-12)