Practical 8: Spatial Filtering — Smoothing & Sharpening
Run this practical in Google Colab
Aim
To implement and analyze the four canonical spatial-domain filters of digital image processing: the box (averaging) filter and median filter for noise smoothing, and the Laplacian and Sobel operators for sharpening and edge detection.
2. Description / Theory
Theory:
- Box (mean): \( \hat f(x,y) = \frac{1}{mn} \sum_{(s,t)\in S_{xy}} f(s,t) \) — linear smoothing.
- Median: \( \hat f(x,y) = \mathrm{median}\{ f(s,t) : (s,t)\in S_{xy} \} \) — non-linear, preserves edges, ideal for impulse noise.
- Laplacian: \( \nabla^2 f \approx f_{N}+f_{S}+f_{E}+f_{W} - 4f(x,y) \) — second derivative for sharpening: \( g = f - \nabla^2 f \).
- Sobel: first-derivative gradient via \(G_x, G_y\); magnitude \(|G| \approx |G_x|+|G_y|\); direction \(\theta = \arctan(G_y/G_x)\).
# ----------------------------------------------------------------------
# Install dependencies (uncomment for Google Colab)
# !pip install opencv-python-headless matplotlib numpy scipy
import cv2
import numpy as np
import matplotlib.pyplot as plt
import os
# === AUTO-DOWNLOAD DATASET (works in Google Colab and locally) ===
import os, urllib.request, zipfile
CHAPTER = "CH03" # Chapter 3: Spatial Filtering
DATASET_PATH = f"datasets/{CHAPTER}/"
DOWNLOAD_BASE = "https://www.imageprocessingplace.com/downloads_V3/dip3e_downloads/dip3e_book_images"
if not os.path.exists(DATASET_PATH) or not any(f.endswith('.tif') for f in os.listdir(DATASET_PATH)):
zip_name = f"DIP3E_{CHAPTER}_Original_Images.zip"
url = f"{DOWNLOAD_BASE}/{zip_name}"
print(f"Downloading {CHAPTER} dataset from imageprocessingplace.com...")
urllib.request.urlretrieve(url, "chapter.zip")
os.makedirs(DATASET_PATH, exist_ok=True)
with zipfile.ZipFile("chapter.zip", "r") as z:
for f in z.namelist():
if f.lower().endswith(".tif"):
fname = os.path.basename(f)
if fname:
with z.open(f) as src, open(os.path.join(DATASET_PATH, fname), "wb") as dst:
dst.write(src.read())
os.remove("chapter.zip")
print(f"Downloaded {len([f for f in os.listdir(DATASET_PATH) if f.endswith('.tif')])} images")
else:
print(f"Dataset ready: {len([f for f in os.listdir(DATASET_PATH) if f.endswith('.tif')])} images")
images = sorted([f for f in os.listdir(DATASET_PATH) if f.endswith('.tif')])
print(f"\nAvailable images ({len(images)}):")
for i, name in enumerate(images, 1):
print(f" {i:2d}. {name}")
def box_filter(img, k=3):
"""Manual box filter via 2D correlation with a uniform kernel."""
p = k // 2
kernel = np.ones((k, k), dtype=np.float64) / (k * k)
padded = np.pad(img.astype(np.float64), p, mode='edge')
out = np.zeros_like(img, dtype=np.float64)
h, w = img.shape
for i in range(h):
for j in range(w):
out[i, j] = np.sum(padded[i:i+k, j:j+k] * kernel)
return np.clip(out, 0, 255).astype(np.uint8)
test_pattern_file = "Fig0333(a)(test_pattern_blurring_orig).tif"
tp = cv2.imread(os.path.join(DATASET_PATH, test_pattern_file), cv2.IMREAD_GRAYSCALE)
print(f"Test pattern: {test_pattern_file}, shape={tp.shape}")
kernel_sizes = [3, 5, 9, 15, 35]
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.ravel()
axes[0].imshow(tp, cmap='gray', vmin=0, vmax=255)
axes[0].set_title('Original\nFig0333(a)'); axes[0].axis('off')
for idx, k in enumerate(kernel_sizes, start=1):
out = box_filter(tp, k=k)
axes[idx].imshow(out, cmap='gray', vmin=0, vmax=255)
axes[idx].set_title(f'Box {k}x{k}'); axes[idx].axis('off')
plt.suptitle('Box Filter: Effect of Kernel Size on Test Pattern', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
def median_filter(img, k=3):
"""Manual median filter via sliding window."""
p = k // 2
padded = np.pad(img, p, mode='edge')
out = np.zeros_like(img)
h, w = img.shape
for i in range(h):
for j in range(w):
out[i, j] = np.median(padded[i:i+k, j:j+k])
return out.astype(np.uint8)
salt_pepper_file = "Fig0335(a)(ckt_board_saltpep_prob_pt05).tif"
if salt_pepper_file not in images:
# Some distributions name it differently
candidates = [f for f in images if 'ckt' in f.lower() and ('saltpep' in f.lower() or 'saltpr' in f.lower())]
salt_pepper_file = candidates[0] if candidates else images[0]
print(f"Substituted: {salt_pepper_file}")
noisy = cv2.imread(os.path.join(DATASET_PATH, salt_pepper_file), cv2.IMREAD_GRAYSCALE)
print(f"Noisy image: {salt_pepper_file}, shape={noisy.shape}")
fig, axes = plt.subplots(1, 5, figsize=(20, 5))
axes[0].imshow(noisy, cmap='gray', vmin=0, vmax=255)
axes[0].set_title('Original\n(Salt & Pepper, p=0.05)'); axes[0].axis('off')
for idx, k in enumerate([3, 5, 7, 9], start=1):
out = median_filter(noisy, k=k)
axes[idx].imshow(out, cmap='gray', vmin=0, vmax=255)
axes[idx].set_title(f'Median {k}x{k}'); axes[idx].axis('off')
plt.suptitle('Median Filter: Removing Impulse Noise', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
k = 3
box_out = box_filter(noisy, k=k)
med_out = median_filter(noisy, k=k)
removed_box = np.abs(noisy.astype(int) - box_out.astype(int)).astype(np.uint8)
removed_med = np.abs(noisy.astype(int) - med_out.astype(int)).astype(np.uint8)
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes[0, 0].imshow(noisy, cmap='gray', vmin=0, vmax=255); axes[0, 0].set_title('Noisy input'); axes[0, 0].axis('off')
axes[0, 1].imshow(box_out, cmap='gray', vmin=0, vmax=255); axes[0, 1].set_title(f'Box {k}x{k}'); axes[0, 1].axis('off')
axes[0, 2].imshow(med_out, cmap='gray', vmin=0, vmax=255); axes[0, 2].set_title(f'Median {k}x{k}'); axes[0, 2].axis('off')
axes[1, 0].axis('off')
axes[1, 1].imshow(removed_box, cmap='hot'); axes[1, 1].set_title('Removed by Box (noise + edges)'); axes[1, 1].axis('off')
axes[1, 2].imshow(removed_med, cmap='hot'); axes[1, 2].set_title('Removed by Median (mostly noise)'); axes[1, 2].axis('off')
plt.suptitle('Box vs Median: Edge Preservation Under Impulse Noise', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
print(f"Mean energy removed by box = {removed_box.mean():.2f}")
print(f"Mean energy removed by median = {removed_med.mean():.2f}")
lap4 = np.array([[ 0, -1, 0],
[-1, 4, -1],
[ 0, -1, 0]], dtype=np.float64)
lap8 = np.array([[-1, -1, -1],
[-1, 8, -1],
[-1, -1, -1]], dtype=np.float64)
moon_file = "Fig0338(a)(blurry_moon).tif"
if moon_file not in images:
candidates = [f for f in images if 'moon' in f.lower()]
moon_file = candidates[0] if candidates else images[0]
print(f"Substituted: {moon_file}")
moon = cv2.imread(os.path.join(DATASET_PATH, moon_file), cv2.IMREAD_GRAYSCALE)
print(f"Moon image: {moon_file}, shape={moon.shape}")
lap4_resp = cv2.filter2D(moon.astype(np.float64), ddepth=-1, kernel=lap4)
lap8_resp = cv2.filter2D(moon.astype(np.float64), ddepth=-1, kernel=lap8)
sharp4 = np.clip(moon.astype(np.float64) - lap4_resp, 0, 255).astype(np.uint8)
sharp8 = np.clip(moon.astype(np.float64) - lap8_resp, 0, 255).astype(np.uint8)
def to_disp(a):
a = a - a.min()
return (255.0 * a / max(a.max(), 1e-9)).astype(np.uint8)
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes[0, 0].imshow(moon, cmap='gray', vmin=0, vmax=255); axes[0, 0].set_title('Original (blurry moon)'); axes[0, 0].axis('off')
axes[0, 1].imshow(to_disp(lap4_resp), cmap='gray'); axes[0, 1].set_title('Laplacian response (4-neighbour)'); axes[0, 1].axis('off')
axes[0, 2].imshow(sharp4, cmap='gray', vmin=0, vmax=255); axes[0, 2].set_title('Sharpened: f - lap4(f)'); axes[0, 2].axis('off')
axes[1, 0].imshow(moon, cmap='gray', vmin=0, vmax=255); axes[1, 0].set_title('Original (blurry moon)'); axes[1, 0].axis('off')
axes[1, 1].imshow(to_disp(lap8_resp), cmap='gray'); axes[1, 1].set_title('Laplacian response (8-neighbour)'); axes[1, 1].axis('off')
axes[1, 2].imshow(sharp8, cmap='gray', vmin=0, vmax=255); axes[1, 2].set_title('Sharpened: f - lap8(f)'); axes[1, 2].axis('off')
plt.suptitle('Laplacian Sharpening: 4-neighbour vs 8-neighbour kernel', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
sx = np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]], dtype=np.float64)
sy = np.array([[-1, -2, -1],
[ 0, 0, 0],
[ 1, 2, 1]], dtype=np.float64)
lens_file = "Fig0342(a)(contact_lens_original).tif"
if lens_file not in images:
candidates = [f for f in images if 'lens' in f.lower() or 'contact' in f.lower()]
lens_file = candidates[0] if candidates else images[0]
print(f"Substituted: {lens_file}")
lens = cv2.imread(os.path.join(DATASET_PATH, lens_file), cv2.IMREAD_GRAYSCALE)
print(f"Edge-detection image: {lens_file}, shape={lens.shape}")
Gx = cv2.filter2D(lens.astype(np.float64), ddepth=-1, kernel=sx)
Gy = cv2.filter2D(lens.astype(np.float64), ddepth=-1, kernel=sy)
Gmag = np.abs(Gx) + np.abs(Gy)
Gmag_disp = np.clip(255.0 * Gmag / max(Gmag.max(), 1e-9), 0, 255).astype(np.uint8)
Gtheta = np.arctan2(Gy, Gx) # radians, range [-pi, pi]
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes[0, 0].imshow(lens, cmap='gray', vmin=0, vmax=255); axes[0, 0].set_title('Original'); axes[0, 0].axis('off')
axes[0, 1].imshow(np.abs(Gx), cmap='gray'); axes[0, 1].set_title('|Gx| (vertical edges)'); axes[0, 1].axis('off')
axes[0, 2].imshow(np.abs(Gy), cmap='gray'); axes[0, 2].set_title('|Gy| (horizontal edges)'); axes[0, 2].axis('off')
axes[1, 0].imshow(Gmag_disp, cmap='gray'); axes[1, 0].set_title('|G| = |Gx|+|Gy|'); axes[1, 0].axis('off')
th_disp = ((Gtheta + np.pi) / (2 * np.pi) * 255).astype(np.uint8)
axes[1, 1].imshow(th_disp, cmap='hsv'); axes[1, 1].set_title('Gradient direction \u03B8'); axes[1, 1].axis('off')
thresh = (Gmag_disp > 64).astype(np.uint8) * 255
axes[1, 2].imshow(thresh, cmap='gray', vmin=0, vmax=255); axes[1, 2].set_title('Thresholded edges'); axes[1, 2].axis('off')
plt.suptitle('Sobel Edge Detection: Gx, Gy, magnitude, direction, threshold', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()Part 1: Box Filter (Smoothing)
Apply the box filter at multiple kernel sizes (3, 5, 9, 15, 35) on the standard test pattern. Larger kernels produce stronger blur because they low-pass at lower cut-off frequencies.
Part 2: Median Filter on Salt-and-Pepper Noise
Apply the median filter at multiple kernel sizes to the circuit-board image corrupted with 5% salt-and-pepper noise.
Part 3: Box vs Median — Edge Preservation
Apply both filters at the same kernel size to the same noisy input. The box filter blurs noise and edges together; the median filter removes noise while preserving edges.
Part 4: Laplacian Sharpening
Compute \(\nabla^2 f\) with the 4-neighbour and 8-neighbour kernels and form the sharpened image \(g = f - \nabla^2 f\). Demonstrated on the blurry-moon image.
Part 5: Sobel Edge Detection
Compute \(G_x\), \(G_y\), gradient magnitude \(|G|\) and direction \(\theta\) on the contact-lens image, then threshold the magnitude to obtain a clean edge map.
Analysis Questions
- As the box-filter kernel grows from 3×3 to 35×35, what happens to the image, and why does this happen at \(O(1/k)\) in the spatial-frequency sense?
- Why does the median filter remove salt-and-pepper noise without blurring edges, while the box filter blurs both?
- The Laplacian is a second-derivative operator. Why does subtracting its response sharpen the image rather than darken it? How does the 8-neighbour kernel differ from the 4-neighbour version?
- Sobel uses two separate kernels and combines their magnitudes. Why is it preferred over the Laplacian for edge detection, even though both involve derivatives?