Image Pyramid

An image pyramid is a stack of copies of one image at decreasing resolutions, tapering from the full-resolution image at the bottom to a few pixels at the top. Pyramids let an algorithm look at an image at several scales for little more memory than the image itself. They are used for estimating large motions, detecting objects and features of any size, blending and compressing images, and building the multi-scale features of modern networks.

Definition

An image pyramid is a multi-scale representation of an image, formed by repeatedly applying a low-pass filter and subsampling. In the usual dyadic pyramid, each level has half the width and height of the level below. A Gaussian pyramid stores the smoothed images themselves; a Laplacian pyramid stores the differences between consecutive levels, which isolate the detail present at one scale and not at the next.

Intuition

Seen from far away, a photograph loses its fine detail, and a Gaussian pyramid imitates this. Smoothing removes detail too fine to survive at the next resolution, and subsampling discards pixels that have become redundant once that detail is gone. Large structures shrink to a few pixels at the coarse levels, where a fixed-size operator can find them.

The Laplacian pyramid records what each step threw away. Subtracting a blurred, re-enlarged copy of a level from the level itself leaves only the detail at that scale: fine texture in the bottom band, larger edges and blobs in the bands above. Adding the bands back together, from the top down, restores the original exactly.

Formal Definition

Gaussian pyramid

Let $g_0$ be the image. Burt and Adelson (1983) define each further level with an operation they call REDUCE: convolve with a small kernel $w$ and subsample by 2,

$$g_{l+1}(i, j) = \sum_{m=-2}^{2} \sum_{n=-2}^{2} w(m, n)\, g_l(2i + m,\; 2j + n), \qquad l = 0, \dots, N - 1.$$

Their kernel, the “generating kernel”, is separable, $w(m, n) = \hat{w}(m)\,\hat{w}(n)$, with five taps

$$\hat{w} = \left[\tfrac{1}{4} - \tfrac{a}{2},\; \tfrac{1}{4},\; a,\; \tfrac{1}{4},\; \tfrac{1}{4} - \tfrac{a}{2}\right].$$

The constraints that fix this form are normalization, symmetry, and equal contribution: every pixel of a level contributes the same total weight to the level above. With $a = 0.4$, repeated REDUCE steps are equivalent to convolving $g_0$ with functions that closely resemble Gaussians, which gives the pyramid its name. A sampled Gaussian or the binomial kernel $\tfrac{1}{16}[1, 4, 6, 4, 1]$ (the case $a = 0.375$) is common in practice.

Laplacian pyramid

The reverse operation, EXPAND, upsamples a level to the size of the one below by interpolating between its samples, in practice by inserting zeros and convolving with the same kernel scaled by 4. The Laplacian pyramid is

$$L_l = g_l - \operatorname{EXPAND}(g_{l+1}), \quad l = 0, \dots, N - 1, \qquad L_N = g_N.$$

Each $L_l$ is a band-pass image, equivalent to convolving the original with the difference of two Gaussian-like functions of neighboring scales. Burt and Adelson noted that this resembles applying a Laplacian operator whose scale doubles from level to level, hence the name. The original is recovered by reversing the construction:

$$g_N = L_N, \qquad g_l = L_l + \operatorname{EXPAND}(g_{l+1}), \quad l = N - 1, \dots, 0.$$

Reconstruction is exact for any kernel, because each band stores exactly what EXPAND fails to predict.

Properties

• Memory. Level $l$ has about $4^{-l}$ times as many samples as the image, so the whole pyramid has $\sum_{l=0}^{N} 4^{-l} < 1/(1 - \tfrac{1}{4}) = \tfrac{4}{3}$ times as many: about a third more than the image. The Laplacian pyramid is overcomplete by the same factor of $4/3$, whereas an orthonormal wavelet transform has exactly as many coefficients as pixels.
• Octave spacing. Each level halves the band limit and the sampling rate together; consecutive levels are one octave apart in scale.
• Localization. Laplacian coefficients are localized in both space and spatial frequency, so a coefficient describes the detail at a given place and scale.
• Shift sensitivity. Subsampling makes the coarse levels depend on how the image is aligned with the sampling grid: shifting the image by one pixel does not shift the coarse levels by half a pixel exactly.

