Edge-Preserving Smoothing: Bilateral & Guided Filters

"I average with pixels from my side of the edge. The pixels across the boundary seem lovely, truly, but we are 80 gray levels apart and it would never work."

A Diplomatically Selective Bilateral Filter

The smoothers of Section 3.2 obeyed the shift-invariance contract from Section 3.1: same weights everywhere, no exceptions. This section breaks the contract deliberately, and it is worth being precise about why mere parameter tuning could never be enough.

1. Why Every Linear Filter Must Blur Edges Intermediate

Suppose a filter is linear and shift-invariant, applies positive weights summing to one, and is asked to process a perfect step edge: a row of dark pixels at value 50 meeting bright pixels at value 150. Any window that straddles the boundary contains pixels from both populations, and a positive weighted average of values from both sides must land strictly between them. The output therefore contains intermediate values, 80, 110, 130, where the input had none: the step has become a ramp. No choice of positive weights avoids this; wider noise-fighting windows only widen the ramp. This is not an implementation defect but a small theorem: linear smoothing and edge preservation are mathematically incompatible.

The escape requires weights that depend on image content, so that a window straddling an edge averages only the pixels on the center's own side. Content-dependent weights make the filter nonlinear (we saw a first preview with the median in Section 3.2), and they cost us the tidy LSI theory of Section 3.1, including the single-kernel representation and the frequency-domain analysis of Chapter 4. What we buy with that sacrifice is the subject of this section.

2. The Bilateral Filter: Two Gaussians, One Veto Intermediate

The bilateral filter (Tomasi and Manduchi, 1998) multiplies two notions of closeness. A spatial Gaussian $G_{\sigma_s}$ weights neighbors by distance, exactly as in Section 3.2. A range Gaussian $G_{\sigma_r}$ weights them by how similar their intensity is to the center pixel's:

$$ B(x) = \frac{1}{W(x)} \sum_{y \in \Omega} G_{\sigma_s}\!\big(\|x - y\|\big)\; G_{\sigma_r}\!\big(|I(x) - I(y)|\big)\; I(y) $$

where $W(x)$ sums the combined weights so they normalize to one at each pixel. The range term is the veto: a neighbor whose intensity differs from the center by much more than $\sigma_r$ contributes essentially nothing, however close it sits spatially. Figure 3.5.1 shows the consequence at an edge: the effective kernel is a Gaussian with one side amputated, so each pixel averages only within its own region and the step survives smoothing untouched.

Parameter intuition: $\sigma_s$ plays the familiar Gaussian role (how far to look), while $\sigma_r$ defines what counts as "the same region" in gray levels. Set $\sigma_r$ well below the contrast of edges you must preserve and above the noise amplitude you must remove; the filter then lives in the comfortable gap between the two. As $\sigma_r \to \infty$ the veto never fires and the bilateral degenerates into the plain Gaussian; as $\sigma_r \to 0$ it refuses to average at all.

import cv2

img = cv2.imread("portrait.jpg")

# d=9: window diameter. sigmaColor: range veto (gray levels).
# sigmaSpace: spatial reach (pixels).
smooth = cv2.bilateralFilter(img, d=9, sigmaColor=40, sigmaSpace=5)

# Compare against a plain Gaussian of similar spatial reach:
gauss = cv2.GaussianBlur(img, (9, 9), 2)

edges_b = cv2.Laplacian(cv2.cvtColor(smooth, cv2.COLOR_BGR2GRAY), cv2.CV_32F).var()
edges_g = cv2.Laplacian(cv2.cvtColor(gauss,  cv2.COLOR_BGR2GRAY), cv2.CV_32F).var()
print(f"edge energy: bilateral {edges_b:.0f} vs gaussian {edges_g:.0f}")
# Representative output:
# edge energy: bilateral 142 vs gaussian 61
# The bilateral result keeps over twice the Laplacian energy:
# noise is gone from flat regions but the edges stayed sharp.

Two costs come with the magic. First, speed: weights must be recomputed at every pixel, so the naive bilateral is $O(k^2)$ per pixel with no separability rescue (fast approximations exist, and OpenCV's implementation is heavily optimized, but it remains far costlier than GaussianBlur). Second, a signature artifact: on smooth gradients punctuated by noise, the range veto can quantize the gradient into terraces, the staircasing or "cartoon" look, which heavy-handed beautify filters exhibit in the wild.

3. The Guided Filter: Edge Awareness at Box-Filter Prices Advanced

The guided filter (He, Sun, and Tang, 2013) achieves a similar goal through entirely different machinery, and its cost analysis explains its industrial dominance. The idea: assume that within each small window $\omega_k$, the output $q$ is a linear function of a guidance image $G$ (often the input itself): $q_i = a_k G_i + b_k$ for $i \in \omega_k$. Solving the least-squares fit per window gives closed forms:

$$ a_k = \frac{\operatorname{cov}_{\omega_k}(G, I)}{\operatorname{var}_{\omega_k}(G) + \epsilon}, \qquad b_k = \bar{I}_k - a_k \bar{G}_k $$

Where the guide is flat (variance far below $\epsilon$), $a_k \to 0$ and the filter outputs the local mean: full smoothing. Where the guide has a strong edge (variance far above $\epsilon$), $a_k \to 1$ and the output tracks the guide: the edge passes through intact. The regularizer $\epsilon$ thus plays exactly the role $\sigma_r^2$ plays for the bilateral: the variance threshold separating "texture to smooth" from "structure to keep". Every quantity in these formulas is a local mean, computable with box filters, and Section 3.6 shows box filters run in $O(1)$ per pixel regardless of radius. The punchline: edge-preserving smoothing whose cost does not depend on the window size at all.

import cv2
import numpy as np

def box(x, r):
    """Mean filter with window (2r+1)x(2r+1); O(1) per pixel."""
    return cv2.boxFilter(x, -1, (2 * r + 1, 2 * r + 1))

def guided_filter(G, I, r=8, eps=0.01):
    """He et al. guided filter. G: guide, I: input, both float32 in [0,1]."""
    mean_G, mean_I = box(G, r), box(I, r)
    corr_GI = box(G * I, r)
    var_G   = box(G * G, r) - mean_G * mean_G        # local variance of guide
    cov_GI  = corr_GI - mean_G * mean_I              # local covariance
    a = cov_GI / (var_G + eps)                       # edge -> a~1, flat -> a~0
    b = mean_I - a * mean_G
    return box(a, r) * G + box(b, r)                 # average the models

img = cv2.imread("portrait.jpg", cv2.IMREAD_GRAYSCALE)
f = img.astype(np.float32) / 255.0
out = guided_filter(f, f, r=8, eps=0.02)             # self-guided smoothing
cv2.imwrite("portrait_guided.png", (out * 255).astype(np.uint8))

4. Choosing an Edge-Preserving Smoother Intermediate

Table 3.5.1 places this section's filters alongside their Section 3.2 ancestors. The pattern to internalize: as you move down the table, the filter consults the image content more and the fixed kernel less.

A second decision that arrives immediately in practice: what to do when the noise is too strong for one pass. Iterating the bilateral filter converges toward flat cartoon regions (sometimes desired for stylization, rarely for photography); the better-behaved route for serious denoising is the patch-based and learned methods of Chapter 7, which generalize "similar pixels" to "similar neighborhoods". A third option, anisotropic diffusion (Perona and Malik, 1990), reformulates edge-aware smoothing as a PDE that diffuses heat everywhere except across edges; it is historically important and conceptually elegant, but in production pipelines the one-shot guided filter almost always wins on predictability and speed.