13.4 Pyramid and wavelet denoising🔗⧉

A sharp boundary between a wall and a window produces large coefficients in a few wavelet bands. In a smooth part of the wall, most detail coefficients are small. Noise adds fluctuations throughout the transform. We can therefore reduce noise by shrinking detail coefficients before reconstructing the image.

The transforms are introduced in Chapter 6.6, Pyramids and wavelets. After the first section, the derivations use an orthonormal wavelet transform and additive white Gaussian noise for the basic derivation. This assumption matters: an orthonormal transform preserves the covariance $\sigma^2 I$, and Gaussianity makes the transformed noise coefficients independent. White noise need not be Gaussian in general.

13.4.1 The method: shrink the bands of a pyramid🔗⧉

Recall how a Laplacian pyramid is built (Chapter 6.6, Pyramids and wavelets; Burt & Adelson 1983 (Laplacian pyramid)). Blur the image and halve its resolution, repeatedly, to get the Gaussian pyramid: the same photograph at 1, 1/2, 1/4 of its size, each level smoother than the last. Each Laplacian level is the difference between one Gaussian level and the next coarser one upsampled back to its size, so it holds the detail that the blur removed, roughly one octave of spatial frequency. The last level is the small, blurred residual. Collapsing reverses the construction exactly: upsample the residual, add the next level, and repeat. Whatever values the detail levels hold, the collapse adds them back in.

That exactness is what makes the pyramid a denoising tool. In a smooth region of the photograph, a detail level is close to zero, and what it does hold is mostly noise. At an edge or in texture, it holds large values. Noise, spread over all frequencies, adds small fluctuations to every level. So the method is three steps:

levels, residual = laplacian_pyramid(y, number_of_levels)
for j, band in enumerate(levels):
    t = k * noise_std_of_band(j)          # threshold follows that band's noise level
    levels[j] = shrink(band, t)           # small values toward zero, large ones kept
x_est = collapse(levels, residual)        # residual left as it is

The shrinkage curve, often called coring in this setting, is flat at zero for small inputs, the dead zone where noise lives, and passes large inputs through (Figure 13.4.2⛶ fullscreen). Edges survive because their coefficients are large relative to the noise; flat regions lose their grain because theirs are not (Figure 13.4.1⛶ fullscreen). The rest of this chapter makes each step precise: which curve to use, how to set the threshold of each band, and why a better statistical model of the coefficients shrinks less at real structure.

A Laplacian pyramid is overcomplete (it stores about 4/3 as many values as the image) and its bands are not orthogonal, so noise in its levels is correlated and its strength differs from band to band. An orthonormal wavelet transform removes both complications, which is why the derivations below use it. The Laplacian pyramid remains the easier transform to build and to reason about, and its lack of subsampling artifacts is an advantage we return to with overcomplete transforms.

