Bilateral Filter

The bilateral filter smooths an image while keeping its strong edges sharp. Like a Gaussian filter, it replaces each pixel with a weighted average of its neighbors, but a neighbor’s weight depends on how similar its value is as well as on how close it is, so pixels across an edge receive almost no weight. Simple and controlled by two intuitive parameters, it became a standard tool for denoising, tone mapping, detail manipulation, and edge-aware upsampling.

The name comes from Tomasi and Manduchi (1998). As the survey by Paris, Kornprobst, Tumblin and Durand (2009) recounts, the same idea appeared earlier in the nonlinear Gaussian filters of Aurich and Weule (1995) and in the SUSAN framework of Smith and Brady (1997). This article describes the method as it is used today.

Problem

Linear smoothing assumes that nearby pixels can be averaged, which fails at edges. The goal is to average within regions of similar intensity, removing noise and fine texture, but not across the boundaries between them.

Inputs and Outputs

Input: a grayscale or color image, a spatial standard deviation $\sigma_s$ in pixels, and a range standard deviation $\sigma_r$ in intensity units. A joint variant also takes a second guidance image that defines the edges.
Output: a smoothed image of the same size. Decomposition applications also use the residual (input minus output), which holds the removed detail and noise.

Intuition

Consider a pixel just on the bright side of a step edge. A Gaussian filter averages it with dark neighbors too, and the edge softens. The bilateral filter asks each neighbor how far away it is and how different its value is. Bright neighbors pass both tests; dark ones fail the second and are effectively ignored. Each side of the edge is smoothed separately, and the step survives. In a flat region every neighbor passes, and the filter acts as an ordinary Gaussian blur: its effective kernel adapts to the image, a full disk in flat regions, a half-disk beside a straight edge.

Algorithm

For each pixel $p$ :

For each neighbor $q$ within a window of radius about $2\sigma_s$ to $3\sigma_s$ , compute a spatial weight from the distance $\|p - q\|$ and a range weight from the difference $|I_p - I_q|$ .
Multiply the two weights.
Accumulate the weighted sum of neighbor values and the sum of the weights.
Divide the weighted sum by the sum of weights.

All outputs are computed from the original input in a single, non-iterative pass, so the parameters have a direct meaning.

Mathematical Formulation

The bilateral filter of a grayscale image $I$ is

\mathrm{BF}[I]_p = \frac{1}{W_p} \sum_{q \in \mathcal{S}} G_{\sigma_s}\!\left(\|p - q\|\right)\, G_{\sigma_r}\!\left(|I_p - I_q|\right)\, I_q ,

with the normalization

W_p = \sum_{q \in \mathcal{S}} G_{\sigma_s}\!\left(\|p - q\|\right)\, G_{\sigma_r}\!\left(|I_p - I_q|\right),

where $G_\sigma(x) = \exp\!\left(-x^2 / 2\sigma^2\right)$ and $\mathcal{S}$ is the set of pixel positions. $G_{\sigma_s}$ is the spatial (domain) kernel and $G_{\sigma_r}$ the range kernel. Because the range weights depend on the image, $W_p$ differs from pixel to pixel, and dividing by it makes the weights at every pixel sum to one; a constant image is left unchanged.

Tomasi and Manduchi defined the filter with general closeness and similarity functions, with Gaussians as the main example. For color images they measured similarity as Euclidean distance in CIE-Lab, a color space in which such distances roughly track perceived differences. Filtering the color vector jointly averages only perceptually similar colors, whereas filtering red, green, and blue separately can create colors at edges that were not in the image.

Two limits show what each kernel contributes. As $\sigma_r \to \infty$ , the range weights become constant and the filter reduces to Gaussian smoothing. As $\sigma_s \to \infty$ , the spatial weights become constant and the filter becomes a pure range filter; Tomasi and Manduchi showed that this merely remaps gray levels according to the image histogram, without any spatial smoothing. Edge-preserving smoothing needs both.

Parameters

Spatial $\sigma_s$ . The neighborhood size, and so the scale of the features removed; it must scale with image resolution. Denoising typically uses small values, decomposing an image into a smooth base and a detail layer large ones.
Range $\sigma_r$ . The contrast of the edges preserved: differences much smaller than $\sigma_r$ are averaged, much larger ones are not. It must scale with the intensity range ( $[0, 1]$ versus $[0, 255]$ ). For denoising it is set in proportion to the noise level; Paris et al. report a recommendation by Liu and colleagues of $\sigma_r \approx 1.95\,\sigma_n$ for a local noise standard deviation $\sigma_n$ .

