Gaussian Filter

The Gaussian filter smooths an image by replacing every pixel with a weighted average of its neighbors, the weights falling off with distance as a Gaussian. It is the default blur of image processing: it suppresses noise before differentiation, removes detail too fine to survive subsampling in an image pyramid, and generates the scale spaces behind multi-scale feature detectors. A single parameter, the standard deviation $\sigma$, sets how large the removed details are.

Problem

Derivatives amplify pixel noise, subsampling turns fine detail into aliases, and detectors meant for large structures are distracted by texture. The Gaussian filter removes variation below a chosen spatial scale while leaving coarser structure largely intact, with a linear, shift-invariant operation that is cheap to compute and easy to analyze.

Inputs and Outputs

• Input: a grayscale or multichannel image, the standard deviation $\sigma$ in pixels (optionally one per axis), a truncation radius, and a rule for handling the border.
• Output: an image of the same size, each channel filtered independently.

Intuition

Picture each pixel exchanging some of its value with its neighbors: a lot with the nearest, less with those farther away, and almost nothing beyond a few $\sigma$. Isolated specks, which differ from everything around them, are pulled toward the local average; broad regions, where neighbors already agree, barely change; edges soften, because pixels on either side are averaged together.

Algorithm

1. Choose $\sigma$ and a radius $r = \lceil k\sigma \rceil$, typically with $k = 3$ or $4$.
2. Sample the 1D Gaussian at the integers $-r, \dots, r$ and divide by the sum of the samples, so the weights sum to one.
3. Convolve every row of the image with the 1D kernel, using the chosen border extension (reflect, replicate, wrap, or zero).
4. Convolve every column of the result with the same kernel.
5. For multichannel images, repeat for each channel.

Mathematical Formulation

The kernel

The 1D Gaussian kernel with standard deviation $\sigma$ is

$$g_\sigma(x) = \frac{1}{\sqrt{2\pi}\,\sigma} \exp\!\left(-\frac{x^2}{2\sigma^2}\right),$$

the density of a zero-mean Gaussian distribution, used here as a weighting function rather than as a model of anything random. The isotropic 2D kernel is

$$G_\sigma(x, y) = \frac{1}{2\pi\sigma^2} \exp\!\left(-\frac{x^2 + y^2}{2\sigma^2}\right) = g_\sigma(x)\, g_\sigma(y),$$

and the filtered image is the convolution $L = G_\sigma * f$.

Sampling, truncation, and normalization

The continuous kernel has infinite extent, so a digital filter samples it at integer offsets and truncates it at a radius $r$. In one dimension about 0.27% of the Gaussian’s mass lies beyond $3\sigma$, and about 0.006% beyond $4\sigma$; a radius of $\lceil 3\sigma \rceil$ is a common compromise, and SciPy’s gaussian_filter truncates at $4\sigma$ by default. The truncated samples are then divided by their sum,

$$h[i] = \frac{g_\sigma(i)}{\sum_{j=-r}^{r} g_\sigma(j)}, \qquad i = -r, \dots, r,$$

so that a constant image passes through unchanged rather than slightly darkened or brightened. For $\sigma$ below about one pixel, the sampled kernel approximates the continuous one poorly: most of its weight sits on the center tap, and its variance no longer equals $\sigma^2$.

Separability

Because $G_\sigma(x, y) = g_\sigma(x)\, g_\sigma(y)$, the 2D filter factors into a 1D convolution along rows and one along columns, which reduces the cost per pixel from $(2r+1)^2$ to $2(2r+1)$ multiplications (convolution explains why separable kernels factor this way). Up to a constant factor, the Gaussian is the only continuous, integrable 2D function that is both rotationally symmetric and a product of a function of $x$ and a function of $y$.

Frequency response

The Fourier transform of a Gaussian is a Gaussian. With frequency $u$ in cycles per pixel,

$$\hat{g}_\sigma(u) = \exp\!\left(-2\pi^2 \sigma^2 u^2\right),$$

a low-pass response that equals one at zero frequency and decreases smoothly and monotonically, with standard deviation $1/(2\pi\sigma)$. A wide kernel in space is narrow in frequency, and the Gaussian minimizes the product of the two spreads (see Fourier transform).

The smooth response means no ringing. An ideal low-pass filter has a sinc kernel whose negative lobes produce overshoot and ripples around edges. The Gaussian kernel is positive everywhere, so a filtered step stays monotonic and within the input’s range. The box filter, also positive, avoids overshoot, but the sidelobes of its sinc-shaped frequency response pass some fine detail with reversed contrast.

Cascades

Convolving two Gaussians gives a Gaussian whose variance is the sum of the two:

$$g_{\sigma_1} * g_{\sigma_2} = g_{\sqrt{\sigma_1^2 + \sigma_2^2}}.$$

This is the rule that variances of independent Gaussian variables add, seen through their densities: blurring with $\sigma_1$ and then $\sigma_2$ equals one blur with $\sqrt{\sigma_1^2 + \sigma_2^2}$. To go from an image already blurred at $\sigma_1$ to one blurred at $\sigma_2 > \sigma_1$, filter with $\sqrt{\sigma_2^2 - \sigma_1^2}$. Scale-space pyramids use this to compute each level incrementally from the previous one.

