Image Pyramids

A Gaussian pyramid is a sequence of images $G_0, G_1, \dots, G_L$ where $G_0$ is the original image and each subsequent level is obtained by low-pass filtering $G_{\ell}$ with a Gaussian-like kernel and then downsampling by a factor of 2 in each dimension. Formally,

where $h$ is a small separable smoothing kernel (classically the 5-tap "generating kernel" of Burt and Adelson) and $*$ denotes 2D convolution. Each level is a scale-space sample of the image: $G_\ell$ approximates the original scene as it would appear blurred by a Gaussian of standard deviation roughly $\sigma \cdot 2^{\ell}$ and viewed at $1/2^{\ell}$ the linear resolution. This connects directly to the continuous notion of scale-space, in which an image is embedded in a one-parameter family $I(x,y;\sigma) = I_0 * G_\sigma$; the Gaussian pyramid is a discretized, computationally efficient version of that family, sampled at octave (doubling) intervals rather than continuously in $\sigma$. A Laplacian pyramid stores what smoothing removes at each step, rather than the smoothed image itself. Define

for $\ell = 0, \dots, L-1$, and let $L_L = G_L$ be the coarsest residual (the "DC" band). Each $L_\ell$ is approximately a band-pass filtered version of the original image — the frequencies present in $G_\ell$ but removed by the blur-and-downsample step that produced $G_{\ell+1}$. Because the name "Laplacian" pyramid comes from the fact that $G_\ell - G_\ell * h$ approximates the response of a Laplacian-of-Gaussian operator, this construction is a fast, approximate multi-scale edge/texture decomposition — a discretized relative of a wavelet transform, but with two crucial practical differences: it is over-complete by a small, controlled factor (each level after upsampling is compared at the same resolution as the level above it), and its filters are simple, separable, and identical at every level, which is what makes it fast. The two pyramids are dual views of one representation: the Gaussian pyramid is redundant but easy to search over scale (each level is a plausible standalone image, useful for detecting objects or matching features at that scale); the Laplacian pyramid is compact and additive (summing the bands reconstructs the original), useful for compression and for editing operations — such as blending — that should be done independently at each spatial-frequency band. Both rely on the smoothing and resampling machinery introduced under Gaussian and Median Filtering and connect forward to general multi-resolution Image Resampling and Warping.

Building the Gaussian pyramid. Start from $G_0 = I$, the input image. At each level, convolve with a small low-pass kernel to prevent aliasing, then keep every other row and column. The classic Burt–Adelson kernel is a separable 5-tap filter with weights approximating a Gaussian, e.g. $\tfrac1{16}[1,4,6,4,1]$, applied along rows then columns (an equivalent 2D kernel of $\sigma \approx 1$). The filter must remove frequencies above the new Nyquist limit; halving the resolution halves the Nyquist frequency, so blurring first is not optional — it is what the sampling theorem demands.

Building the Laplacian pyramid. For each level $\ell$, upsample $G_{\ell+1}$ back to the resolution of $G_\ell$ (insert zeros between samples, then convolve with $4h$ so that interpolation, not mere zero-stuffing, fills in the new samples — the factor of 4 compensates for the energy lost by zero insertion in 2D), and subtract from $G_\ell$. The residual $L_\ell$ captures detail lost between consecutive Gaussian levels.

