Lesson 14: Nonlinear Filters¶
Every filter so far (box, Gaussian, Sobel, Laplacian, ...) has been a convolution: a fixed, linear combination of neighboring pixel values. But linear filters are limited in what they can do — they can't reason about which neighboring pixels are "outliers" or which are "on the other side of an edge." This lesson looks at filters that break that linearity: median filtering, grayscale erosion/dilation as min/max filters, and the edge-preserving bilateral filter.
import numpy as np
import cv2
import matplotlib.pyplot as plt
Median filtering: robust to outliers¶
A median filter replaces each pixel with the median (not the mean) of its neighborhood. Whereas a single wildly wrong pixel value (outlier) can drastically affect the mean, it barely moves a median, if at all — so median filtering is especially effective against salt-and-pepper noise (isolated pixels randomly slammed to black or white), which linear smoothing handles poorly.
clean = np.zeros((120, 120), dtype=np.uint8)
cv2.rectangle(clean, (25, 25), (95, 95), 200, -1)
rng = np.random.default_rng(2)
salt_pepper = clean.copy()
hit = rng.random(clean.shape) < 0.06
salt_pepper[hit & (rng.random(clean.shape) < 0.5)] = 0
salt_pepper[hit & (rng.random(clean.shape) >= 0.5)] = 255
median_result = cv2.medianBlur(salt_pepper, 5)
gaussian_result = cv2.GaussianBlur(salt_pepper, (5, 5), sigmaX=1.5)
fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))
for ax, im, title in zip(axes, [salt_pepper, gaussian_result, median_result],
['Salt-and-pepper noise', 'Gaussian blur', 'Median filter']):
ax.imshow(im, cmap='gray')
ax.set_title(title, fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
print('mean abs error vs. clean image:')
print(f' Gaussian blur : {np.abs(gaussian_result.astype(int) - clean.astype(int)).mean():.2f}')
print(f' Median filter : {np.abs(median_result.astype(int) - clean.astype(int)).mean():.2f}')
Gaussian blur smears each corrupted pixel's extreme value into its neighbors, leaving faint speckle everywhere. The median filter simply outvotes an isolated bad pixel with its many good neighbors, removing nearly all of the noise while keeping edges sharp.
Min and max filters: grayscale erosion and dilation¶
Lesson 3 introduced erosion and dilation on binary images. On a grayscale image, they generalize naturally: erosion replaces each pixel with the minimum value in its neighborhood, and dilation with the maximum. Both are nonlinear (not a weighted sum), and both are edge-preserving in a different way than blurring — they shift edges rather than blurring them.
def min_filter(image, ksize):
r = ksize // 2
padded = cv2.copyMakeBorder(image, r, r, r, r, cv2.BORDER_REPLICATE)
out = np.zeros_like(image)
for i in range(image.shape[0]):
for j in range(image.shape[1]):
out[i, j] = padded[i:i + ksize, j:j + ksize].min()
return out
small = clean[::4, ::4] # smaller image so the manual pixel loop stays fast
mine = min_filter(small, 3)
reference = cv2.erode(small, np.ones((3, 3), np.uint8))
print('manual min filter == cv2.erode?', np.array_equal(mine, reference))
The problem linear blur can't solve: smoothing without crossing edges¶
Gaussian blur treats every neighbor equally regardless of how different its value is — that's exactly what makes it blur across edges. The bilateral filter fixes this by weighting each neighbor by two factors: how close it is spatially (like a Gaussian blur) and how similar its intensity is to the center pixel:
$$I'(p) = \frac{1}{W_p}\sum_{q \in N(p)} G_{\sigma_s}(\|p - q\|) \cdot G_{\sigma_r}(|I(p) - I(q)|) \cdot I(q)$$
where $W_p$ normalizes the weights to sum to 1. Neighbors that are spatially close and have similar intensity get high weight (smoothed together); neighbors that are spatially close but have very different intensity (i.e. on the other side of an edge) get down-weighted almost to zero, so the edge survives.
def bilateral_filter(image, radius, sigma_spatial, sigma_range):
image_f = image.astype(np.float64)
ys, xs = np.mgrid[-radius:radius + 1, -radius:radius + 1]
spatial_weight = np.exp(-(xs**2 + ys**2) / (2 * sigma_spatial**2))
padded = cv2.copyMakeBorder(image_f, radius, radius, radius, radius, cv2.BORDER_REFLECT101)
out = np.zeros_like(image_f)
h, w = image.shape
for i in range(h):
for j in range(w):
patch = padded[i:i + 2 * radius + 1, j:j + 2 * radius + 1]
range_weight = np.exp(-(patch - image_f[i, j])**2 / (2 * sigma_range**2))
weight = spatial_weight * range_weight
out[i, j] = (weight * patch).sum() / weight.sum()
return out
Sanity check against OpenCV¶
cv2.bilateralFilter uses the same idea with some internal approximations for speed, so we expect a close but not pixel-exact match.
step = np.zeros((40, 40), dtype=np.uint8)
step[20:, :] = 200
noisy_step = np.clip(step.astype(np.float64) + rng.normal(0, 10, step.shape), 0, 255).astype(np.uint8)
mine = bilateral_filter(noisy_step, radius=3, sigma_spatial=3, sigma_range=25)
cv_result = cv2.bilateralFilter(noisy_step, d=7, sigmaColor=25, sigmaSpace=3).astype(np.float64)
print(f'max abs difference = {np.abs(mine - cv_result).max():.2f} gray levels')
print(f'mean abs difference = {np.abs(mine - cv_result).mean():.2f} gray levels')
Edge preservation, side by side¶
The real comparison that matters: on a noisy step edge, does the filter smooth the noise while keeping the edge sharp, or does it blur the edge along with the noise?
big_step = np.zeros((150, 150), dtype=np.uint8)
big_step[75:, :] = 200
noisy_big = np.clip(big_step.astype(np.float64) + rng.normal(0, 15, big_step.shape), 0, 255).astype(np.uint8)
gaussian_smoothed = cv2.GaussianBlur(noisy_big, (0, 0), sigmaX=4)
bilateral_smoothed = cv2.bilateralFilter(noisy_big, d=9, sigmaColor=40, sigmaSpace=15)
fig, axes = plt.subplots(2, 3, figsize=(10, 6))
for ax, im, title in zip(axes[0], [noisy_big, gaussian_smoothed, bilateral_smoothed],
['Noisy step edge', 'Gaussian blur', 'Bilateral filter']):
ax.imshow(im, cmap='gray', vmin=0, vmax=255)
ax.set_title(title, fontsize=9)
ax.axis('off')
col = 75
for ax, im, title in zip(axes[1], [noisy_big, gaussian_smoothed, bilateral_smoothed],
['profile: noisy', 'profile: Gaussian', 'profile: bilateral']):
ax.plot(im[:, col])
ax.set_ylim(0, 255)
ax.set_title(title, fontsize=9)
plt.tight_layout()
plt.show()
The Gaussian-blurred profile shows a gradual ramp across many rows — the edge has been visibly softened. The bilateral-filtered profile stays flat on each side (noise removed) but still jumps sharply near row 75 — the edge survived because pixels across the boundary were too different in intensity to be averaged together, no matter how spatially close they were.
Exercise¶
- Apply a median filter (instead of Gaussian or bilateral) to
noisy_bigand add its profile to the comparison. How does it compare to the bilateral filter at preserving the edge? - In
bilateral_filter, setsigma_rangevery large (e.g. 1000). What should this reduce to, and does your result confirm it? - Bilateral filtering is much slower than Gaussian blur because the weights must be recomputed at every pixel (they depend on local intensity, not just position). Time the from-scratch
bilateral_filteragainstcv2.bilateralFilteron a 100x100 image with%timeit, and describe in words what OpenCV's implementation likely does differently to be faster.