Gaussian and Median Filtering

A digital image $I(x,y)$ is a discrete sampling of a scene corrupted by acquisition noise, quantization, and sometimes transmission artifacts. Smoothing (or denoising) seeks an estimate $\hat{I}$ of the true underlying signal that suppresses noise while retaining meaningful structure such as edges and textures. This is fundamentally a local estimation problem: at each pixel, combine the values in a neighborhood into a single output value, trusting that nearby pixels carry correlated information about the true scene. The Gaussian filter performs this combination via Linear Filtering and Convolution with a kernel derived from the 2D Gaussian density,

The output is the convolution $\hat{I} = I * G_\sigma$, a weighted average where nearer pixels receive exponentially larger weight. Because the kernel is a genuine probability density (it integrates to 1, after discretization and renormalization), the filter cannot create new intensity extrema; it can only average existing values, which is precisely why it is a low-pass filter — it attenuates high spatial frequencies (fine detail, noise) while passing low frequencies (smooth regions) largely unchanged. The median filter replaces the linear weighted average with an order statistic: at each pixel, sort the intensities inside a window (typically $k \times k$, $k$ odd) and output the middle value. Because a median is a robust statistic — it is insensitive to any single extreme value in the window, unlike a mean — the median filter can remove impulsive noise (isolated corrupted pixels with arbitrary intensity, so-called "salt-and-pepper" noise) while leaving step edges essentially untouched, since one side of an edge always dominates the window and supplies the median. This nonlinearity means the median filter has no convolution representation, no frequency response in the usual sense, and no superposition property — but it buys edge preservation that no linear filter can match.