Reconstruction. Because $L_\ell = G_\ell - \text{UPSAMPLE}(G_{\ell+1})$, we can recover $G_\ell$ exactly from $G_{\ell+1}$ and $L_\ell$ by rearranging: $G_\ell = L_\ell + \text{UPSAMPLE}(G_{\ell+1})$. Starting from the coarsest level $G_L = L_L$ and working upward reconstructs every level, including the original image $G_0$. ```

# ---- Gaussian pyramid construction ---- function BUILD_GAUSSIAN_PYRAMID(I, num_levels): G[0] = I for l = 0 to num_levels - 1: blurred = CONVOLVE(G[l], KERNEL_5TAP) # low-pass, separable G[l+1] = DOWNSAMPLE_BY_2(blurred) # keep even rows & cols return G # G[0..num_levels]

# ---- Laplacian pyramid construction ---- function BUILD_LAPLACIAN_PYRAMID(G, num_levels): for l = 0 to num_levels - 1: expanded = UPSAMPLE_BY_2(G[l+1]) # zero-insert expanded = CONVOLVE(expanded, 4 * KERNEL_5TAP) # interpolate expanded = CROP_TO_SIZE(expanded, G[l]) # handle odd dimensions L[l] = G[l] - expanded L[num_levels] = G[num_levels] # coarsest band = DC residual return L # L[0..num_levels]

# ---- Exact reconstruction ---- function RECONSTRUCT(L, num_levels): R = L[num_levels] for l = num_levels - 1 down to 0: expanded = UPSAMPLE_BY_2(R) expanded = CONVOLVE(expanded, 4 * KERNEL_5TAP) expanded = CROP_TO_SIZE(expanded, L[l]) R = L[l] + expanded return R # equals original image G[0]

# ---- Coarse-to-fine search skeleton (e.g. optical flow / template match) ---- function COARSE_TO_FINE_SEARCH(pyramidA, pyramidB, num_levels): estimate = ZERO_MOTION for l = num_levels down to 0: estimate = UPSAMPLE_ESTIMATE(estimate, scale=2) # refine previous level's guess estimate = REFINE(pyramidA[l], pyramidB[l], estimate) return estimate ``` Two details matter for correctness. First, when the image dimensions are odd, downsampling and upsampling are not perfect inverses of each other in size; `CROP_TO_SIZE` (or equivalently padding rules fixed in advance) must be applied consistently during both construction and reconstruction, or the residual will not exactly cancel. Second, the boundary handling of the convolution (reflect, replicate, or zero-pad) must be identical in the blur used for downsampling and the blur used for upsampling-interpolation, otherwise reconstruction accumulates edge artifacts.

Time. Let the base image have $n$ pixels. Level $\ell$ has approximately $n / 4^{\ell}$ pixels (a factor of 2 fewer rows and 2 fewer columns). Convolving with the fixed 5-tap separable kernel costs $O(1)$ work per pixel (5 multiply-adds per pass, two passes for separability — a constant independent of image size), so building level $\ell$ costs $O(n/4^{\ell})$. Summing over all levels,

a convergent geometric series bounded by $\tfrac43 n$ regardless of how many levels $L$ are built. Building the entire pyramid — Gaussian or Laplacian — is therefore $O(n)$ time, the same order as a single full-resolution blur. This is the key efficiency argument: extending an algorithm to "search over all scales via a pyramid" multiplies the cost of a single-scale pass by at most a constant factor of about $4/3$, not by the number of levels or by $\log n$.

Space. The same geometric series bounds total storage: the whole pyramid occupies at most $\tfrac43 n$ pixels of extra memory (Gaussian pyramid) or exactly $\tfrac43 n$ for the Laplacian pyramid stored in place of the Gaussian one (the Laplacian levels are the same sizes as the corresponding Gaussian levels, and the pyramid needs no additional top-level storage beyond the base image). Both pyramids are $O(n)$-space structures — mild overcompleteness, not an explosion.

Number of levels. Levels can be built until the image reaches a minimal size (e.g. $1\times1$), giving $L = O(\log n)$ levels total, but note this does not change the total-time bound above — the geometric shrinkage guarantees the sum converges however many levels exist.

Correctness of reconstruction (Laplacian pyramid is a perfect decomposition). By construction,

Define the reconstruction operator applied bottom-up starting from $R_L := G_L$ and $R_\ell := L_\ell + \text{EXPAND}(R_{\ell+1})$. We show by downward induction on $\ell$ that $R_\ell = G_\ell$ for all $\ell$. Base case ($\ell = L$): $R_L = G_L$ by definition. Inductive step: assume $R_{\ell+1} = G_{\ell+1}$. Then

using the inductive hypothesis in the middle step. This holds exactly, in exact arithmetic, using any fixed, deterministic EXPAND operator (linear or not), because the construction and reconstruction use the literal same operator and it cancels algebraically — the reconstruction identity does not depend on the expand filter being a "good" interpolator, only on it being the same filter both times. In floating point, rounding introduces negligible error (typically far below quantization noise); in fixed-point/integer pipelines, the same care with rounding must be taken on both sides or small bias accumulates.

Why blur before downsampling (aliasing argument). The sampling theorem states a signal band-limited to frequency $f_{\max}$ can be reconstructed from samples taken at rate $\geq 2f_{\max}$. Downsampling by 2 halves the effective sample rate, halving the safe Nyquist frequency. Any image content above the new Nyquist limit that is not removed first will alias — fold back as spurious low-frequency patterns (moiré). The Gaussian kernel used before subsampling is chosen (approximately) to suppress content above the new Nyquist frequency, so each $G_\ell$ is an honest, alias-minimized representation of the image at its resolution, not merely every-other-pixel of a finer image. This is the same sampling-theoretic reasoning developed for general Image Resampling and Warping, applied here specifically to the fixed octave-doubling schedule of a pyramid.