Because the weights multiply, no averaging occurs where either is near zero: increasing $\sigma_s$ does not blur an edge whose contrast is well above $\sigma_r$ .

Complexity

A brute-force implementation over a window of radius $r \propto \sigma_s$ costs $O(N r^2)$ for $N$ pixels. The filter is neither a convolution nor separable, so neither the FFT nor row-and-column factoring applies directly. Tomasi and Manduchi precomputed the range weights in a lookup table, one entry per intensity difference; even so, Paris et al. report that a direct implementation can take minutes for a megapixel image with a large kernel.

Fast approximations exploit the observation that the filter becomes linear in a higher-dimensional space:

Piecewise-linear approximation. Durand and Dorsey (2002) fixed a few intensity levels, computed for each a Gaussian convolution of the image weighted by the range kernel centered at that level, and interpolated between the two levels nearest each pixel’s value, on a downsampled image. This cut the time for a megapixel image to about a second.
Bilateral grid. Paris and Durand (2006) showed that bilateral filtering a grayscale image is exactly a 3D Gaussian convolution over position and intensity of a homogeneous representation of the image, followed by sampling at each pixel’s own intensity (“slicing”) and a division. Because the convolved volume is band-limited, it can be computed on a coarse 3D grid, which they named the bilateral grid. Chen, Paris and Durand (2007) made the grid a primary data structure, ran it on graphics hardware in real time for high-definition video, and used it for other edge-aware operations such as local histogram equalization and edge-aware painting.
Permutohedral lattice. Color bilateral filtering is a 5D problem (two positions, three colors), where regular grids grow expensive. Adams, Baek and Davis (2010) sampled the space on the permutohedral lattice, whose simplices make splatting and slicing cheap, and obtained a high-dimensional Gaussian filter whose cost is linear in the number of pixels and polynomial ( $O(d^2)$ per pixel) in the dimension $d$ .

Implementation

The example filters a noisy vertical step edge (height 0.6) by brute force and compares it with a Gaussian filter of the same $\sigma_s = 2$ , at a low and a high noise level. It reports the noise remaining in a flat region, the jump between the two columns either side of the edge, and the overall error against the clean image.

import numpy as np
from scipy import ndimage


def bilateral(img, sigma_s, sigma_r, truncate=3.0):
    """Brute-force bilateral filter of a grayscale image (reflect padding)."""
    r = int(np.ceil(truncate * sigma_s))
    pad = np.pad(img, r, mode="reflect")
    num = np.zeros_like(img)
    den = np.zeros_like(img)
    h, w = img.shape
    for dy in range(-r, r + 1):
        for dx in range(-r, r + 1):
            spatial = np.exp(-(dx**2 + dy**2) / (2 * sigma_s**2))
            neighbor = pad[r + dy:r + dy + h, r + dx:r + dx + w]
            weight = spatial * np.exp(-(neighbor - img)**2 / (2 * sigma_r**2))
            num += weight * neighbor
            den += weight
    return num / den


def report(name, out, clean):
    flat = out[:, 8:24]                       # left region, away from the edge
    contrast = out[:, 32].mean() - out[:, 31].mean()   # jump across the edge
    rmse = np.sqrt(np.mean((out - clean)**2))
    print(f"{name:26s} flat-region std {flat.std():.4f}   "
          f"edge jump {contrast:.3f}   RMSE {rmse:.4f}")


rng = np.random.default_rng(0)
clean = np.full((64, 64), 0.2)
clean[:, 32:] = 0.8                           # vertical step edge of height 0.6
for sigma_n in [0.05, 0.3]:
    noisy = clean + rng.normal(0, sigma_n, clean.shape)
    print(f"noise sigma {sigma_n}:")
    report("input", noisy, clean)
    report("Gaussian, sigma 2", ndimage.gaussian_filter(noisy, 2, mode="reflect",
                                                        truncate=3.0), clean)
    for sigma_r in [0.1, 0.3, 10.0]:
        report(f"bilateral, sigma_r {sigma_r}", bilateral(noisy, 2, sigma_r), clean)

Output:

