From 310e11929b912228c566d002071faa7de735a626 Mon Sep 17 00:00:00 2001 From: Peter Jacobson Date: Fri, 12 Jun 2026 00:25:16 +1000 Subject: [PATCH] Compensate the data window in fourier_filter; centre soft-border mask for odd sizes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit FFT physics review pass (2026-06-12), pinned by new tests in tests/test_fft_physics.py: - fourier_filter windowed the data before the transform (GUI default: Hann) but never divided the window back out after the inverse, so every radial FFT filter result carried the window envelope — 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. fft_soft_border has always compensated its taper; fourier_filter now does the same, with the compensation gain capped at 20x so near-zero edge weights cannot amplify residual spectral leakage into spikes. The window="none" path is bit-identical to before (pinned). Note: replaying saved states that used a window now produces the corrected (envelope-free) output — same correctness-over-replay precedent as the gradient-slope and qPlus-setpoint fixes. - fft_soft_border's radial mask used index arithmetic centred at N/2.0, which is identical to the fftfreq convention for even sizes (pinned by a bit-equality test) but sat half a pixel off the fftshift DC bin for odd ones — a tiny low-pass could miss DC and shift the image mean. The mask now uses the same fftfreq-based cycles/pixel convention as fourier_filter. Verified clean in the same pass: fft_magnitude (axes, coherent gain, ROI mean handling), inverse_fft (raw unwindowed transform, conjugate-exact masks, imaginary-residual reporting), mains q-axes, the FFT viewer's q-axis construction and radial profile (per-axis physical d=, isotropic q-map), measurements/fft_points (explicit d= with honest cycles/pixel vs cycles/unit labels), and line_periodicity (rfftfreq with physical spacing, window-gain-normalised power, DC excluded). One suspicion was disproved by test: the soft-border half-axis normalisation is exactly equivalent to cycles/pixel (Δx/cx = 2f), not elliptical. Co-Authored-By: Claude Fable 5 --- probeflow/processing/filters.py | 22 +++++++-- tests/test_fft_physics.py | 80 +++++++++++++++++++++++++++++++++ 2 files changed, 98 insertions(+), 4 deletions(-) diff --git a/probeflow/processing/filters.py b/probeflow/processing/filters.py index 716c687..c772063 100644 --- a/probeflow/processing/filters.py +++ b/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": diff --git a/tests/test_fft_physics.py b/tests/test_fft_physics.py index 2a0a26d..a3cc811 100644 --- a/tests/test_fft_physics.py +++ b/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)