Gaussian filtering. The 2D Gaussian kernel is separable: $G_\sigma(x,y) = g_\sigma(x)\, g_\sigma(y)$ where $g_\sigma$ is the 1D Gaussian. This algebraic fact follows directly from the exponential splitting $e^{-(x^2+y^2)/2\sigma^2} = e^{-x^2/2\sigma^2}\cdot e^{-y^2/2\sigma^2}$. Consequently, a 2D convolution can be computed as two sequential 1D convolutions — one along rows, one along columns — which is both simpler to implement and dramatically faster. The kernel is truncated to a finite support, typically a half-width of $r = \lceil 3\sigma \rceil$ (capturing more than 99.7% of the Gaussian's mass), and normalized so its discrete entries sum to 1 (preserving average brightness). ``` function GAUSSIAN_KERNEL_1D(sigma): r = ceil(3 sigma) k = array of size (2r + 1) for i in -r .. r: k[i + r] = exp(-(i i) / (2 sigma sigma)) normalize k so sum(k) == 1 return k function CONVOLVE_1D_HORIZONTAL(I, k, r):

# I: H x W image, k: kernel of half-width r O = array H x W, zero-initialized for y in 0 .. H-1: for x in 0 .. W-1: acc = 0 for i in -r .. r: xi = clamp(x + i, 0, W - 1) # border handling acc += I[y][xi] * k[i + r] O[y][x] = acc return O function CONVOLVE_1D_VERTICAL(I, k, r):

# symmetric to the horizontal pass, sums over rows O = array H x W, zero-initialized for y in 0 .. H-1: for x in 0 .. W-1: acc = 0 for i in -r .. r: yi = clamp(y + i, 0, H - 1) acc += I[yi][x] * k[i + r] O[y][x] = acc return O function GAUSSIAN_BLUR(I, sigma): k = GAUSSIAN_KERNEL_1D(sigma) r = (length(k) - 1) / 2 T = CONVOLVE_1D_HORIZONTAL(I, k, r) O = CONVOLVE_1D_VERTICAL(T, k, r) return O ``` Border handling (`clamp`, reflect, or wrap) is an implementation choice; clamping (replicate edge pixels) is simplest and avoids introducing artificial dark/bright borders.

Median filtering. The naive algorithm collects the $k^2$ neighborhood values for every pixel, sorts them, and picks the middle element. ``` function MEDIAN_FILTER(I, k):

# k odd, window is k x k, r = (k-1)/2 r = (k - 1) / 2 O = array H x W for y in 0 .. H-1: for x in 0 .. W-1: window = [] for dy in -r .. r: for dx in -r .. r: yi = clamp(y + dy, 0, H - 1) xi = clamp(x + dx, 0, W - 1) window.append(I[yi][xi]) sort(window) O[y][x] = window[length(window) / 2] # middle element return O ``` For a fixed small window (e.g. $3\times3$ or $5\times5$), this is fast enough in practice. For large windows, the classic Huang running-median algorithm exploits the fact that as the window slides one column at a time, only one column of pixels enters and one leaves; a per-column histogram (256 bins for 8-bit images) is updated incrementally, and the median is tracked via a running count rather than resorting from scratch. ``` function HUANG_MEDIAN_FILTER(I, k): r = (k - 1) / 2 O = array H x W for y in 0 .. H-1: hist = zeros(256) for dy in -r .. r: for dx in -r .. r: yi = clamp(y + dy, 0, H - 1) hist[I[yi][clamp(dx, 0, W-1)]] += 1 median_idx = FIND_MEDIAN(hist, kk) O[y][0] = median_idx for x in 1 .. W-1: for dy in -r .. r: yi = clamp(y + dy, 0, H - 1) out_col = clamp(x - r - 1, 0, W - 1) in_col = clamp(x + r, 0, W - 1) hist[I[yi][out_col]] -= 1 hist[I[yi][in_col]] += 1 median_idx = FIND_MEDIAN(hist, kk) # incremental scan, amortized O(1) O[y][x] = median_idx return O ``` `FIND_MEDIAN` tracks a running cumulative count and a pointer into the histogram; because the histogram changes by only $2k$ entries per column step, the median pointer moves only a small, bounded amount on average, giving an amortized near-constant update per pixel.

Gaussian filter — time and space. For an $H \times W$ image and a kernel with half-width $r$ (so kernel diameter $n = 2r+1$):

• Direct 2D convolution: each output pixel costs $O(n^2)$ multiply-adds, so total time is $O(H W n^2)$.
• Separable implementation: each output pixel of the horizontal pass costs $O(n)$, and likewise for the vertical pass, giving total time $O(H W n)$ — a factor of $n$ speedup, which matters enormously since $n = O(\sigma)$ and large $\sigma$ (needed for heavy denoising) makes $n^2$ prohibitive.
• Space: $O(HW)$ for the image buffers, plus $O(n)$ for the 1D kernel; the intermediate buffer $T$ adds another $O(HW)$, so total space is $O(HW)$ (constant-factor overhead).

Why separability is correct, not just faster: 2D convolution is associative and the 2D Gaussian factors as an outer product $G_\sigma(x,y) = g_\sigma(x) g_\sigma(y)$. Substituting into the convolution integral,

which is exactly "convolve each row with $g_\sigma$, then convolve each resulting column with $g_\sigma$." The inner integral is the horizontal pass; the outer is the vertical pass. No approximation is involved — the separable computation produces the exact same result as direct 2D convolution (up to floating-point rounding), not merely a fast approximation. Why Gaussian smoothing is the "right" linear filter: under the model that each pixel is the true signal plus i.i.d. Gaussian noise, a local weighted average is the maximum-likelihood estimate of the mean, and the Gaussian weights are the unique kernel that is simultaneously separable, rotationally symmetric, and infinitely smooth — properties needed for consistent multi-scale (scale-space) analysis, since Gaussian blurring is the unique linear, shift-invariant smoothing that does not create new local extrema as $\sigma$ increases (the scale-space causality property). This is why the Gaussian, rather than a box or triangular filter, underlies the construction of Image Pyramids.

Median filter — time and space. For a $k \times k$ window ($n = k^2$ samples per pixel):

• Naive (collect + sort per pixel): each pixel costs $O(n \log n)$ for sorting, so total time is $O(HW k^2 \log(k^2)) = O(HW k^2 \log k)$.
• Selection-based (median-of-medians or quickselect instead of full sort): $O(n)$ expected (or worst-case with median-of-medians) per pixel, giving $O(HW k^2)$ total.
• Huang's running-histogram method: maintaining a per-column histogram and updating it incrementally as the window slides costs $O(k)$ to add/remove a column and amortized small cost to relocate the median pointer, giving $O(HW k)$ total for 8-bit images (histogram bins bounded by a constant, 256) — a substantial improvement for large $k$.
• Space: $O(HW)$ for image buffers; $O(256)$ or $O(\text{bit-depth range})$ for the histogram in Huang's method, versus $O(k^2)$ for the naive per-pixel window buffer.

Why the median is correct (robust) where the mean is not: the breakdown point of an estimator is the fraction of arbitrarily corrupted samples it can tolerate before producing an arbitrarily wrong output. The mean has breakdown point $0$: a single sample driven to $+\infty$ drags the mean to $+\infty$. The median has breakdown point $\tfrac{1}{2}$: up to (just under) half the samples in the window can be arbitrarily corrupted and the median remains bounded by the true values, because the median is defined purely by rank, not magnitude. This is exactly why median filtering removes salt-and-pepper noise (where a small fraction of pixels take extreme values) essentially perfectly for reasonable noise densities, while a mean/Gaussian filter only attenuates such noise proportionally to kernel weight. Edge preservation, argued precisely: consider a 1D step edge with true values $a$ on the left and $b > a$ on the right, and a window of odd length $k$ centered at position $x$. If more than half the window lies on one side of the edge, the median equals exactly $a$ or exactly $b$ — the true value — with zero blur, regardless of $k$. A Gaussian filter, in contrast, always produces a weighted blend of $a$ and $b$ whenever any nonzero-weight kernel taps straddle the edge, producing a smeared transition whose width scales with $\sigma$. This is the formal sense in which median filtering is "edge-preserving" and Gaussian filtering is not.

Honest tradeoffs. The Gaussian filter is linear (compositions of Gaussian blurs are themselves Gaussian, $G_{\sigma_1} * G_{\sigma_2} = G_{\sqrt{\sigma_1^2+\sigma_2^2}}$), has a clean frequency response (a Gaussian in the spatial domain transforms to a Gaussian in frequency, so it is a smooth low-pass filter with no ringing), and is separable and fast — but it blurs every edge and cannot remove impulsive noise without an impractically small kernel. The median filter preserves edges and kills outliers, but it is nonlinear (no frequency-domain shortcut, no simple composition rule), can erode thin structures (a line thinner than the window's "majority" width gets removed, since the majority of window entries then come from the background), and its fast implementation (Huang's algorithm) is more complex to code than a straightforward separable convolution.

Take a $5\times5$ 8-bit image with a clean vertical step edge (left half = 20, right half = 200), and inject one salt-noise pixel (value 255) at position $(2,2)$ (0-indexed, row, col): ``` 20 20 20 200 200 20 20 20 200 200 20 20 255 200 200 20 20 20 200 200 20 20 20 200 200 ```

Gaussian filter, $3\times3$ kernel from $\sigma = 1$ (unnormalized weights, rounded, then normalized to sum 1):

Applying $K$ centered at $(2,2)$ (using clamped borders where needed), the weighted sum over the $3\times3$ neighborhood ``` 20 20 200 20 255 200 20 20 200 ``` is

The Gaussian output at $(2,2)$ is about 124 — a value that is neither the clean background (20) nor the clean foreground (200) nor the noise spike (255); the noise has been blended into the result and the edge nearby is softened over roughly a $\pm 3\sigma$ band.

Median filter, same $3\times3$ window at $(2,2)$: sort the nine values ``` 20 20 20 20 20 20 200 200 255 ``` sorted: $20,20,20,20,20,20,200,200,255$. The median (5th of 9) is 20. The noise spike is completely eliminated, and the output at $(2,2)$ correctly reflects that the majority of the neighborhood is background. Away from the spike — say at pixel $(2,3)$, whose $3\times3$ window is entirely on the "200" side except for one column touching the edge — the median filter also reproduces a clean value near 200, whereas the Gaussian filter at the same location would show a smooth, gradual transition spanning several pixels across the true edge.

This trace demonstrates concretely: the linear Gaussian filter spreads the single bad sample's influence to its neighbors (since convolution is a weighted sum that cannot ignore any input), while the nonlinear median filter's rank-based selection discards the outlier outright.

The Gaussian kernel's mathematical roots trace to Carl Friedrich Gauss's work on the normal distribution and least-squares estimation, but its systematic use as an image-smoothing and multiscale analysis tool emerged from 20th-century signal processing and the computer vision community's development of scale-space theory, where Gaussian blurring was shown to be the unique kernel satisfying causality (no new spurious extrema) across scales. The median filter's development is more directly tied to robust statistics and digital signal processing.

Nearly every camera image-signal processor (ISP) applies denoising — often a combination of Gaussian-like and median-like filters — before demosaicing, sharpening, and compression. Photo and video editing software (Adobe Photoshop, GIMP, mobile camera pipelines) expose both "Gaussian Blur" and "Median" or "Despeckle" filters directly as user-facing tools. In computer vision pipelines, Gaussian smoothing is the mandatory precursor to gradient-based operators such as Sobel or Canny edge detection (since numerical differentiation amplifies high-frequency noise, a topic developed further in Image Gradients), and it is the core recursive operation used to build Image Pyramids for scale-invariant keypoint detectors such as SIFT and for multiresolution blending and coarse-to-fine optical flow. Median filtering is used in astronomical imaging to remove cosmic-ray hits, in medical imaging to suppress impulsive sensor artifacts, and in classic salt-and-pepper noise removal for scanned documents and old photographs. Both filters also serve as pedagogical anchors: the Gaussian illustrates the power and limits of linear, shift-invariant systems, while the median illustrates why nonlinear, order-statistic-based methods are indispensable whenever real-world noise is not well modeled by a Gaussian distribution.

Ask the tutor: How does the bilateral filter modify the Gaussian kernel to preserve edges, and why does that break separability?