1#!/usr/bin/env python3
2"""
3Image Signal Processing: 2D DFT, Frequency-Domain and Spatial Filtering
4=========================================================================
5
6Images are 2-dimensional signals. Every technique from 1-D signal processing
7has a 2-D counterpart, and the connections are direct and elegant.
8
92-D Discrete Fourier Transform (DFT)
10--------------------------------------
11For an MรN image f(m, n):
12
13 F(u, v) = sum_{m=0}^{M-1} sum_{n=0}^{N-1}
14 f(m, n) * exp(-j*2*pi*(u*m/M + v*n/N))
15
16 - DC component F(0,0) = mean pixel value * M*N
17 - After fftshift: low frequencies are at the centre
18 - Convolution theorem: filtering in spatial domain โก multiplication in freq domain
19
202-D Filtering
21--------------
22Frequency domain:
23 F_filtered(u, v) = H(u, v) * F(u, v) (element-wise)
24 then ifft2 โ spatial result
25
26Spatial domain:
27 g(m, n) = f(m, n) ** h(m, n) (2-D convolution)
28
29Both are equivalent for linear shift-invariant (LSI) filters. The frequency-
30domain approach is often more intuitive for designing ideal brick-wall filters,
31while the spatial domain is faster for small kernels.
32
33Author: Educational example for Signal Processing
34License: MIT
35"""
36
37import numpy as np
38import matplotlib.pyplot as plt
39from scipy import ndimage
40
41
42# ============================================================================
43# SYNTHETIC TEST IMAGE
44# ============================================================================
45
46def make_test_image(size=256):
47 """
48 Create a synthetic test image combining:
49 - Checkerboard pattern (broad spectrum, high-frequency edges)
50 - Gaussian-blurred circles (low-frequency smooth features)
51 - Additive Gaussian noise
52
53 Args:
54 size (int): Image size (size ร size pixels)
55
56 Returns:
57 ndarray: Float64 image in [0, 1], shape (size, size)
58 """
59 img = np.zeros((size, size), dtype=float)
60
61 # Checkerboard: 16ร16 pixel squares
62 sq = size // 16
63 for i in range(16):
64 for j in range(16):
65 if (i + j) % 2 == 0:
66 img[i*sq:(i+1)*sq, j*sq:(j+1)*sq] = 0.6
67
68 # Smooth circles on top
69 y, x = np.ogrid[:size, :size]
70 cx, cy = size // 2, size // 2
71 for r, amp in [(60, 0.4), (30, -0.3), (15, 0.3)]:
72 mask = (x - cx) ** 2 + (y - cy) ** 2 < r ** 2
73 img[mask] += amp
74
75 # Gaussian blur on the circles to make them smooth
76 img = ndimage.gaussian_filter(img, sigma=3)
77
78 # Add noise
79 rng = np.random.default_rng(42)
80 img += 0.05 * rng.standard_normal(img.shape)
81
82 # Clip to [0, 1]
83 img = np.clip(img, 0, 1)
84 return img
85
86
87# ============================================================================
88# 2-D DFT UTILITIES
89# ============================================================================
90
91def compute_2d_spectrum(img):
92 """
93 Compute 2-D FFT and return the log-magnitude spectrum (shifted to centre).
94
95 Args:
96 img (ndarray): 2-D real image
97
98 Returns:
99 F (ndarray): Complex DFT coefficients (shifted)
100 mag (ndarray): Log-magnitude spectrum for display
101 """
102 F_raw = np.fft.fft2(img)
103 F = np.fft.fftshift(F_raw)
104 mag = np.log1p(np.abs(F)) # log(1 + |F|) avoids log(0)
105 return F, mag
106
107
108def freq_axes(img):
109 """
110 Return normalised frequency axes (cycles/pixel) for a shifted 2-D DFT.
111
112 Returns:
113 fu, fv (ndarray): 1-D frequency vectors, each in [-0.5, 0.5].
114 """
115 M, N = img.shape
116 fu = np.fft.fftshift(np.fft.fftfreq(M))
117 fv = np.fft.fftshift(np.fft.fftfreq(N))
118 return fu, fv
119
120
121# ============================================================================
122# FREQUENCY-DOMAIN FILTERS
123# ============================================================================
124
125def ideal_lowpass(img, cutoff):
126 """
127 Ideal (brick-wall) circular lowpass filter in the frequency domain.
128
129 H(u, v) = 1 if sqrt(u^2 + v^2) <= cutoff, else 0
130
131 Very sharp cutoff โ significant Gibbs ringing in spatial domain.
132
133 Args:
134 img (ndarray): Input image
135 cutoff (float) : Cutoff radius in normalised frequency [0, 0.5]
136
137 Returns:
138 filtered (ndarray): Filtered image (real part of IFFT)
139 H (ndarray): Filter mask (shifted), for display
140 """
141 F_raw = np.fft.fft2(img)
142 F = np.fft.fftshift(F_raw)
143
144 fu, fv = freq_axes(img)
145 UU, VV = np.meshgrid(fu, fv, indexing='ij')
146 dist = np.sqrt(UU ** 2 + VV ** 2)
147 H = (dist <= cutoff).astype(float)
148
149 F_filtered = np.fft.ifftshift(F * H)
150 filtered = np.real(np.fft.ifft2(F_filtered))
151 return filtered, H
152
153
154def gaussian_lowpass(img, sigma_freq):
155 """
156 Gaussian lowpass filter in the frequency domain.
157
158 H(u, v) = exp(-(u^2 + v^2) / (2*sigma^2))
159
160 Smooth roll-off โ no ringing. The spatial-domain kernel is also Gaussian.
161
162 Args:
163 img (ndarray): Input image
164 sigma_freq (float) : Standard deviation in normalised frequency
165
166 Returns:
167 filtered (ndarray): Filtered image
168 H (ndarray): Gaussian filter mask (shifted)
169 """
170 F_raw = np.fft.fft2(img)
171 F = np.fft.fftshift(F_raw)
172
173 fu, fv = freq_axes(img)
174 UU, VV = np.meshgrid(fu, fv, indexing='ij')
175 dist2 = UU ** 2 + VV ** 2
176 H = np.exp(-dist2 / (2 * sigma_freq ** 2))
177
178 F_filtered = np.fft.ifftshift(F * H)
179 filtered = np.real(np.fft.ifft2(F_filtered))
180 return filtered, H
181
182
183def ideal_highpass(img, cutoff):
184 """
185 Ideal highpass filter: H_hp = 1 - H_lp.
186
187 Removes low-frequency content; enhances edges and fine details.
188
189 Args:
190 img (ndarray): Input image
191 cutoff (float) : Cutoff radius in normalised frequency
192
193 Returns:
194 filtered (ndarray): Filtered image
195 H (ndarray): Highpass filter mask
196 """
197 F_raw = np.fft.fft2(img)
198 F = np.fft.fftshift(F_raw)
199
200 fu, fv = freq_axes(img)
201 UU, VV = np.meshgrid(fu, fv, indexing='ij')
202 dist = np.sqrt(UU ** 2 + VV ** 2)
203 H = (dist > cutoff).astype(float)
204
205 F_filtered = np.fft.ifftshift(F * H)
206 filtered = np.real(np.fft.ifft2(F_filtered))
207 return filtered, H
208
209
210# ============================================================================
211# SPATIAL-DOMAIN FILTERS
212# ============================================================================
213
214def sobel_edges(img):
215 """
216 Sobel edge detection using 3ร3 derivative kernels.
217
218 Gx detects horizontal changes, Gy detects vertical changes.
219 Edge magnitude: G = sqrt(Gx^2 + Gy^2)
220
221 Args:
222 img (ndarray): Grayscale image
223
224 Returns:
225 edges (ndarray): Edge magnitude image
226 """
227 Kx = np.array([[-1, 0, 1],
228 [-2, 0, 2],
229 [-1, 0, 1]], dtype=float)
230 Ky = Kx.T # transposing Kx gives the vertical kernel
231
232 Gx = ndimage.convolve(img, Kx)
233 Gy = ndimage.convolve(img, Ky)
234 edges = np.hypot(Gx, Gy)
235 return edges
236
237
238def gaussian_blur_spatial(img, sigma):
239 """
240 Gaussian blur via direct spatial convolution (SciPy implementation).
241
242 The 2-D Gaussian kernel:
243 h(m, n) = (1/(2*pi*sigma^2)) * exp(-(m^2+n^2)/(2*sigma^2))
244
245 Args:
246 img (ndarray): Input image
247 sigma (float) : Blur radius in pixels
248
249 Returns:
250 ndarray: Blurred image
251 """
252 return ndimage.gaussian_filter(img, sigma=sigma)
253
254
255# ============================================================================
256# PLOTTING HELPERS
257# ============================================================================
258
259def show_spectrum_pair(ax_img, ax_spec, image, title_img, title_spec, cmap='gray'):
260 """Display an image and its log-magnitude spectrum side-by-side on given axes."""
261 ax_img.imshow(image, cmap=cmap, vmin=0, vmax=1)
262 ax_img.set_title(title_img)
263 ax_img.axis('off')
264
265 _, mag = compute_2d_spectrum(image)
266 ax_spec.imshow(mag, cmap='hot')
267 ax_spec.set_title(title_spec)
268 ax_spec.axis('off')
269
270
271# ============================================================================
272# DEMO SECTIONS
273# ============================================================================
274
275def demo_spectrum(img):
276 """Visualise the test image and its 2-D DFT magnitude spectrum."""
277 print("=" * 60)
278 print("SECTION 1: 2-D DFT Magnitude Spectrum")
279 print("=" * 60)
280
281 F, mag = compute_2d_spectrum(img)
282 dc = np.abs(np.fft.fft2(img)[0, 0])
283 print(f" Image size : {img.shape}")
284 print(f" Mean pixel : {img.mean():.3f}")
285 print(f" DC magnitude F(0,0) = {dc:.1f} (โ mean ร {img.shape[0]} ร {img.shape[1]} = "
286 f"{img.mean()*img.shape[0]*img.shape[1]:.1f})")
287 print(f" Log spectrum range : [{mag.min():.2f}, {mag.max():.2f}]")
288
289 fig, axes = plt.subplots(1, 3, figsize=(14, 5))
290 fig.suptitle("Test Image and its 2-D DFT Magnitude Spectrum",
291 fontsize=12, fontweight='bold')
292
293 axes[0].imshow(img, cmap='gray', vmin=0, vmax=1)
294 axes[0].set_title("(a) Synthetic Test Image\n(checkerboard + circles + noise)")
295 axes[0].axis('off')
296
297 axes[1].imshow(mag, cmap='hot')
298 axes[1].set_title("(b) Log-Magnitude Spectrum |F(u,v)|\n"
299 "(DC at centre; bright = high energy)")
300 axes[1].axis('off')
301
302 # Annotate the spectrum image
303 cx, cy = np.array(img.shape) // 2
304 axes[1].plot(cy, cx, 'c+', markersize=12, markeredgewidth=2,
305 label='DC (u=v=0)')
306 axes[1].legend(loc='upper right', fontsize=8)
307
308 # 1-D cross-section through the spectrum
309 mid = mag.shape[0] // 2
310 axes[2].plot(mag[mid, :], 'C0', linewidth=1.2)
311 axes[2].set_title("(c) Horizontal Cross-section of Spectrum\n"
312 "at v=0; symmetric about centre")
313 axes[2].set_xlabel("Column index (u)")
314 axes[2].set_ylabel("Log magnitude")
315 axes[2].grid(True, alpha=0.3)
316
317 plt.tight_layout()
318 plt.savefig("15_spectrum.png", dpi=120)
319 print(" Saved: 15_spectrum.png")
320 plt.show()
321
322
323def demo_lowpass_filters(img):
324 """Compare ideal and Gaussian lowpass filters in the frequency domain."""
325 print("\n" + "=" * 60)
326 print("SECTION 2: 2-D Lowpass Filters (frequency domain)")
327 print("=" * 60)
328
329 cutoff = 0.08 # normalised frequency (0.5 = Nyquist)
330 sigma_g = 0.05
331
332 lp_ideal, H_ideal = ideal_lowpass(img, cutoff=cutoff)
333 lp_gauss, H_gauss = gaussian_lowpass(img, sigma_freq=sigma_g)
334
335 print(f" Ideal LP cutoff : {cutoff} (cycles/pixel)")
336 print(f" Gaussian LP sigma : {sigma_g} (cycles/pixel)")
337 print(f" Ideal LP output range : [{lp_ideal.min():.3f}, {lp_ideal.max():.3f}]")
338 print(f" Gaussian LP output range: [{lp_gauss.min():.3f}, {lp_gauss.max():.3f}]")
339
340 fig, axes = plt.subplots(2, 4, figsize=(16, 8))
341 fig.suptitle("2-D Lowpass Filters: Ideal vs Gaussian\n"
342 "Top row: filter masks | Bottom row: filtered images",
343 fontsize=12, fontweight='bold')
344
345 # Row 0: filter masks
346 axes[0, 0].imshow(img, cmap='gray', vmin=0, vmax=1)
347 axes[0, 0].set_title("Original Image")
348 axes[0, 0].axis('off')
349
350 _, orig_mag = compute_2d_spectrum(img)
351 axes[0, 1].imshow(orig_mag, cmap='hot')
352 axes[0, 1].set_title("Original Spectrum")
353 axes[0, 1].axis('off')
354
355 axes[0, 2].imshow(H_ideal, cmap='gray')
356 axes[0, 2].set_title(f"Ideal LP Mask\n(cutoff={cutoff})")
357 axes[0, 2].axis('off')
358
359 axes[0, 3].imshow(H_gauss, cmap='gray')
360 axes[0, 3].set_title(f"Gaussian LP Mask\n(ฯ={sigma_g})")
361 axes[0, 3].axis('off')
362
363 # Row 1: filtered results
364 lp_ideal_clipped = np.clip(lp_ideal, 0, 1)
365 lp_gauss_clipped = np.clip(lp_gauss, 0, 1)
366
367 axes[1, 0].imshow(lp_ideal_clipped, cmap='gray', vmin=0, vmax=1)
368 axes[1, 0].set_title("Ideal LP Result\n(note: Gibbs ringing)")
369 axes[1, 0].axis('off')
370
371 _, mag_ideal = compute_2d_spectrum(lp_ideal_clipped)
372 axes[1, 1].imshow(mag_ideal, cmap='hot')
373 axes[1, 1].set_title("Ideal LP Spectrum\n(high freqs removed)")
374 axes[1, 1].axis('off')
375
376 axes[1, 2].imshow(lp_gauss_clipped, cmap='gray', vmin=0, vmax=1)
377 axes[1, 2].set_title("Gaussian LP Result\n(smooth, no ringing)")
378 axes[1, 2].axis('off')
379
380 _, mag_gauss = compute_2d_spectrum(lp_gauss_clipped)
381 axes[1, 3].imshow(mag_gauss, cmap='hot')
382 axes[1, 3].set_title("Gaussian LP Spectrum")
383 axes[1, 3].axis('off')
384
385 plt.tight_layout()
386 plt.savefig("15_lowpass_filters.png", dpi=120)
387 print(" Saved: 15_lowpass_filters.png")
388 plt.show()
389
390 return lp_ideal, lp_gauss
391
392
393def demo_highpass_and_edges(img):
394 """
395 Apply highpass filtering in the frequency domain and compare with
396 spatial-domain Sobel edge detection.
397 """
398 print("\n" + "=" * 60)
399 print("SECTION 3: Highpass Filter vs Spatial Sobel Edge Detection")
400 print("=" * 60)
401
402 hp_img, H_hp = ideal_highpass(img, cutoff=0.05)
403 sobel_img = sobel_edges(img)
404
405 # Normalise for display
406 hp_display = np.abs(hp_img)
407 hp_display = hp_display / hp_display.max()
408 sobel_display = sobel_img / sobel_img.max()
409
410 print(f" HP filter output max : {hp_img.max():.4f}")
411 print(f" Sobel edges output max: {sobel_img.max():.4f}")
412
413 fig, axes = plt.subplots(2, 3, figsize=(14, 9))
414 fig.suptitle("Highpass Filtering (Frequency Domain) vs Sobel Edge Detection (Spatial)",
415 fontsize=12, fontweight='bold')
416
417 axes[0, 0].imshow(img, cmap='gray', vmin=0, vmax=1)
418 axes[0, 0].set_title("(a) Original Image")
419 axes[0, 0].axis('off')
420
421 axes[0, 1].imshow(H_hp, cmap='gray')
422 axes[0, 1].set_title("(b) Ideal HP Mask\n(zeros at centre = DC blocked)")
423 axes[0, 1].axis('off')
424
425 axes[0, 2].imshow(hp_display, cmap='gray', vmin=0, vmax=1)
426 axes[0, 2].set_title("(c) Ideal HP Output\n(edges + fine detail)")
427 axes[0, 2].axis('off')
428
429 # Spectrum of HP output
430 _, mag_hp = compute_2d_spectrum(hp_display)
431 axes[1, 0].imshow(mag_hp, cmap='hot')
432 axes[1, 0].set_title("(d) HP Output Spectrum\n(DC gap visible at centre)")
433 axes[1, 0].axis('off')
434
435 axes[1, 1].imshow(sobel_display, cmap='gray', vmin=0, vmax=1)
436 axes[1, 1].set_title("(e) Sobel Edge Magnitude\n(spatial 3ร3 kernel)")
437 axes[1, 1].axis('off')
438
439 # Overlay comparison
440 axes[1, 2].imshow(img, cmap='gray', vmin=0, vmax=1, alpha=0.6)
441 axes[1, 2].imshow(sobel_display, cmap='Reds', alpha=0.5, vmin=0.1, vmax=1)
442 axes[1, 2].set_title("(f) Sobel Edges Overlaid on Original")
443 axes[1, 2].axis('off')
444
445 plt.tight_layout()
446 plt.savefig("15_highpass_edges.png", dpi=120)
447 print(" Saved: 15_highpass_edges.png")
448 plt.show()
449
450
451def demo_freq_vs_spatial_blur(img):
452 """
453 Demonstrate equivalence of frequency-domain Gaussian LP filter and
454 spatial-domain Gaussian blur, and show the effect of increasing blur.
455 """
456 print("\n" + "=" * 60)
457 print("SECTION 4: Frequency Domain vs Spatial Domain Blur Comparison")
458 print("=" * 60)
459
460 # Gaussian LP in frequency domain (sigma_freq โ sigma_spatial via uncertainty)
461 # Relationship: sigma_spatial * sigma_freq โ 1/(2*pi)
462 sigma_freq = 0.04
463 sigma_spatial = 1.0 / (2 * np.pi * sigma_freq) # theoretical equivalence
464 print(f" Gaussian LP sigma_freq={sigma_freq:.3f} โ "
465 f"spatial sigmaโ{sigma_spatial:.1f} pixels")
466
467 blurred_freq, _ = gaussian_lowpass(img, sigma_freq=sigma_freq)
468 blurred_spatial = gaussian_blur_spatial(img, sigma=sigma_spatial)
469
470 diff = np.abs(blurred_freq - blurred_spatial)
471 print(f" Max pixel difference (freq vs spatial): {diff.max():.6f}")
472 print(" (small difference confirms near-equivalence for this image size)")
473
474 # Show progression of spatial blur
475 sigmas = [1, 4, 10]
476
477 fig, axes = plt.subplots(2, 3, figsize=(14, 9))
478 fig.suptitle("Frequency Domain vs Spatial Domain Gaussian Blur",
479 fontsize=12, fontweight='bold')
480
481 # Top row: spatial blur at increasing sigma
482 for ax, sigma in zip(axes[0], sigmas):
483 blurred = gaussian_blur_spatial(img, sigma=sigma)
484 ax.imshow(np.clip(blurred, 0, 1), cmap='gray', vmin=0, vmax=1)
485 ax.set_title(f"Spatial Blur ฯ={sigma} px")
486 ax.axis('off')
487
488 # Bottom row: their spectra
489 for ax, sigma in zip(axes[1], sigmas):
490 blurred = gaussian_blur_spatial(img, sigma=sigma)
491 _, mag = compute_2d_spectrum(blurred)
492 ax.imshow(mag, cmap='hot')
493 ax.set_title(f"Spectrum (ฯ={sigma} px)\n"
494 f"Larger ฯ โ narrower spectrum")
495 ax.axis('off')
496
497 plt.tight_layout()
498 plt.savefig("15_blur_comparison.png", dpi=120)
499 print(" Saved: 15_blur_comparison.png")
500 plt.show()
501
502
503# ============================================================================
504# MAIN
505# ============================================================================
506
507if __name__ == "__main__":
508 print("Image Signal Processing: 2D DFT and Filtering")
509 print("=" * 60)
510 print("Key relationships:")
511 print(" Convolution theorem : f*h <-> FยทH (spatial โ frequency)")
512 print(" Large ฯ spatial โ narrow Gaussian in frequency (smooth image)")
513 print(" Ideal LP (brick wall) โ sinc ripple in spatial domain (Gibbs)")
514 print(" Gaussian LP โ Gaussian in spatial domain (no ringing)")
515 print()
516
517 size = 256
518 print(f"Creating synthetic {size}ร{size} test image...")
519 img = make_test_image(size=size)
520 print(f" Image stats: min={img.min():.3f}, max={img.max():.3f}, "
521 f"mean={img.mean():.3f}")
522
523 demo_spectrum(img)
524 demo_lowpass_filters(img)
525 demo_highpass_and_edges(img)
526 demo_freq_vs_spatial_blur(img)
527
528 print("\nDone. Four PNG files saved.")
529 print("\nKey takeaways:")
530 print(" - 2-D DFT: low freqs at centre (after fftshift), "
531 "high freqs at edges")
532 print(" - Ideal LP: sharp mask โ Gibbs ringing in spatial domain")
533 print(" - Gaussian LP: smooth mask โ smooth blur, no ringing")
534 print(" - Highpass = 1 - Lowpass: emphasises edges and fine detail")
535 print(" - Frequency and spatial domain filtering are theoretically "
536 "equivalent (convolution theorem)")