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)