Image Resampling and Warping

A digital image is a discrete grid of samples $I[i,j]$ of an underlying continuous function $f(x,y)$. A geometric transform is a mapping $T: \mathbb{R}^2 \to \mathbb{R}^2$ that sends a coordinate $(x,y)$ in one image to a coordinate $(x', y') = T(x,y)$ in another. Common transforms include translation, rotation, similarity, affine, and projective (homography) maps; the latter three are commonly written as $2 \times 3$ or $3 \times 3$ matrices acting on homogeneous coordinates $(x, y, 1)^T$.

Resampling is the process of producing a new discrete grid of samples of $f$ evaluated at a new set of locations — either because we changed resolution (scaling) or because we composed $f$ with a transform $T$ (warping). Since $f$ is only known at integer grid points, we must reconstruct a continuous approximation $\hat f(x,y) = \sum_{i,j} I[i,j] \, h(x-i, y-j)$ using an interpolation kernel $h$, evaluate $\hat f$ at the required non-integer location, and — if the local scale factor is less than one (minification) — prefilter to remove frequencies that would otherwise alias. There are two algorithmic strategies for applying $T$ to an image. Forward warping iterates over source pixels, computes $T(x,y)$, and scatters (splats) the value into the destination; this can leave holes where no source pixel lands, and can double-write when multiple source pixels map to the same destination cell, requiring splatting kernels or z-buffering for occlusion. Inverse warping iterates over destination pixels $(x', y')$, computes the source location $T^{-1}(x', y')$, and gathers (interpolates) a value from the source; every destination pixel is filled exactly once, with no holes, which is why inverse warping is the standard choice whenever $T^{-1}$ is available in closed form (true for all affine and projective maps).

The inverse-warping algorithm has three cleanly separated stages: compute the inverse transform, loop over destination pixels, and for each one reconstruct a value from the (possibly prefiltered) source image.

Step 1 — Invert the transform. For an affine map $T(\mathbf{x}) = A\mathbf{x} + \mathbf{b}$, the inverse is $T^{-1}(\mathbf{x}') = A^{-1}(\mathbf{x}' - \mathbf{b})$, requiring one 2×2 (or 3×3 homogeneous) matrix inversion done once, not per pixel.

Step 2 — Determine the local scale factor. Compute the Jacobian of $T^{-1}$ at each destination pixel (or, for affine/projective maps, from the singular values of $A^{-1}$). If the minimum singular value is less than 1, the map is shrinking the image locally and a wider reconstruction/prefiltering kernel is needed to avoid aliasing.

Step 3 — Gather and interpolate. For each destination pixel, map back to a source coordinate, take the (possibly scaled) interpolation kernel footprint around that coordinate, and compute a weighted sum of the neighboring source pixels. ``` function INVERSE_WARP(src, H, out_width, out_height, kernel, antialias): H_inv = INVERT(H) # invert 3x3 homography once dst = new image(out_width, out_height) for y_dst in 0 .. out_height - 1: for x_dst in 0 .. out_width - 1:

# map destination pixel center to source coordinates (xs, ys, w) = H_inv * (x_dst + 0.5, y_dst + 0.5, 1) xs = xs / w ys = ys / w if not IN_BOUNDS(src, xs, ys): dst[x_dst, y_dst] = BACKGROUND continue if antialias: scale = LOCAL_MINIFICATION_SCALE(H_inv, x_dst, y_dst) sigma = max(1.0, 1.0 / scale) # widen kernel when shrinking else: sigma = 1.0 dst[x_dst, y_dst] = RECONSTRUCT(src, xs, ys, kernel, sigma) return dst function RECONSTRUCT(src, xs, ys, kernel, sigma):

# gather weighted sum over kernel support, scaled by sigma radius = kernel.support sigma x0 = floor(xs - radius) x1 = ceil(xs + radius) y0 = floor(ys - radius) y1 = ceil(ys + radius) accum = 0.0 weight_sum = 0.0 for j in y0 .. y1: for i in x0 .. x1: wx = kernel.weight((xs - i) / sigma) / sigma wy = kernel.weight((ys - j) / sigma) / sigma w = wx wy accum += w * CLAMP_FETCH(src, i, j) weight_sum += w return accum / weight_sum # normalize (handles boundary truncation) ``` The kernel abstraction covers the standard interpolation families:

• Nearest neighbor: $h(t) = 1$ if $|t| < 0.5$, else $0$. Support radius $0.5$. Zero computation beyond a round, but produces blocky, jagged edges — it reconstructs a piecewise-constant signal, the worst approximation to a smooth image.
• Bilinear: $h(t) = \max(0, 1-|t|)$, a triangle (tent) filter, applied separably in $x$ and $y$. Support radius $1$. This is the convolution of two box filters, equivalent to linear interpolation along each axis, hence "bilinear."
• Bicubic: piecewise cubic polynomial (e.g., the Catmull-Rom or Mitchell-Netravali family), support radius $2$, matches value and first derivative continuity at sample points, producing sharper results with controlled ringing.
• Windowed sinc (e.g., Lanczos): $h(t) = \mathrm{sinc}(t)\cdot \mathrm{sinc}(t/a)$ for $|t| < a$, an apodized approximation to the ideal (infinite-support) sinc reconstruction filter dictated by the sampling theorem. Best frequency response, most expensive, can ring on hard edges.

Antialiasing during minification works by widening the reconstruction kernel (equivalently, prefiltering the continuous reconstructed signal with a low-pass filter matched to the new, coarser sampling rate) before resampling — this is exactly the `sigma` scaling in the pseudocode above. An engineering alternative, used pervasively in GPUs and texture mapping, is mipmapping: precompute a pyramid of progressively half-scale, prefiltered copies of the source image, then at each destination pixel select (or trilinearly blend between) the pyramid level whose sample spacing matches the local minification factor. This turns an expensive wide-kernel convolution into a cheap fixed-size lookup at precomputed resolutions.

Time complexity. For an output image of $N = W_{out} \times H_{out}$ pixels and an interpolation kernel with fixed support radius $r$ (so a $k \times k$ neighborhood with $k = O(r)$), the inverse-warp loop performs $O(N \cdot k^2)$ work: one $O(k^2)$ gather-and-weight computation per destination pixel, plus $O(1)$ amortized cost for the one-time matrix inversion. Nearest neighbor is $k=1$, bilinear $k=2$, bicubic $k=4$, a typical windowed sinc $k \in \{6, 8\}$. This is linear in image size for any fixed kernel, which is optimal — you cannot produce $N$ output pixels in sublinear time. When antialiasing widens the kernel adaptively under strong minification, $k$ can grow proportionally to the minification factor $1/s$, making brute-force prefiltering $O(N/s^2)$ in the worst case (shrinking by $10\times$ in both axes costs $100\times$ the naive per-pixel work). Mipmapping avoids this: building the pyramid costs $O(N_{src})$ total (a geometric series of half-resolution prefiltering passes summing to at most $\tfrac{4}{3}N_{src}$), and each destination lookup afterward is $O(1)$ (or $O(k^2)$ for a fixed small trilinear-blend kernel), independent of the minification factor — trading a bounded amount of extra memory, $O(\tfrac{4}{3} N_{src})$, for a guarantee against worst-case aliasing cost.

Space complexity. A basic inverse warp needs only $O(N)$ space for the output buffer (the source is read-only and can be streamed or memory-mapped). Mipmap-based antialiasing adds $O(N_{src})$ extra space for the pyramid.

Correctness argument. The central claim is: inverse warping with a properly scaled reconstruction kernel produces, for each destination pixel, the value that a band-limited continuous reconstruction of the source, transformed by $T$, would have at that location — up to the kernel's approximation error. The argument has two parts. No holes, no double-writes (coverage). Because the destination loop enumerates every integer pixel $(x', y')$ exactly once and always computes a defined source coordinate via $T^{-1}$ (assuming $T$ is invertible, true for similarity/affine/projective maps with nonzero determinant), every output pixel receives exactly one write. This is an immediate structural invariant of the algorithm — contrast with forward warping, where the image of the integer destination grid under $T$ is not guaranteed to cover all integer source... rather destination... coordinates, since $T$ can be expansive (magnification creates gaps) or contractive (minification creates collisions), forcing hole-filling or splatting heuristics that inverse warping sidemsteps entirely. Approximation quality is governed by the sampling theorem. By the Nyquist–Shannon theorem (from The Fourier View of Images), a signal sampled at rate $f_s$ is exactly recoverable via sinc interpolation only if it is band-limited to frequencies below $f_s/2$. Real images are not perfectly band-limited (edges have infinite-extent spectra), so any finite-support kernel is an approximation to ideal sinc reconstruction, and its Fourier transform (the kernel's frequency response) determines the error: a good kernel has a Fourier transform close to a box (unit gain in-band, zero out-of-band). Nearest neighbor and bilinear have frequency responses that fall off slowly (a $\mathrm{sinc}$ and $\mathrm{sinc}^2$ respectively in the reciprocal domain) and pass significant energy above Nyquist, which manifests as blur (attenuated in-band high frequencies) and some residual aliasing; windowed sinc kernels approximate the ideal box response far more closely at higher computational cost. When the local transform is a minification (scale factor $s<1$), the effective output sampling rate is $s \cdot f_s$, so the reconstructed continuous signal must be low-pass filtered to below $s f_s / 2$ before resampling — skipping this step lets frequencies above the new Nyquist limit fold back (alias) into low frequencies, which is precisely why the `sigma` widening step in the pseudocode is not an optional nicety but a mathematical requirement for correctness under minification.