Parameters

• Standard deviation $\sigma$. The scale of the smoothing, in pixels: structures much smaller than $\sigma$ are averaged away, much larger ones survive. For denoising it trades noise reduction against detail; before subsampling by a factor $s$ it must grow roughly in proportion to $s$ (sampling and aliasing works through $\sigma = s/2$).
• Truncation radius. $3\sigma$ leaves errors of the order of $10^{-3}$ of the signal range, $4\sigma$ of the order of $10^{-5}$ (see Implementation). Where a library takes a kernel size and $\sigma$ separately, a size small relative to $\sigma$ silently truncates the Gaussian toward a box.
• Border mode. As for any convolution, reflection or replication avoids the dark border that zero padding produces.

Complexity

For an image of $N$ pixels and radius $r$:

• Direct 2D convolution: $O(N r^2)$.
• Separable convolution: $O(N r)$, the usual choice for small and moderate $\sigma$.
• FFT-based convolution: $O(N \log N)$ regardless of $\sigma$; it treats the image as periodic.
• Recursive (IIR) filtering: $O(N)$, independent of $\sigma$. Deriche’s recursive filters (from 1987) approximated the Gaussian by a sum of a causal and an anti-causal recursive filter. Young and van Vliet (1995) instead used a cascade of the two, run forward and then backward along each row and column; an order-$N$ filter costs $2N$ multiplications per pixel per dimension, six for the third-order filter they recommend, whatever $\sigma$ is. Recursive filters are approximations, least accurate for small $\sigma$, and need care at the borders.
• Repeated box filters: a box filter costs $O(N)$ with running sums, whatever its width. Repeated box filtering converges to a Gaussian by the central limit theorem, and Wells (1986) synthesized Gaussian filters from cascades of uniform filters. Three passes of width $w$ give variance $3(w^2 - 1)/12$, so integer widths reach only certain $\sigma$.

Implementation

The example builds a sampled kernel, checks separability, verifies that variances add in a cascade, measures the error caused by truncating at $k\sigma$, and compares three box filters with a Gaussian of the same variance.

import numpy as np
from scipy import ndimage, signal


def gaussian_kernel(sigma, truncate=3.0):
    """Sampled, normalized 1D Gaussian with radius ceil(truncate * sigma)."""
    r = int(np.ceil(truncate * sigma))
    x = np.arange(-r, r + 1)
    g = np.exp(-x**2 / (2 * sigma**2))
    return g / g.sum()


# 1. A sampled kernel: radius 3 sigma, normalized to sum 1.
g = gaussian_kernel(1.0)
print("sigma=1 kernel:", np.round(g, 4), f"sum = {g.sum():.6f}")

# 2. Separable filtering equals filtering with the 2D kernel.
rng = np.random.default_rng(0)
img = rng.random((128, 128))
g = gaussian_kernel(2.0)
G2 = np.outer(g, g)
full = signal.convolve2d(img, G2, mode="valid")
sep = signal.convolve2d(signal.convolve2d(img, g[None, :], mode="valid"),
                        g[:, None], mode="valid")
print(f"separable vs 2D: max difference {np.abs(full - sep).max():.1e}, "
      f"multiplications per pixel {2 * g.size} vs {G2.size}")

# 3. Cascade: blurring with sigma1 then sigma2 equals one blur with
#    sqrt(sigma1^2 + sigma2^2).
g1, g2 = gaussian_kernel(2.0, 5), gaussian_kernel(1.5, 5)
cascade = np.convolve(g1, g2)                       # combined kernel
direct = gaussian_kernel(2.5, 5)                    # sqrt(4 + 2.25) = 2.5


def variance(k):
    x = np.arange(k.size) - k.size // 2
    return np.sum(x**2 * k)


r = (cascade.size - direct.size) // 2
print(f"kernel variances: {variance(g1):.4f} + {variance(g2):.4f} = "
      f"{variance(cascade):.4f} (2.5^2 = 6.25)")
print(f"max difference from the sigma 2.5 kernel: "
      f"{np.abs(cascade[r:-r] - direct).max():.1e}")

# 4. Truncation: the mass a radius-k*sigma kernel leaves out, and the
#    resulting error on an image, relative to a wide (6 sigma) reference.
sigma = 3.0
xx = np.arange(-60, 61)
dense = np.exp(-xx**2 / (2 * sigma**2)); dense /= dense.sum()
ref = ndimage.convolve1d(img, gaussian_kernel(sigma, 6), axis=1, mode="reflect")
for k in [1, 2, 3, 4]:
    lost = dense[np.abs(xx) > np.ceil(k * sigma)].sum()
    out = ndimage.convolve1d(img, gaussian_kernel(sigma, k), axis=1, mode="reflect")
    print(f"radius {k} sigma: mass outside {lost:.2e}, "
          f"max error on image {np.abs(out - ref).max():.1e}")