noise sigma 0.05:
input                      flat-region std 0.0502   edge jump 0.591   RMSE 0.0499
Gaussian, sigma 2          flat-region std 0.0080   edge jump 0.119   RMSE 0.0507
bilateral, sigma_r 0.1     flat-region std 0.0143   edge jump 0.595   RMSE 0.0141
bilateral, sigma_r 0.3     flat-region std 0.0083   edge jump 0.489   RMSE 0.0142
bilateral, sigma_r 10.0    flat-region std 0.0080   edge jump 0.120   RMSE 0.0507
noise sigma 0.3:
input                      flat-region std 0.3052   edge jump 0.614   RMSE 0.3012
Gaussian, sigma 2          flat-region std 0.0479   edge jump 0.120   RMSE 0.0693
bilateral, sigma_r 0.1     flat-region std 0.2866   edge jump 0.601   RMSE 0.2827
bilateral, sigma_r 0.3     flat-region std 0.1745   edge jump 0.481   RMSE 0.1749
bilateral, sigma_r 10.0    flat-region std 0.0473   edge jump 0.121   RMSE 0.0689

At low noise, the Gaussian filter reduces flat-region noise sixfold but shrinks the edge jump from 0.6 to 0.12, so its overall error is no better than the input’s. The bilateral filter with $\sigma_r = 0.1$ , twice the noise level, keeps the edge and cuts the error by more than a factor of three; it leaves more flat-region noise because the range kernel also down-weights strongly noisy neighbors. With $\sigma_r = 0.3$ it smooths flat regions as well as the Gaussian but softens the edge, and with $\sigma_r = 10$ it is the Gaussian. At high noise, half the edge height, no $\sigma_r$ separates noise from edge: a narrow range kernel treats the noise as edges and barely smooths, and the plain Gaussian has the lowest error. The range kernel works only when edges are much stronger than the noise.

Properties and Behavior

Nonlinear but shift-invariant. The weights depend on the image, so filtering a sum of images is not the sum of the filtered images. Shifting the image, or adding a constant to its intensities, shifts the output accordingly.
Staircasing. The filter is best at step-like edges. It tends to break smooth gradients into piecewise-constant pieces, the staircase effect analyzed by Buades, Coll and Morel, and can sharpen soft edges into steps.
Gradient reversal. When the residual detail layer is amplified, as in detail enhancement or tone mapping, the filter can produce reversed profiles next to edges. He, Sun and Tang (2013) trace this to edge pixels with few similar neighbors, whose weighted average becomes unstable.
Impulse noise survives. An outlier differs from all its neighbors, so the range kernel gives them little weight and the outlier stays. Paris et al. describe computing the range weights from a median-filtered copy of the image to handle this case.

Limitations

The bilateral filter is a simple, flexible denoiser rather than a state-of-the-art one: it compares single pixel values, which are themselves noisy. Its brute-force cost grows with $\sigma_s^2$ , and fast versions are approximations whose accuracy and memory depend on $\sigma_s$ , $\sigma_r$ , and the range dimension.

Variants

Joint (cross) bilateral filter. Petschnigg and colleagues (2004) and Eisemann and Durand (2004) independently computed the range weights from a second image: a flash photograph supplies reliable edges for smoothing a noisy no-flash photograph of the same scene. Kopf and colleagues (2007) used it for joint bilateral upsampling: a high-resolution image guides the upsampling of a solution computed at low resolution.
Guided filter. He, Sun and Tang (2010, 2013) proposed an alternative based on a local linear model of the output in terms of a guidance image. It preserves edges in a similar way, avoids gradient reversal, and has an exact algorithm whose cost is linear in the number of pixels, regardless of kernel size and intensity range.
Relation to anisotropic diffusion. Perona–Malik diffusion also smooths less across strong gradients, but iteratively, as a partial differential equation. Durand and Dorsey showed that a bilateral filter restricted to the four nearest neighbors is a discrete Perona–Malik scheme.
Neighborhood and patch-based filters. Yaroslavsky’s neighborhood filter is the bilateral filter with a box spatial kernel. Non-local means compares whole patches instead of single pixel values, which makes the weights more robust to noise; Adams and colleagues treat it and the bilateral filter alike as high-dimensional Gaussian filters.

Uses in Computer Vision

Tone mapping. Durand and Dorsey (2002) split the log intensity of a high-dynamic-range image into a bilateral-filtered base layer and a residual detail layer. Compressing only the base and adding the detail back reduces contrast without the halos a Gaussian base layer produces at high-contrast edges.
Detail enhancement, stylization, and guided upsampling. Amplifying the residual exaggerates texture, iterating the filter drives the image toward piecewise-constant, cartoon-like regions, and joint bilateral filtering upsamples low-resolution results such as depth maps from range sensors.
Segmentation refinement. Krähenbühl and Koltun (2011) built a fully connected conditional random field whose pairwise term includes a bilateral appearance kernel over position and color, with mean-field inference computed by permutohedral-lattice filtering. The first DeepLab model (Chen and colleagues, 2014) applied such a CRF to a convolutional network’s output to sharpen segment boundaries.