Relation to scale space and wavelets

A Gaussian scale space is the continuous family of images obtained by convolving the original with Gaussians of every standard deviation $\sigma$, all kept at full resolution. A Gaussian pyramid approximately samples this family at octave steps of $\sigma$ and, because each level is band-limited, subsamples it, so it trades the continuous family’s redundancy for compactness. Methods that need finer scale steps sample several scales per octave and subsample only once per octave.

Wavelet transforms also decompose an image into band-pass components at octave-spaced scales, but critically sampled wavelets have no redundancy. The price is aliasing within the subbands. Simoncelli and Freeman’s comparison lists the Laplacian pyramid as overcomplete by $4/3$ and dyadic wavelets by 1.

Steerable pyramids (Simoncelli and Freeman, 1995) split each band further by orientation. Their basis functions are directional derivative operators, from which image gradients and higher derivatives can be read at every scale; the transform is self-inverting (its inverse is its transpose), and the subbands are essentially free of aliasing. With $k$ orientation bands the representation is overcomplete by $4k/3$.

Examples

The code builds Gaussian and Laplacian pyramids with Burt and Adelson’s kernel, checks the memory cost and the reconstruction, and shows what subsampling does without smoothing.

import numpy as np
from scipy import ndimage

# Burt and Adelson's 5-tap generating kernel with a = 0.4: [0.05, 0.25, 0.4, 0.25, 0.05].
a = 0.4
KERNEL = np.array([0.25 - a / 2, 0.25, a, 0.25, 0.25 - a / 2])


def blur(image):
    """Separable smoothing with the generating kernel; mirror the borders."""
    out = ndimage.convolve1d(image, KERNEL, axis=0, mode="reflect")
    return ndimage.convolve1d(out, KERNEL, axis=1, mode="reflect")


def reduce(image):
    """Smooth, then keep every second row and column."""
    return blur(image)[::2, ::2]


def expand(image, shape):
    """Upsample to `shape` by inserting zeros, then interpolate with the same kernel."""
    up = np.zeros(shape)
    up[::2, ::2] = image
    return 4.0 * blur(up)  # the factor 4 restores the mean lost to the inserted zeros


def gaussian_pyramid(image, levels):
    pyr = [image]
    for _ in range(levels - 1):
        pyr.append(reduce(pyr[-1]))
    return pyr


def laplacian_pyramid(image, levels):
    g = gaussian_pyramid(image, levels)
    lap = [g[k] - expand(g[k + 1], g[k].shape) for k in range(levels - 1)]
    return lap + [g[-1]]  # the top level is the coarsest Gaussian image


def reconstruct(lap):
    image = lap[-1]
    for band in reversed(lap[:-1]):
        image = band + expand(image, band.shape)
    return image


rng = np.random.default_rng(0)
image = ndimage.gaussian_filter(rng.standard_normal((257, 343)), 1.5)  # odd sizes on purpose

g = gaussian_pyramid(image, 5)
lap = laplacian_pyramid(image, 5)
print("level shapes:", [level.shape for level in g])
samples = sum(level.size for level in lap)
print(f"samples in pyramid / samples in image: {samples / image.size:.4f}")
print(f"max reconstruction error: {np.abs(reconstruct(lap) - image).max():.2e}")

# Aliasing: vertical stripes with a period of 2.5 pixels, subsampled by 2.
x = np.arange(250)
stripes = np.tile(np.cos(2 * np.pi * x / 2.5), (250, 1))
naive = stripes[::2, ::2]
smoothed = reduce(stripes)
k = np.abs(np.fft.rfft(naive[0]))[1:].argmax() + 1
print(f"naive subsampling: amplitude {np.abs(naive).max():.3f}, "
      f"apparent period {naive.shape[1] / k:.1f} samples = {2 * naive.shape[1] / k:.1f} original pixels")
print(f"smoothing first:   amplitude {np.abs(smoothed[4:-4, 4:-4]).max():.3f}")

Output:

