10.1 Bilateral filtering🔗⧉

A blur mixes two sides of an edge and makes it soft. That is useful for some operations, but harmful when we want to remove noise while preserving sharpness or to smooth texture while keeping the boundary. The big idea of the bilateral filter is to only take into account pixels that are similar enough. This notion of affinity between pixels will come back over and over for edge-preserving techniques.

We illustrate bilateral filtering for the task of denoising, but it also underlies many algorithms because it enables the factoring of images into large-scale structure and local detail while preserving strong edges. This proves powerful for applications such as high-dynamic-range tone mapping (Chapter 18.3, Two-scale and bilateral tone mapping), upsampling low-resolution results, and combining images with different illumination, where the edges come from a second image (#Cross / joint bilateral).

10.1.1 Bilateral filter🔗⧉

Consider an input signal composed of a step edge with noise (Figure 10.1.1⛶ fullscreen). The goal of bilateral filtering is to remove the noise by averaging nearby pixels without blurring the edge.

Start from the Gaussian blur. It replaces each pixel with a weighted average of its neighbors, where the weight only depends on distance. It is a spatial Gaussian $g_s(\lVert p - q\rVert)$ of the distance between the center pixel $p$ and a neighbor $q$: nearby pixels count for a lot, distant ones taper off. It is easiest to watch this on a 1-D scanline (a single row of the image), where the picture becomes a wiggly curve of value against position (Figure 10.1.1⛶ fullscreen).

The trouble, once more, is the edge. A step edge on the scanline, with values low on the left and high on the right, gets smeared by the Gaussian, because right at the step the averaging window reaches across the boundary and blends the two plateaus together. The same bleed, in two dimensions and on a compressed base layer, is what paints halos along high-contrast edges in tone mapping (Chapter 18.3, Two-scale and bilateral tone mapping).

The fix is a single extra factor. Keep the spatial Gaussian, but multiply in a second Gaussian, this one on the value difference between the center pixel and the neighbor, $g_r(\lvert I_p - I_q\rvert)$. This range weight asks "are these two pixels similar in value?" and down-weights neighbors that are too different. A neighbor that is close in space but far in value, one on the other side of an edge, now earns a weight near zero, however spatially near it sits. The filter, all by itself, declines to average across the edge.

A practical wrinkle follows. The range weight asks whether two pixels are similar. But on a badly corrupted image the raw values are an unreliable judge: noise can make two same-region pixels look different, or two different ones look alike, so the weights themselves come out noisy. The weights need not be read from the original image at all. We can compute $g_r$ from a pre-cleaned version (a quick median filter, say, or any cheap denoise) and then apply those steadier weights when averaging the original pixels. Decoupling what decides the affinity from what actually gets averaged is exactly the move the cross / joint bilateral (below) makes with a separate guide image; using a median-filtered copy of the same image as that guide is just the in-place version of the trick.

Putting the two Gaussians together gives the bilateral filter (Tomasi & Manduchi 1998). The weight a neighbor $q$ receives at pixel $p$ is the product of the spatial and range terms,

and the output is the normalized weighted average of the neighbors,

In words: "average my neighbors, but only those that are both nearby and similar to me." The two widths do different jobs. $\sigma_s$, in pixels, sets how large a region we pool over, the reach of the spatial Gaussian. $\sigma_r$, in intensity units, sets what counts as "too different": the value gap beyond which a neighbor stops counting. Because $\sigma_r$ is measured in the same units as $I$, we must fix the encoding before we can quote it: for HDR tone mapping the filter runs on log luminance (Durand & Dorsey 2002, where a range $\sigma_r \approx 0.4$ in $\log_{10}$ is typical); other tasks may use linear values or a specified output encoding. There is no single required domain for every bilateral application.

As code, the brute-force bilateral is a double loop with a per-pixel normalization, and one subtlety distinguishes it from convolution: you may precompute the spatial weight $g_s$ (it is the same window everywhere), but the range weight $g_r$ must be recomputed for every pair, because it depends on the value difference at that pixel. The range function can be tabulated for quantized differences; it is the resulting image-dependent kernel that cannot be reused as one convolution. Specify spatial truncation and boundary handling, include the center sample, and accumulate in a sufficiently wide type.

precompute g_s over the window          # spatial weights, same for every pixel
for each output pixel p:
    sum ← 0;  W ← 0                      # accumulator and normalizer
    for each neighbour q in the window around p:
        w ← g_s(‖p − q‖) · g_r(|I[p] − I[q]|)   # range term recomputed per pair
        sum ← sum + w · I[q]
        W   ← W + w
    I_bf[p] ← sum / W                    # per-pixel normalization

For a color image the range distance $|I_p - I_q|$ is the 3-D distance in RGB (or YUV); everything else is unchanged.

The filter also travels under other names in editing software. Photoshop's Surface Blur is a bilateral filter: its radius is the spatial width and its threshold is the range width, over a box spatial kernel rather than a Gaussian. The name is worth knowing: those two sliders are the two widths of this section, and a tool many readers already use is a direct way to feel what each one does.

The 1-D picture carries straight to the 2-D image: the bilateral smooths flat texture and noise inside each region while holding the edges, so a portrait emerges with clean skin and intact eyelashes, and an HDR base layer stops cleanly at the silhouette instead of ramping across it (Figure 10.1.2⛶ fullscreen).

The range Gaussian $g_r$ is more than a clever weight; it is a similarity, an affinity between two pixels: a number that says how much these two pixels belong together, computed from their value difference. The bilateral averages pixels in proportion to that affinity. Every later method in this chapter, and several beyond it, is a different answer to one question: how should we measure whether two pixels belong together?

Debugging: hand-checkable limits🔗⧉

The bilateral has two limits you can predict by hand, and together they make an excellent correctness test for an implementation (Figure 10.1.7⛶ fullscreen). Push the range width to each extreme:

• Large range variance ($\sigma_r \to \infty$): the range Gaussian flattens, $g_r \to 1$, so it does nothing: every value difference counts as "similar." The weight reduces to the spatial Gaussian alone, and the bilateral collapses to a plain Gaussian blur.
• Tiny range variance ($\sigma_r \to 0$): only neighbors with the same value retain weight in the limit. They may include many pixels, but averaging equal values leaves the output unchanged. Averaging those equal values changes nothing, and the output is the identity: the image comes back untouched, edges and all.

A middle $\sigma_r$ interpolates between these two: enough to smooth away noise and texture, not so much that it bleeds across real edges. If your implementation does not reproduce a Gaussian at large $\sigma_r$ and the identity at small $\sigma_r$, it is wrong. These are the first two things to check (and the recipe in the debugging sidebar above).

The same two knobs guide tuning. Pick $\sigma_s$ from the spatial scale of the structure you mean to remove (how large a region should be pooled). Pick $\sigma_r$ from the contrast of the edges you mean to preserve: set it just below the value jump of a real edge, so the filter smooths everything gentler than that and stops at anything sharper. And state the encoding first: the same numeric $\sigma_r$ means something completely different in linear light, in gamma, and in log luminance.

10.1.2 Robustness🔗⧉

With the Gaussians defined above and the center included, $w(p,p)=1$, so $W_p\geq1$. The normalized bilateral is a convex average and cannot blow up from division by a near-zero normalizer. An isolated value can still be poorly denoised because it assigns little weight to other samples. This is a limitation of the affinity, not proof that normalization is numerically invalid.

The deeper weakness is that the range term reads the noisy image itself. At an outlier the center's own corrupted value decides which neighbors count as similar, so the pixel that most needs help is the one whose weights are least trustworthy: it finds nothing that resembles it and averages with almost nothing but itself. A hot pixel survives the filter for the same reason an edge does.

A cheap repair is to compute the range term from a pre-filtered copy instead of from the input. Run a small median filter first, then use that image as the reference for the value differences while still averaging the original samples. The median removes impulses without moving an edge (Chapter 13.3, Gaussian, median, and bilateral filtering), so the weights are decided by something close to the clean signal while the output remains an average of real measurements. Mechanically this is the cross-bilateral construction of the next section, with the guide derived from the input rather than captured separately.

The cost is that the guide's failures become the filter's. A median erases thin lines, small dots and isolated highlights, and a structure the guide has removed has no range weight left to defend it: the bilateral will happily average it away. Check the guide on the structures you care about before trusting the result, and keep the median's window small enough that it removes impulses rather than detail.

An alternative discussed by Aubry et al. is an unnormalized correction Aubry et al. 2014. Its scale must be defined. With these unnormalized spatial weights, introduce a step size $\eta$:

Each neighbor contributes a weighted difference from the center. A sufficient condition for this update to remain a convex combination is $0\leq\eta\sum_{q\ne p}w(p,q)\leq1$. One practical convention normalizes the spatial kernel to sum to one and uses $0\leq\eta\leq1$. Without a scale convention, adding more neighbors can cause overshoot. The update generally differs from the normalized bilateral; neither is a universal replacement for the other.

This unnormalized form is also the bilateral cousin of the local Laplacian filters: working in deltas from the input, scale by scale, is exactly the move that makes the multiscale local Laplacian fast and robust (Paris / Aubry; cross-ref Chapter 10.4, Local Laplacian filters and Chapter 6.6, Pyramids and wavelets).

10.1.3 Cross / joint bilateral🔗⧉

So far the affinity has been read off the same image we are smoothing (the range weight uses $I$ itself). But there is no reason it must. The cross (or joint) bilateral filter reads the range weight from a different image, a guide $J$, while still averaging the target $I$:

The averaging runs over the target $I$; the affinity comes from the guide $J$. Why want that? Because sometimes the image you wish to smooth is too noisy or too dark for its own edges to be trusted. If you read the affinity off $I$, you are reading off the noise. If you hold a cleaner signal that shares the same scene geometry, borrow its edges instead.

The headline application is flash / no-flash photography (Petschnigg et al. 2004; Eisemann & Durand 2004). Shoot a dim scene twice: once without flash (the right mood and color, but dark and noisy) and once with flash (clean and well-exposed, but flat, harsh lighting). After registration, many object boundaries correspond. Motion, occlusion, flash shadows and specularities can introduce mismatched edges, so a practical pipeline needs confidence masks or fallback behavior. So smooth the noisy no-flash image using the flash image as the guide: the flash image's clean edges drive the range weight, telling the filter where it may and may not average, while the values being averaged come from the ambient no-flash shot. The result keeps the no-flash mood and kills the no-flash noise (Figure 10.1.8⛶ fullscreen). We develop the full flash/no-flash pipeline later (forward-ref → Flash / no flash, computational illumination).

The full flash/no-flash story, including Petschnigg's detail transfer, Eisemann & Durand's intrinsic-relighting pipeline, flash-shadow removal, and their end-to-end results, belongs with the rest of computational illumination, where capturing the scene under two lights is the whole subject (Flash / no flash).

The cross bilateral is the first instance of a more general joint / guided idea (filter one image using structure taken from another) that culminates in the guided filter (later in the part).

10.1.4 Bilateral grid🔗⧉

The bilateral filter is non-linear: the weights depend on the pixel values, so they differ at every output pixel, and the whole thing is not a convolution. That is why you cannot reuse your fast Gaussian-blur code, and why a brute-force bilateral is a slow double loop. The bilateral grid (Chen, Paris & Durand 2007) is the observation that you can turn it back into an ordinary convolution: add a dimension (Figure 10.1.9⛶ fullscreen).

The key picture: the product of a spatial Gaussian (in $x, y$) and a range Gaussian (in $I$) is just a single 3-D Gaussian in the space $(x, y, I)$. So lift every pixel out of the 2-D plane into a 3-D space by using its value as a third coordinate. Pixel $(x, y)$ with value $I(x, y)$ goes to grid position $(x, y, I(x,y))$. In this lifted space the nonlinear 2-D bilateral becomes a plain, separable 3-D Gaussian blur, that is, an ordinary linear convolution, and the edge-preserving behavior is automatic: two pixels on opposite sides of an edge have very different values, so they land in different layers of the grid, which reduces their interaction according to the range bandwidth. A Gaussian still has nonzero tails, and discretization adds approximation error.

The pipeline has three steps (splat, blur, slice):

• Splat. Build a coarse 3-D grid indexed by $(x, y, I)$. For each pixel, accumulate two things into its cell: its value and a weight (a count of how many pixels landed there). Many cells receive no pixels and stay empty; the per-cell weight is exactly the bookkeeping that will reproduce the bilateral's normalizer.
• Blur. Convolve the grid with a plain separable 3-D Gaussian, an ordinary blur along each axis in turn (the $[1\ 4\ 6\ 4\ 1]/16$ stencil if the Gaussian is one cell wide). Blur both the accumulated values and the weights. Spatial smoothing and range smoothing both fall out of this single convolution; mixing across layers is attenuated according to the range kernel.
• Slice. Read the result back out at each pixel's own location $(x, y, I(x,y))$ by trilinear interpolation, and divide the interpolated value by the interpolated weight: that division is the bilateral normalizer $1/W_p$, recovered for free.

The whole pipeline is short, and every step is an ordinary operation (accumulate, blur, look up):

# splat — build the coarse grid (downsample by σ_s in space, σ_r in range)
for each pixel (x, y):
    c ← (x/σ_s, y/σ_s, I[x,y]/σ_r)      # the cell this pixel lands in
    splat(grid, c, I[x,y])              # distribute to surrounding vertices
    splat(weight, c, 1)                 # same interpolation weights

# blur — plain separable 3-D Gaussian, on BOTH grid and weight
blur grid, weight along x, then y, then I    # e.g. [1 4 6 4 1]/16 per axis

# slice — read back at each pixel's own (x, y, I), normalized
for each pixel (x, y):
    c ← (x/σ_s, y/σ_s, I[x,y]/σ_r)
    out[x, y] ← interpolate(grid, c) / interpolate(weight, c)      # trilinear interpolation at c
return out

Why is this fast? Because the grid is coarse. The right resolution is set by the filter's own widths (a grid cell about $\sigma_s$ wide in space and $\sigma_r$ deep in range), so the grid is heavily downsampled in both space and value, often by large factors. A small grid, a separable blur, a trilinear read: the whole thing runs in real time on a GPU. This is the pipeline that Halide schedules for live bilateral filtering (see Chapter 29.3, Halide at a high level: Decoupling Algorithms from Schedules).

On a single scanline the grid has only two axes, position and intensity, so each step of the pipeline can be followed with its numbers and its cost counted against the brute-force filter (Figure 10.1.11⛶ fullscreen).

The grid also makes the affinity lesson geometric: the third axis is the affinity, drawn as a coordinate. Two pixels far apart in value are literally far apart in the grid, so a local blur, which only mixes nearby cells, reduces their contribution according to distance and the blur support. The bilateral's "do these pixels belong together?" becomes, in the grid, simply "are these cells close?"

Because the grid is coarse and the blur is a plain convolution, the whole edge-aware filter runs fast enough to be interactive, and the authors put that speed to work on a range of applications.

A caveat: the curse of dimensionality. The grid buys speed by coarsening space and range, but its memory grows with the dimension of the cells. A grayscale bilateral needs a 3-D grid (two spatial axes plus one range axis); a color bilateral, or a cross-bilateral keyed on several guide channels, needs a range axis per channel, a 5-D grid for RGB, or more dimensions for additional guide channels whose cell count grows exponentially. And most of those cells are empty: real images occupy only a thin sheet in the space (exactly what the 3-D grid above makes visible: the occupied cells hug the image's intensity surface), so a dense array is enormously wasteful. Andrew Adams and colleagues answer this with the permutohedral lattice (Adams, Baek & Davis 2010 (permutohedral lattice)), a sparse, simplex-based tessellation of high-dimensional space that splats, blurs and slices without allocating a dense Cartesian grid exponential in the feature dimension; runtime still depends on dimension, samples and occupied vertices, making high-dimensional (color, cross-, even non-local) bilateral filtering practical. Read the other way, the bilateral grid is just the low-dimensional, dense special case of this general high-dimensional Gaussian filtering.

How much the grid saves depends on both widths; timing the two filters on the same image shows where it overtakes the brute force (Figure 10.1.16⛶ fullscreen). The brute force visits a window whose area grows as $\sigma_s^2$, so its cost climbs with the spatial width. The grid's cost does the opposite: a wider $\sigma_s$ means fewer, larger cells, and once the grid is small, what remains is the splat and the slice, which visit every pixel once whatever the widths. With both filters compiled from Halide and run on one thread, the grid overtakes the brute force at a spatial width of only one or two pixels, and at $\sigma_s=8$ it is about 67 times faster on the machine we measured, at a cost of a small error concentrated at strong edges, where the coarse cells blur most.

10.1.5 Bilateral upsampling and learning🔗⧉

Many operations are too expensive to run on a full-resolution photograph: a stereo matcher, an optimization, a large network. A common answer is to run them on a small copy and bring the result back to full resolution. Ordinary interpolation blurs the result across edges, so a depth map or a mask bleeds from the foreground onto the background. The full-resolution image is still available, though, and its edges say where the result should change. The bilateral filter's range weight is the tool for using them, and this section follows that idea from simple upsampling to layers inside neural networks.

Joint bilateral upsampling. Kopf and colleagues upsample a low-resolution solution $S$ with the full-resolution image as the guide $J$ Kopf et al. 2007 (joint bilateral upsampling). Each output pixel $p$ averages the low-resolution values near its position $p_\downarrow$ in the small image, weighted by distance in the small image and by similarity in the guide:

It is the cross bilateral filter of the previous section with the two inputs at different resolutions. A value from the other side of an edge in the guide gets almost no weight, so the result changes where the photograph does. The guide is not a measurement of the quantity being upsampled, though: a painted stripe is a color edge with no depth edge, and a foreground object the color of its background has a depth edge the guide cannot see. Chapter 29.1, Algorithmic speedups gives the procedure in pseudocode and these failure cases in detail.

Upsampling an operator rather than an image. Chen and colleagues upsample the effect of an operator instead of its output Chen et al. 2016 (bilateral guided upsampling). They run the operator on a small copy of the input and fit, in a coarse bilateral grid, a local affine color transform per cell: the one that best maps the small input to the small output. The grid is then sliced at full resolution, with each pixel's own position and intensity, and the sliced transform is applied to the full-resolution input. Because the transform is applied to the original pixels, the output keeps all of their detail, and the operator only has to get the large-scale changes right. This works for tone adjustments, color grading and many stylizations, anything a smoothly varying color transform can express, and fails for operations that move or create pixels, such as warps or added texture.

Learning the grid: HDRnet. HDRnet goes one step further and predicts the grid of affine coefficients with a network instead of fitting it to an operator's output (Gharbi et al. 2017, Deep Bilateral Learning). A convolutional network runs on a low-resolution copy of the image and predicts, per grid cell, the coefficients of a local affine color transform. At full resolution the grid is sliced at each pixel's position and at a guide value, itself a learned per-pixel function of the color, and the sliced transform is applied to the pixel. The expensive network stays at low resolution, where it can see the whole image, and the edge-aware slice carries its decision back to every pixel. The bilateral grid becomes a learned, edge-preserving upsampling layer, fast enough to run in a phone's viewfinder (Figure 7.3.4⛶ fullscreen, Figure 7.3.5⛶ fullscreen). Chapter 7.3, Deep learning shows more of its results.

The bilateral solver. Upsampling averages the low-resolution values; some problems instead need to clean them. A network's depth map or segmentation is often noisy, has holes, and comes with a confidence that varies from pixel to pixel. Barron and Poole's fast bilateral solver finds the image $x$ that stays close to a target $t$ where the confidence $c$ is high and is smooth under the bilateral affinity $\hat W$ of the guide everywhere Barron & Poole 2016 (fast bilateral solver):

Solved pixel by pixel this would be a huge system, but the affinity is a bilateral filter, so the problem can be splatted into a bilateral grid, solved there with a small linear system, and sliced back. The solution is differentiable with respect to the target and the confidence, so the solver can sit inside a network and be trained through. It has been used to refine stereo depth, to upsample depth and semantic segmentation, and for colorization.

Trainable edge-aware layers. Two lines of work made the bilateral filter itself a network layer. Jampani, Kiefel and Gehler learned the filter's weights on the permutohedral lattice, the high-dimensional generalization of the grid Adams, Baek & Davis 2010 (permutohedral lattice), instead of fixing them to a Gaussian Jampani, Kiefel & Gehler 2016. Wu and colleagues made the guided filter of Chapter 10.6, Guided image filtering differentiable, with a guide image learned by a small network, so that a network can run at low resolution and be upsampled by a guided filter trained end to end with it Wu et al. 2018 (deep guided filter).

Sharpening segmentations. Semantic segmentation networks produce coarse labels whose boundaries do not follow the objects. A fully connected conditional random field, solved with bilateral filtering in its inner loop, snaps the labels to the image's edges: pixels with similar color and position are encouraged to share a label Krähenbühl & Koltun 2011 (dense CRF). Zheng and colleagues unrolled that inference as layers of a recurrent network so that it could be trained together with the network that produces the labels Zheng et al. 2015 (CRF-as-RNN).

Bilateral grids in neural rendering. Reconstructions of a scene from many photographs, such as neural radiance fields and Gaussian splatting, assume that the same point has the same color in every view. A phone's processing breaks that assumption: exposure, white balance and local tone mapping change from frame to frame, and the reconstruction absorbs the differences as floating artifacts. Wang and colleagues give each view a small learned bilateral grid of affine color transforms, which models that processing and is discarded after training, so the reconstructed scene stays consistent Wang et al. 2024 (BilaRF).

The pattern is the same throughout. A computation too expensive or too coarse for the full image runs at low resolution or on low-quality data, and an edge-aware structure built from the full-resolution image carries the result back, keeping its edges where the photograph's are. What the guide cannot supply is information it does not contain.

10.1.6 Anisotropic diffusion and robust estimators🔗⧉

The bilateral filter is one member of a family, and the other two members were derived from different starting points and arrive at nearly the same arithmetic. Seeing them together explains what the range weight is really doing.

Read the filter as an estimator first. Averaging a neighborhood is what you do when you believe every neighbor is a measurement of the same value corrupted by Gaussian noise, which is the quadratic penalty $\rho(x)=x^2$. A quadratic grows without bound, so one very different neighbor, the pixel across an edge, pulls the estimate as far as it likes. A robust error norm saturates instead: past some difference a neighbor stops mattering more the further away it is. Edge preservation and outlier rejection are then the same operation, and the bilateral range weight is a particular choice of robust norm (Figure 10.1.18⛶ fullscreen).

A weighted-difference update resembles one discrete diffusion step. For a symmetric pairwise energy, its force is proportional to the penalty derivative $\psi(t)=\rho'(t)$, while the weight is proportional to $\psi(t)/t$. Weight and influence are therefore different curves: a Gaussian weight yields $\psi(t)=t\exp(-t^2/(2\sigma_r^2))$, up to scale. Perona–Malik diffusion uses a local gradient-dependent conductance Perona & Malik 1990; a finite-neighborhood bilateral update is related but not the same discretization. Some penalties saturate, while others, such as the Lorentzian, grow without bound despite a diminishing influence Black et al. 1998.

Perona and Malik started from diffusion rather than from estimation Perona & Malik 1990. Blurring with a Gaussian is the same as running the heat equation $\partial_\tau I = \operatorname{div}(\nabla I)$ for a time $\tau$ that grows with the blur's variance: intensity flows from brighter to darker neighbors until they agree, across edges as readily as within regions. Their change was to let the flow depend on the local gradient,

with a conductance $g$ near one for small gradients and near zero for large ones, for example $g(s)=1/(1+(s/\kappa)^2)$ or $g(s)=e^{-(s/\kappa)^2}$ with an edge scale $\kappa$. Heat then flows freely inside a region and barely crosses a strong edge, so iterating the equation smooths regions while boundaries survive. Two refinements matter in practice. With these conductances the flux $g(s)\,s$ decreases once $s$ exceeds a value near $\kappa$, so diffusion can steepen an edge rather than merely preserve it, and the continuous equation is ill-posed in that regime; practical schemes measure the gradient on a slightly smoothed image and take explicit time steps small enough to stay stable. A bilateral pass is not a discretization of this equation: it averages over a whole neighborhood in one large step, where diffusion takes many small nearest-neighbor steps.

A tensor can express directional diffusion. Diffusion-tensor imaging estimates physical water diffusion from measurements; it is not an image-smoothing algorithm that recovers anatomy by running Perona–Malik.

The practical consequence is that these are three vocabularies for one design decision: how much should a dissimilar neighbor count? The bilateral answers with a range Gaussian, Perona and Malik with a conductance that falls with the gradient Perona & Malik 1990, robust statistics with an influence function. Chapter 13.3, Gaussian, median, and bilateral filtering compares their behavior on noise, and Chapter 10.4, Local Laplacian filters takes the same question to a multiscale decomposition.