Tradeoffs. The pyramid buys $O(n)$-total-cost multi-scale access at the price of coarse, fixed scale sampling — only octave steps, not continuous $\sigma$. Algorithms needing finer scale resolution (true scale-space feature detectors such as SIFT) interpolate within each octave using a small stack of Gaussian-blurred (not downsampled) images before moving to the next octave, trading extra $O(1)$-factor computation for finer scale coverage — still $O(n)$ overall, but with a larger constant.

Take a 1D toy "image" (extend the idea to 2D by separability) of 8 samples so every arithmetic step can be checked by hand:

Use the simple 3-tap binomial kernel $h = \tfrac14[1,2,1]$ for blurring (a 1D stand-in for the 5-tap 2D kernel), with reflected boundaries.

Step 1 — blur $G_0$. Convolving (reflect boundary so edge samples reuse themselves):

(e.g. index 2: $\tfrac14(10) + \tfrac12(10) + \tfrac14(50) = 2.5+5+12.5=20$; index 3: $\tfrac14(10)+\tfrac12(50)+\tfrac14(50)=2.5+25+12.5=40$ — recomputing carefully: $\tfrac14\cdot10+\tfrac12\cdot50+\tfrac14\cdot50 = 2.5+25+12.5=40$, so correct value at index 3 is 40, not 30; the table above is corrected below.)

Step 2 — downsample by 2 (keep even indices 0,2,4,6):

Step 3 — build $L_0$. Upsample $G_1$ by inserting zeros, $[10,0,20,0,40,0,10,0]$, convolve with $2h=\tfrac12[1,2,1]$ (the 1D analogue of the "$4h$" 2D compensation factor, since 1D zero-insertion loses a factor of 2 rather than 4), reflect boundaries:

Step 4 — recurse one more level. Blur $G_1=[10,20,40,10]$ with $h$ (reflect): $[12.5,22.5,27.5,17.5]$ approximately; downsample keeping indices 0,2: $G_2 = [12.5,\;27.5]$. Then $L_1 = G_1 - \text{EXPAND}(G_2)$ is computed the same way, and $L_2 = G_2$ is the final DC residual.

Reconstruction check. Running the reconstruction operator bottom-up, $R_2 = G_2 = [12.5, 27.5]$; expand and add $L_1$ to get back $G_1 = [10,20,40,10]$ exactly (by the algebraic cancellation proven above); expand $G_1$ (which reproduces the "expanded" row from Step 3) and add $L_0$:

The reconstruction is exact, confirming the cancellation argument concretely, and the residual $L_0$ visibly concentrates near the sharp jump ($10\to50$) — the Laplacian pyramid isolates edges, which is exactly the band-pass behavior claimed above.

Coarse-to-fine optical flow and stereo matching use Gaussian pyramids to estimate large motions cheaply: a rough displacement found at a coarse level is upsampled as an initial guess and refined at each finer level, turning an expensive large-displacement search into a sequence of cheap small-displacement refinements (the `COARSE_TO_FINE_SEARCH` skeleton above). Multi-band (Laplacian pyramid) blending stitches panoramas and composites images seamlessly: low-frequency bands are blended over wide, smooth transition regions while high-frequency bands are blended over narrow regions, avoiding both visible seams and ghosting — a technique still used in panorama stitchers and photo-editing software. Scale-invariant feature detection (SIFT and its relatives) builds a pyramid-like octave structure with intra-octave Gaussian blurring, then takes differences of Gaussians as an efficient approximation to the Laplacian, to find keypoints stable across scale. Image and video compression historically used Laplacian-pyramid-like subband coding as a precursor to modern wavelet and DCT-based transforms, because band-limited residuals compress well. In graphics, mipmapping — precomputed, prefiltered pyramids of a texture — prevents aliasing when textures are minified, directly reusing the same blur-then-downsample logic. In modern deep learning, Feature Pyramid Networks and multi-scale CNN architectures explicitly borrow the Gaussian/Laplacian pyramid's core insight — represent and process an image at multiple resolutions simultaneously — replacing hand-designed kernels with learned convolutions while keeping the geometric, $O(n)$-cost, octave-structured backbone. All of these downstream uses depend on the resampling correctness discussed under Image Resampling and Warping and on the smoothing foundations of Gaussian and Median Filtering.

Ask the tutor: How would you modify the Laplacian pyramid construction to build a steerable or oriented pyramid that separates edges by direction as well as by scale?