Tradeoffs, honestly stated. Higher-order kernels reduce blur and aliasing but cost more compute and can introduce ringing (Gibbs-phenomenon overshoot near sharp edges) because no finite-support, non-negative kernel can have an ideal brick-wall frequency response — the Fourier uncertainty principle forbids a filter that is simultaneously compact in space and perfectly band-limited in frequency. Mipmapping is fast and alias-free at coarse scales but blurs anisotropically-foreshortened regions (e.g., a ground plane receding to the horizon) unless combined with anisotropic filtering, which samples multiple mip levels along the direction of maximum compression.

We rotate a small $4\times 4$ grayscale image by $\theta = 30°$ about its center, using inverse warping with bilinear interpolation, and trace one destination pixel by hand. Source image (row $j$, column $i$, values 0–255):

Center of rotation: $(c_x, c_y) = (1.5, 1.5)$. The forward rotation matrix (image coords, $y$ down) is

Since rotations are orthogonal, $R^{-1} = R(-\theta) = R(\theta)^T$. We compute the source location for destination pixel $(x', y') = (2, 1)$:

Step 1 — center the coordinate:

Step 2 — apply $R^{-1} = R(-30°)$, with $\cos(-30°) = 0.8660$, $\sin(-30°) = -0.5$:

Step 3 — uncenter:

So destination pixel $(2,1)$ maps back to source coordinate $(1.683, 0.817)$, which lies inside the image bounds — good, no background fill needed.

Step 4 — bilinear gather. The four surrounding integer samples are at $i \in \{1,2\}$, $j \in \{0,1\}$:

Fractional offsets: $a = x_s - 1 = 0.683$, $b = y_s - 0 = 0.817$. Bilinear interpolation:

So the rotated output pixel $(2,1)$ gets value ≈ 59.5 (rounded to 60 for an 8-bit image). Repeating this for every destination pixel — including checking bounds for corner pixels like $(0,0)$, which will map outside the source and must be filled with a background value or clamped/extrapolated — completes the warp. Note that no source pixel was "moved"; every output value was independently reconstructed, which is exactly why inverse warping guarantees full, gap-free coverage of the destination grid.

A minification check. Suppose instead the same rotation is combined with a uniform scale of $s = 0.3$ (shrinking to 30% size) mapping into a much smaller destination grid. The local Jacobian of $T^{-1}$ now has singular values $\approx 1/0.3 = 3.33$ in the source-to-destination stretch, meaning each destination pixel corresponds to a $3.33 \times 3.33$ block of source pixels. Using plain bilinear interpolation (support radius 1) here would only look at 4 source pixels while ignoring roughly 11 pixels' worth of high-frequency content in between — a textbook aliasing setup. The corrected antialiased version widens $\sigma \approx 3.33$, effectively box- or tent-filtering over the full $3.33\times3.33$ neighborhood before subsampling, which is the discrete analogue of low-pass filtering to the new Nyquist limit.

The theory underlying resampling — that a discrete signal is a set of samples of a continuous one, recoverable exactly only if band-limited — traces to sampling theory developed for communications engineering and formalized rigorously by Claude Shannon, building on earlier work by Harry Nyquist. Its application to two-dimensional images and computer graphics matured through the 1970s–1980s as digital image processing and computer graphics became practical, with foundational texture-mapping and antialiasing techniques (mipmapping, elliptical weighted averaging) developed to make real-time and offline rendering tractable.

Image resampling and warping underlie nearly every geometric operation in a vision or graphics system. Camera calibration and lens-distortion correction resample a distorted image onto an undistorted grid using a per-pixel inverse map. Stereo rectification warps left and right camera images onto a common epipolar geometry before matching. Panorama stitching software (and phone camera "pano" modes) warp overlapping photos onto a shared cylindrical or spherical canvas using exactly the inverse-warp-plus-bilinear (or better) pipeline described here. GPUs implement mipmapped, anisotropically filtered texture sampling in hardware for real-time rendering in games and simulators. Data augmentation pipelines for training convolutional neural networks apply random rotations, scales, and perspective jitter to training images, and sloppy resampling there can silently bake aliasing artifacts into a training set. Medical imaging registration, satellite image orthorectification, and augmented-reality overlay compositing all reduce, at their geometric core, to the same reconstruct-transform-resample loop.

Ask the tutor: How would you design an adaptive-kernel warping routine that automatically balances aliasing suppression against blur for an arbitrary homography?