# 5. Three passes of a 5-tap box filter approximate a Gaussian with the same
#    variance: 3 * (5^2 - 1) / 12 = 6, i.e. sigma = sqrt(6).
box = np.ones(5) / 5
b3 = np.convolve(np.convolve(box, box), box)        # 13 taps
gb = gaussian_kernel(np.sqrt(6), 6)
c = gb.size // 2
print(f"three 5-tap boxes vs Gaussian sigma={np.sqrt(6):.3f}: "
      f"max difference {np.abs(b3 - gb[c - 6:c + 7]).max():.4f} "
      f"(peak value {gb[c]:.4f})")

Output:

sigma=1 kernel: [0.0044 0.054  0.242  0.3991 0.242  0.054  0.0044] sum = 1.000000
separable vs 2D: max difference 5.6e-16, multiplications per pixel 26 vs 169
kernel variances: 4.0000 + 2.2500 = 6.2500 (2.5^2 = 6.25)
max difference from the sigma 2.5 kernel: 1.5e-08
radius 1 sigma: mass outside 2.41e-01, max error on image 1.4e-01
radius 2 sigma: mass outside 2.95e-02, max error on image 1.6e-02
radius 3 sigma: mass outside 1.47e-03, max error on image 9.1e-04
radius 4 sigma: mass outside 2.84e-05, max error on image 2.0e-05
three 5-tap boxes vs Gaussian sigma=2.449: max difference 0.0109 (peak value 0.1629)

The separable filter matches the 2D one to rounding error at less than a sixth of the cost. The variances of the sampled kernels add, and the cascade equals a single $\sigma = 2.5$ kernel to about $10^{-8}$. Truncating at $\sigma$ or $2\sigma$ visibly changes the result, while $3\sigma$ keeps the error below $10^{-3}$ on an image with values in $[0, 1]$; the mass left out is below the continuous figure because the radius is rounded up. Three box passes come within about 7% of the Gaussian’s peak value.

Properties and Behavior

• Linear and shift-invariant. The filter is a convolution, so it commutes with other linear filters. In particular, differentiating a Gaussian-smoothed image equals convolving with the derivative of the Gaussian, $\partial_x (G_\sigma * f) = (\partial_x G_\sigma) * f$, which is how image gradients are computed at a chosen scale.
• Unique in scale-space theory. Koenderink (1984) required that a family of progressively blurred images should not create new structure as the scale increases, a property he called causality, and combined it with the requirements that all positions and scales be treated alike. He showed that the family must then satisfy the diffusion (heat) equation, $\partial_t L = \tfrac{1}{2} \nabla^2 L$ with $t = \sigma^2$, whose solution is convolution with a Gaussian. Lindeberg (1994) and others derived the same uniqueness from other combinations of axioms, such as linearity, shift invariance, the cascade property, and non-enhancement of local extrema. For discrete signals, Lindeberg’s theory leads to a discrete analogue of the Gaussian kernel rather than to samples of the continuous one; the two differ appreciably only for small $\sigma$.
• A choice of scale, not a model of noise. Gaussian smoothing does not assume that the noise is Gaussian; it assumes that the signal of interest varies more slowly than the noise.

Limitations

• Blurs edges. The filter averages across intensity discontinuities just as it averages within smooth regions, so edges, thin lines, and corners lose contrast and spread out. Edge-preserving filters such as the bilateral filter and the median filter avoid this by making the weights depend on image content, at the cost of linearity.
• Poor against impulse noise. An outlier pixel, such as a salt-and-pepper speck, is spread over a neighborhood rather than removed; rank filters handle it better.
• One scale everywhere. A single $\sigma$ cannot preserve fine detail in one region and smooth heavily in another.

Variants

• Anisotropic Gaussian: a different $\sigma$ per axis, or a general covariance $\Sigma$ in place of $\sigma^2 I$ for oriented smoothing.
• Derivatives of Gaussian: kernels $\partial_x G_\sigma$, $\partial_y G_\sigma$, and higher derivatives, which smooth and differentiate in one step; the Canny edge detector is built on them, and the Sobel operator can be seen as a small approximation.
• Laplacian of Gaussian and difference of Gaussians: band-pass filters that respond to blobs and edges at a chosen scale (see the Laplacian operator).
• Binomial kernels: rows of Pascal’s triangle such as $\tfrac{1}{16}[1, 4, 6, 4, 1]$, integer approximations of the Gaussian used in pyramids.

Uses in Computer Vision

• Denoising and pre-smoothing. It is the baseline denoiser against which edge-preserving methods are compared, and gradients, corner measures, and edge detectors use Gaussian-smoothed derivatives, with $\sigma$ selecting the scale of the structures they respond to.
• Anti-aliasing and scale spaces. It prefilters every level of a Gaussian image pyramid, and blurring at a sequence of $\sigma$ values, computed incrementally with the cascade property, gives the scale space in which detectors such as SIFT look for features.
• Sharpening by unsharp masking. Subtracting a Gaussian-blurred copy isolates fine detail, and adding a multiple of it back, $f + \alpha\,(f - G_\sigma * f)$, sharpens the image; $\sigma$ sets which details are boosted.
• Weighting windows. Gaussian windows weight the sums in local operations such as the structure tensor and window-based motion estimation.