Working band by band also matches how we see. The visual system analyzes an image through channels tuned to spatial frequency and orientation, and its sensitivity to contrast depends on frequency: it peaks at intermediate frequencies and falls for very fine detail (Chapter 3.4, Spatial and temporal vision; Campbell & Robson 1968). The resemblance between band-pass channels and pyramid levels was part of the original motivation for the Laplacian pyramid. For a denoiser it means that errors in different bands are not equally visible. At ordinary viewing distances, the finest level of a high-resolution photograph often falls beyond the peak of that sensitivity, so residual fine grain is less objectionable than the blotchy, low-frequency noise of a coarser band, which a fine-band threshold leaves untouched. Color detail is seen more coarsely still (Spatial and temporal vision#Chrominance spatial vision), which is why chroma bands can be cored harder than luminance bands (Gaussian, median, and bilateral filtering#Denoise color more than brightness). Where the thresholds should be set is therefore partly a statistical question, how much noise each band carries, and partly a perceptual one, how much of what remains a viewer will see.

13.4.2 Sparse in the right basis🔗⧉

Many photographs have wavelet-coefficient histograms concentrated near zero with broader tails than a Gaussian. Large coefficients occur at boundaries and structured detail; numerous small coefficients describe smooth regions. This is statistical sparsity, not a claim that every small coefficient is disposable or every photograph has the same distribution (Figure 13.4.3⛶ fullscreen).

The argument is a claim about a distribution, so look at the distribution. Take a photograph, take one subband of its wavelet transform, and plot the histogram of its coefficients on a logarithmic vertical axis. For most photographs it is a sharp spike at zero with long tails, close to a Laplacian or, better, a hyper-Laplacian with an exponent below one. Do the same for pure Gaussian noise of the level you expect: on the same axis it is a downward parabola, with no spike and no tails. A threshold is a cut through both populations. The fraction of the noise histogram inside the cut is noise removed; the fraction of the image histogram inside it is signal removed with it. Moving the threshold trades one against the other, and the two curves show the trade directly.

Let $c_y=c+n_c$ be a measured coefficient. Under an orthonormal transform of white Gaussian noise, $n_c$ has the original noise variance $\sigma^2$. Signal and noise distributions overlap, so a coefficient cannot be classified perfectly from its magnitude. Shrinkage is an estimate under uncertainty. It will remove some real detail and retain some noise.

For a redundant pyramid, the noise variance depends on each analysis filter's energy, and coefficients can be correlated. If band $j$ uses filter $h_j$, white input noise contributes variance $\sigma_j^2=\sigma^2\sum_q h_j(q)^2$. Use these band-dependent levels or measure them by passing simulated noise through the actual transform, with the same normalization and boundaries.

13.4.3 Thresholding, hard and soft🔗⧉

Hard thresholding retains a coefficient above a magnitude threshold and removes it otherwise:

Soft thresholding also reduces the magnitude of each survivor:

A hard threshold of 3 maps coefficients 2.9 and 3.1 to 0 and 3.1; a soft threshold maps them to 0 and 0.1. Hard thresholding avoids shrinking large coefficients but introduces a discontinuity. Soft thresholding is continuous but reduces contrast. Neither promises untouched edges (Figure 13.4.4⛶ fullscreen).

For $M$ independent Gaussian noise coefficients, the universal threshold $t=\sigma\sqrt{2\log M}$ is a conservative choice intended to reject large pure-noise excursions (Donoho & Johnstone 1994). Its assumptions are clearest in an orthonormal basis. Applying it independently to bands uses each band's coefficient count and noise level; this is a specified implementation choice rather than a guarantee for an arbitrary pyramid.

coeffs = orthonormal_wavelet_transform(y)
v = coeffs.finest_diagonal
sigma_est = median(abs(v - median(v))) / 0.6745
for j in detail_band_indices(coeffs):
    band = coeffs[j]
    threshold = sigma_est * sqrt(2 * log(number_of_coefficients(band)))
    coeffs[j] = soft_threshold(band, threshold)
x_est = inverse_wavelet_transform(coeffs)

The median absolute deviation estimate assumes the chosen high-frequency band is sufficiently dominated by Gaussian noise. Fine texture can inflate it; use a calibrated level when available. We leave the coarse approximation band unchanged in this simple version to avoid suppressing image brightness and broad shading. It still contains noise and can be estimated separately in more elaborate methods.

13.4.4 Shrinkage as a MAP estimate🔗⧉

A Laplace prior $p(c)\propto\exp(-|c|/s)$ and Gaussian coefficient noise give the maximum a posteriori estimate

For a positive solution, differentiating gives $c=c_y-\sigma^2/s$; the negative case gives $c=c_y+\sigma^2/s$. If these would cross zero, the minimum is at zero. Together the cases give soft thresholding with $t=\sigma^2/s$. Higher noise increases shrinkage, while a broader prior permits larger coefficients.

Equivalently, soft thresholding minimizes $\tfrac12(c-c_y)^2+t|c|$, making it the proximal operator of the absolute-value penalty. Chapter 13.7, Learned denoisers uses this operation in an iterative sparse estimator and then learns its parameters. This coefficientwise closed form relies on the separable model; a redundant analysis penalty does not generally reduce to one transform-threshold-inverse pass.

Other priors produce different shrinkage curves. A hyper-Laplacian penalty $|c|^\alpha$ with $0<\alpha<1$ shrinks large coefficients less aggressively, but makes the scalar problem nonconvex. Its solution requires care about thresholds and competing minima; some particular exponents admit specialized closed forms.

This scalar calculation connects several chapters. A data term plus a prior gives the inverse-problem formulation in Chapter 7.1, Imaging as Bayesian inference with probabilistic priors. Soft thresholding is the simple proximal step that more elaborate iterative solvers can replace by a denoiser in Chapter 13.11, Denoising as a universal prior. Unrolling such an iteration produces the trainable stages discussed in Chapter 13.7, Learned denoisers. The progression starts with an estimator whose threshold we can derive, then changes either the prior or the way the estimator is computed.

13.4.5 Ringing, and why the transform is usually overcomplete🔗⧉

A critically sampled wavelet transform depends on the alignment of the image with its sampling grid. Shift an edge by one pixel and its coefficients change; thresholding them can then produce a different edge profile or ringing. Cycle spinning Coifman & Donoho 1995 (cycle spinning) reduces this dependence:

result = zero_image()
for shift in chosen_grid_shifts:
    moved = shift_image(y, shift)
    estimate = wavelet_shrinkage(moved)
    result += undo_shift(estimate, shift)
result /= number_of_shifts

For circular shifts, unshift the entire reconstruction before aggregation; other boundary rules require valid-region handling. Use at least one shift. An undecimated transform provides another way to avoid subsampling dependence, at a cost in storage and computation. Its synthesis and threshold conventions matter, so do not assume every redundant transform is identical to a particular cycle-spinning implementation (Figure 13.4.5⛶ fullscreen).

All shifted estimates use the same noisy input. Their errors are correlated, so averaging them does not earn the $1/\sqrt N$ guarantee of independent captures. The benefit comes from reducing dependence on a particular transform grid and combining estimates with different artifacts.

For cycle spinning, shift the image, transform and denoise it, then undo the shift before combining it with the other results. Different shifts place the same edge differently relative to the wavelet grid. Their artifacts therefore need not coincide, even though the reconstructed scene should align after the shifts are undone. Averaging can reduce that grid dependence at a cost in computation. Averaging without undoing the shifts would instead blur the image, and a finite set of shifts only samples the possible grid alignments.

13.4.6 Coring band by band, in practice🔗⧉

Thresholds can vary by scale and orientation, and that is where the operator makes different decisions about different structures. Fine-scale coefficients carry sensor noise and small texture; coarser coefficients describe broader transitions. Raising only the fine-scale threshold can reduce grain without flattening broad contrast, while raising every threshold produces the smooth, plastic rendering that gives aggressive denoising its reputation. Spatially varying thresholds add a further control where noise or motion confidence changes across the frame. The band noise levels themselves depend on the analysis filters, their normalization and any correlation already in the input, so the finest band does not universally hold almost all of the noise.

13.4.7 Coring and the Gaussian scale mixture🔗⧉

Bayesian coring estimates a coefficient using a fitted distribution rather than selecting a threshold only by convention (Simoncelli & Adelson 1996). A conditional-mean estimator generally produces a smooth shrinkage curve. It differs from the MAP estimator, which selects a posterior mode.

Coefficients also provide context for one another: an edge often produces large magnitudes at nearby positions and scales. The dependence is in the magnitudes, not the values. Neighboring coefficients of a photograph are close to uncorrelated, since their signs flip as often as not, yet their magnitudes are strongly correlated: where one is large, its neighbors tend to be large too Simoncelli & Olshausen 2001. A model that treats coefficients as independent misses exactly that. A Gaussian scale mixture models a neighborhood vector as $c=\sqrt z\,u$, where $u$ is Gaussian and $z$ is an independent positive multiplier. The multiplier represents local activity. BLS-GSM, Bayesian least squares with a Gaussian scale mixture, averages the conditional estimates over uncertainty in that multiplier (Portilla et al. 2003 (BLS-GSM)). A coefficient supported by neighboring activity can then be shrunk less than an isolated coefficient of the same magnitude.

The transform matters as well. A separable wavelet has only horizontal, vertical and a mixed diagonal band, so an edge at an intermediate angle spreads its energy over several bands and several coefficients, each of which is small enough to be shrunk. Oriented transforms such as the steerable pyramid Simoncelli & Freeman 1995 or curvelets Candès & Donoho 2004 (curvelets) have more orientations, so a straight edge concentrates into fewer, larger coefficients that survive shrinkage. They are overcomplete, which costs memory and time, and they help most on images dominated by long edges.

This extends scalar shrinkage to a local statistical model in an oriented pyramid. It still depends on the transform, the assumed noise covariance, and the fitted image statistics; complex or unfamiliar texture can violate those assumptions.

Consider two wavelet coefficients of the same magnitude. One lies among many large coefficients along an edge; the other is isolated in an otherwise quiet neighborhood. A fixed threshold treats them alike. A model that estimates local signal variance can interpret the first as part of real structure and shrink it less, while shrinking the isolated coefficient more. A Gaussian scale mixture expresses this by allowing a shared local scale to modulate a group of coefficients. Estimating that scale uses context across nearby positions, orientations, or levels. The benefit comes from the joint model, not from assuming that every small coefficient is noise or every large coefficient is signal Portilla et al. 2003 (BLS-GSM).