level shapes: [(257, 343), (129, 172), (65, 86), (33, 43), (17, 22)]
samples in pyramid / samples in image: 1.3355
max reconstruction error: 1.11e-16
naive subsampling: amplitude 1.000, apparent period 5.0 samples = 10.0 original pixels
smoothing first:   amplitude 0.026

The five levels hold 1.34 times as many samples as the image, close to the limit of $4/3$, and the image is reconstructed to rounding error even though its sides are odd. The stripes, with a period of 2.5 pixels, are too fine to be represented after halving the sampling rate. Subsampled directly, they do not disappear: they turn into false stripes with a period of 10 pixels at full contrast. Smoothing first reduces them to 2.6% of their amplitude, which is the response of the generating kernel at that frequency.

Common Misconceptions

• “Downsampling means taking every second pixel.” Without smoothing first, detail finer than the new sampling rate folds into false, coarser patterns (aliasing), as the example shows.
• “The Laplacian pyramid compresses the image by itself.” It has a third more samples than the image. Burt and Adelson obtained compression by quantizing the band images, whose values are mostly near zero, and coding them with few bits.
• “A pyramid is the same as scale space.” A pyramid is a subsampled, octave-spaced sampling of a scale space, not the continuous family.

Pitfalls

• Border handling. The kernel reaches past the image edge, and the choice of extension (reflection, replication, zero) changes the values near the border at every level, increasingly so at coarse levels where a few pixels span much of the image. In their spline paper, Burt and Adelson extrapolated by reflection and inversion about the edge pixel. Whatever the choice, use the same one in REDUCE and EXPAND.
• Odd sizes. Halving an odd dimension does not give an integer. Burt and Adelson used images of size $M 2^N + 1$ so that every level aligns exactly; most implementations instead round up and pass the target shape to EXPAND, as the example does. Mismatched shapes are a common source of off-by-one bugs.
• Scaling coordinates. A position or displacement measured at level $l$ must be multiplied by $2^{l}$ to express it at full resolution.

Where It Is Used

• Coarse-to-fine estimation. Halving the resolution halves every displacement, so methods valid only for small motions, such as Lucas–Kanade, are run from the top of a pyramid down, each level refining the estimate passed from the level above. This is the standard way to estimate large motions in classical optical flow estimation and feature tracking; see coarse-to-fine estimation, which also covers its failure on small, fast objects.
• Multi-scale detection. A detector with a fixed window finds objects of one size; scanning the window over every level of a pyramid finds them at all sizes. Lin et al. note that such “featurized image pyramids” were heavily used in the era of hand-engineered features, and that detectors like the deformable part model needed dense scale sampling, such as 10 scales per octave.
• Scale-invariant features. SIFT (Lowe, 2004) finds keypoints as extrema of difference-of-Gaussian images in a scale space organized by octaves: within an octave the image is blurred at several scales, and the next octave starts from a blurred image subsampled by taking every second pixel.
• Image blending. Burt and Adelson’s multiresolution spline (1983) joins two images by blending each Laplacian band separately, with a transition zone proportional to the wavelengths in that band, and then reconstructing. Coarse differences, such as overall brightness, blend over a wide zone, while fine detail meets along a narrow seam. For a region of arbitrary shape, a binary mask is turned into a Gaussian pyramid whose levels serve as per-band blending weights. The same pyramid blending is one of the compositing methods used in panorama stitching.
• Compression and progressive transmission. Quantizing the Laplacian bands compresses the image, and sending the top levels first lets a receiver see a coarse version that sharpens as finer bands arrive (Burt and Adelson, 1983).
• Antialiased resampling. Shrinking an image by a large factor, as in image warping, avoids aliasing by sampling a prefiltered pyramid level; texture mapping does this with mipmaps.
• Feature pyramids in CNNs. A convolutional network’s successive stages already produce feature maps at decreasing resolutions. Feature Pyramid Networks (Lin et al., 2017) add a top-down pathway with lateral connections, so that every scale carries strong semantic features, at marginal extra cost and without running the network on every level of an image pyramid. Optical flow networks such as PWC-Net likewise build learned feature pyramids instead of pyramids of pixels.