BookNotesSlidesP-sets Syllabusdraft · v0.2.49.1
◀▲⌂▶
💬Comments welcome. To leave a note, select any text and click the note / highlight button that pops up — or open the panel with the tab at the top-right (‹). Notes are visible only inside our private review group.
💡 In a hurry? Jump to this chapter’s 2 big lessons ↓

6.5 Resampling and upsampling⧉

Resizing an image is routine: shrink a photo for the web, zoom into a crop, straighten a tilted horizon, or undo lens distortion. It feels like it should be a one-liner. To make a picture twice as big you spread the pixels out; to make it smaller you throw some away. The whole point of this chapter is that both of those instincts are wrong, and that doing the job well, well enough that repeated rotations do not degrade the image badly, depends on everything we built in Chapter 6.2, Convolution and Chapter 6.3, Fourier, and on the sampling theory of Chapter 6.4, Sampling, downsampling, and aliasing.

The recorded pixels sample an image already filtered by optics and pixel integration. Resizing asks for values on a new grid. We first reconstruct a continuous approximation from the available samples, then filter it for the new grid and sample it there. Interpolation chooses that approximation; it cannot identify arbitrary unrecorded detail.

6.5.1 Domain operations: moving pixels around⧉

Recall the three families of image operation from Chapter 8.3, Tone point operations, sorted by what each output pixel is allowed to read. Resampling is the domain family: the underlying continuous image is remapped in position; discrete resampling generally blends its recorded colors. Formally a domain transform applies a function $f$ from the plane to itself (from $\mathbb{R}^2$ to $\mathbb{R}^2$) to the image's coordinates, with the rule that if the point $(x,y)$ carried color $c$ in the input, then $f(x,y)$ carries that same color $c$ in the output. The range (the colors) is untouched; only the domain (the coordinates) is remapped.

That one definition spans a huge amount of ground. When $f$ is a uniform scale it is resizing; a rotation straightens a horizon; an affine or perspective $f$ fixes keystoning or stitches a panorama; and a freely chosen $f$ is a warp: picture the image printed on rubber and stretched. The free-form case, including how a user even specifies such an $f$ (and the morphing it unlocks), gets its own Chapter 11.2, Antialiasing for complex transforms. Here we take $f$ as given and concentrate on the hard part shared by every domain transform: the new positions almost never land on the input grid, so we must invent values between the samples.

6.5.2 Start with scaling up⧉

Start with the friendliest case, scaling an image up, say by an integer factor $k$. The output has $k$ times as many pixels along each axis, so most output pixels sit between input pixels and have no value of their own: we have too few samples and must fill in the gaps. (Downsampling is the opposite complaint, too many input samples crowding into each output pixel, and, surprisingly, it needs a slightly different fix. Hold that thought; the machinery is almost identical.) Upsampling is the cleaner place to build intuition, so we begin there.

Figure 6.5.1 follows one small image through every method of this chapter, one stage per method; its first stage is the input itself.

fig-upsampling-stages
Figure 6.5.1, stage 1. One small image followed through every upsampling method of this chapter, one stage per method. Here it is the 32 × 32 input, each pixel drawn as a square. The inputs are two photograph crops (the iguana and the Prudential tower, used elsewhere in this chapter) and four synthetic patterns: letters, a 45° edge with thin lines, a checkerboard, and a zone plate whose rings approach the Nyquist frequency at the edges. The plot underneath gives the luma of one row; click a picture to choose the row. Every stage resamples the sRGB code values, as most editors do; the section on resampling and encoding below says what that changes.

6.5.3 The naive idea, and why it leaves gaps⧉

Here is the instinct nearly everyone has first: loop over the input pixels, and for each one compute where it lands in the output and write its color there. To scale up by $k$, the input pixel at $(x,y)$ goes to output location $(k\,x,\ k\,y)$:

$$ \text{out}[\,k\,y,\ k\,x\,] \;\leftarrow\; I[\,y,\ x\,]. $$

This is forward mapping: push each input pixel forward to its destination. It is intuitive, and it is broken. Scale up by $k=2$ and you write output pixels $(0,0),(2,0),(4,0),\dots$ while every odd-coordinate pixel is never touched: a regular lattice of black holes (Figure 6.5.2, left). The output is mostly gaps, because spreading $N$ input pixels over $k^2 N$ output positions can fill at most $N$ of them: one cell of every $2\times2$ block, three left empty.

You might try to patch it with an inner loop that paints the little square between four placed pixels. That limps along for an axis-aligned enlargement, but the instant $f$ is a rotation the gaps stop being neat squares: they become slivers at odd angles, and "fill the gap" no longer has an easy shape. Forward mapping fights you. We need a different loop.

fig-resample-forward-inverse
Figure 6.5.2. Forward vs. inverse mapping on a 2× upscale of a small rainbow grid. Left (forward): loop over input pixels and push each to $(2i,2j)$; only one cell of every 2×2 output block gets written, leaving the other three as black gaps (a regular lattice of holes; worse under a rotation, where the holes become tilted slivers). Right (inverse): loop over output pixels and, for each, apply $f^{-1}$ to reach back into the input and sample, so every pixel is filled, no gaps. The arrows trace the two routes for a few pixels: forward arrows start at an input pixel and end where it is written, inverse arrows start at an output pixel and end at the input pixel it reads. The whole chapter loops over the output and reaches back into the input.

The second stage of Figure 6.5.1 runs this naive loop on a real image.

fig-upsampling-stages
Figure 6.5.1, stage 2. Naive forward mapping: each input pixel is written to the one output pixel under its center, at $k$ times its coordinates, and nothing else is written. Each $k\times k$ block of the output gets one pixel and the rest stay black, so at ×4 the 128 × 128 output has 1,024 written pixels and 15,360 holes; the title counts them for any factor. In the row plot the red dots are the only values the output row has. From this stage on the output is drawn to scale: an output pixel takes the same room on the screen whatever the factor, so ×2 comes out a quarter as wide as ×8, and the input appears at that same scale at the far left, beside its close-up. "Fill the frame" stretches every factor to the same size, which shows the detail and hides the factor, because every factor samples the same reconstruction, only more or less densely; the row plot shows that density as one dot per output pixel.

6.5.4 Loop over the output, use the inverse transform⧉

The fix is to turn the loop inside out. Instead of pushing inputs forward, iterate over the output pixels (all of them, so none can be skipped) and for each one ask where did this come from in the input? That means applying the inverse transform $f^{-1}$ to the output coordinate and reading the input there:

$$ \text{out}[\,y,\ x\,] \;\leftarrow\; I\big[\,f^{-1}(x,y)\,\big]. $$

For a scale-up by $k$ the inverse divides by $k$, so output pixel $(x,y)$ reads the input at $(x/k,\ y/k)$. Now every output pixel receives exactly one value and there are no gaps by construction. This is inverse mapping: loop over output pixels and use the inverse transform when it exists. An occlusion, a fold, or a location outside the source also needs an explicit visibility or boundary rule. It is the same "where-from, not where-to" reflex as the flip in convolution (to fill a destination, reach back into the source), and it is the single most important habit in this chapter. (Figure 6.5.2 sets the two loops side by side.)

💡 Big lesson (L6.19) — resample by gathering, rarely by scattering

Loop over output pixels and pull back through $f^{-1}$ to find where each one reads from. Gathering fills every output exactly once; scattering pushes input pixels forward and leaves holes whose shape depends on the transform, a lattice of gaps under a scale and tilted slivers under a rotation. The same reflex governs warping, morphing, texture mapping, and flow-based reprojection later. Scattering earns its place only when the forward map is the thing you actually have, as in splatting a point cloud, and then it needs an explicit accumulation and normalization pass.

There is one immediate wrinkle, and it is the real subject of everything that follows. Inverse mapping is clean, but $f^{-1}(x,y)$ lands at a fractional coordinate: $(x/k,y/k)$ is rarely an integer, and under a rotation it essentially never is. There is no pixel sitting at $(1.3, 4.7)$. So the genuine content of resampling is the question we have been circling: given the input samples, what value should we report at a position between them? Every method below is one answer.

6.5.5 Nearest neighbor⧉

Round the fractional coordinate to the nearest sample and copy its value. At $1.3$, return $I[1]$; at $1.8$, return $I[2]$. Truncating the coordinate instead gives a different, biased convention. Nearest neighbor is appropriate for categorical labels and for deliberately showing the original pixels, but looks blocky on photographs.

It works and it is fast, but the result is conspicuously blocky: a whole range of fractional positions all snap to the same input pixel, so an upscaled image looks like fat rectangular tiles (Figure 6.4.16, nearest panel), each input pixel becoming a hard $k\times k$ block. Occasionally that is precisely what you want (pixel art, or a debug view where you want to see the literal samples), but as a general-purpose resize it is the floor, not the goal. As the challenge for nearest neighbor: try upscaling a face by $8\times$ this way and watch it turn to mosaic. We do much better by interpolating between samples instead of snapping to one.

The third stage of Figure 6.5.1 replaces that scatter with a gather and rounds.

fig-upsampling-stages
Figure 6.5.1, stage 3. Nearest neighbor, by inverse mapping: every output pixel reads the input sample nearest to the position it maps back to, so each input pixel becomes a $k\times k$ block and nothing is left empty. The kernel plot shows the box this amounts to; the row plot is a staircase through the samples.

6.5.6 Linear interpolation, in 1-D first⧉

Drop to one dimension, where the idea is cleanest. An image row is a 1-D array of values, and we want the value at $1.3$, between the samples at $1$ and $2$. The natural move is to blend the two neighbors in proportion to how close the query sits to each. At $1.3$ we are $0.3$ of the way from $1$ toward $2$, so

$$ I[1.3] \;=\; 0.7\cdot I[1] \;+\; 0.3\cdot I[2]. $$

This is linear interpolation, and there is one detail that trips up everyone at first: the weights run opposite to the distances. The neighbor you are $0.3$ away from gets weight $0.7$, not $0.3$, because you want weight $1$ when you sit exactly on a sample and $0$ when you are a full step away. Say it back in words: take a weighted average of the two bracketing samples, weighting each by one minus the distance to it. Geometrically, the reconstructed signal is just the straight line connecting consecutive samples (Figure 6.5.3).

fig-linear-interp-1d
Figure 6.5.3. 1-D linear interpolation. To reconstruct the value at the fractional position $1.3$ between samples $I[1]$ and $I[2]$, blend them in proportion to one minus the distance to each: $I[1.3]=0.7\cdot I[1]+0.3\cdot I[2]$. The reconstructed signal is the straight line joining consecutive samples; a neighbor's weight is $1$ when you sit on it and falls linearly to $0$ a full sample away, the opposite of the raw distance.

6.5.7 Bilinear, in 2-D⧉

In two dimensions we now have four surrounding samples (the corners of the unit square the query falls inside) and two fractional coordinates to interpolate against. The clean construction is to do 1-D linear interpolation twice: first along $x$, then along $y$ (Figure 6.5.4). Blend the top pair of pixels horizontally to get one value, blend the bottom pair horizontally to get another, then blend those two results vertically with the $y$ fraction. The four corner weights come out as products of the $x$- and $y$-blend factors, which is exactly why this is bilinear: linear in each axis, but their product, not a single plane through all four points.

That last point is a real subtlety. Four arbitrary corner values do not in general lie on a flat plane, so a "purely linear" 2-D surface through all four simply does not exist. Bilinear sidesteps this by being linear separately in each direction (the product of two tents), which is well-defined and cheap. The resulting surface is doubly ruled: dead straight along every horizontal and every vertical line, yet curved overall, twisting into a saddle when the four corners disagree (Figure 6.5.5). The only way to get a truly planar piecewise-linear reconstruction is to triangulate: cut each square into two triangles and interpolate within each, since three corners do fix a plane.

fig-bilinear
Figure 6.5.4. Bilinear interpolation as two 1-D passes, in both orders (interactive). The query $(s,t)$ falls inside the unit square spanned by its four nearest input pixels, $c_{00}$ to $c_{11}$; each square is shaded with the interpolated values. On the left we interpolate along $x$ first, on the bottom and top edges (green), $f(s,0)=(1-s)\,c_{00}+s\,c_{10}$ and $f(s,1)=(1-s)\,c_{01}+s\,c_{11}$, then along $y$ between those two (orange). On the right we interpolate along $y$ first, on the left and right edges, $f(0,t)$ and $f(1,t)$, then along $x$. The intermediate values differ, the result does not: both orders expand to $(1-s)(1-t)\,c_{00}+s(1-t)\,c_{10}+(1-s)t\,c_{01}+st\,c_{11}$, four weights that are products of the $x$- and $y$-blend factors, linear in each axis, hence bilinear. Drag the query point in either square, or set the corner values; the defaults are those of the height field that follows.
fig-bilinear-heightfield
Figure 6.5.5. The bilinear surface over one unit cell, drawn as a height field whose four corner samples are the controlling heights. Every line of constant $x$ and every line of constant $y$ is dead straight (the sheet is ruled in both directions), yet the surface as a whole bows into a saddle, proof that bilinear is linear in each axis separately but not one flat plane through the four corners. At a query point $(s, t)$, marked with its construction lines on the unit square under the sheet, the figure draws the two 1-D passes that compute it: interpolate along $x$ on the two edges $t=0$ and $t=1$, then along $y$ between the results, or along $y$ first and $x$ second; both orders reach the same point of the sheet. Drag to orbit; toggles add the two diagonals, which curve unless the corners are coplanar, and shading.

The next stage of Figure 6.5.1 applies bilinear interpolation to the whole image.

fig-upsampling-stages
Figure 6.5.1, stage 4. Bilinear: each output pixel blends its four nearest samples with the tent kernel. The blocks are gone and the edges are soft. In the row plot the curve joins the samples with straight segments and kinks at each.

Choosing the diagonal of each cell from the data is data-dependent triangulation, introduced for image reconstruction by Yu, Morse and Sederberg (2001); the per-cell rule of the sidebar is the pixel-level version of Su and Willis (2004), and the fifth stage of Figure 6.5.1 runs it after bilinear.

fig-upsampling-stages
Figure 6.5.1, stage 5. Data-dependent triangles: each cell of four samples, $a\ b$ above $c\ d$, is cut along the diagonal whose two end samples differ less in luma, and each triangle is interpolated linearly. The counts under the input give how many cells took each cut; the toggle draws the chosen diagonals on the input, leaving out the ties, which are cut a–d. On the 45° edge every cell the edge crosses is cut along it, and the edge stays sharper than with bilinear. A one-pixel line at 45° is the method's own failure: in every cell it crosses, both diagonals join equal samples (line to line, background to background), so the tie-break decides. The line that runs along the tie-break diagonal comes out smooth, and the one that runs across it breaks into beads. Method: Su and Willis (2004), after Yu, Morse and Sederberg (2001).
Exercise

Show that interpolating $y$ first, then $x$ gives exactly the same result as $x$ first, then $y$. (Write both as weighted sums of the four corners and compare the coefficients, the order of the two separable passes does not matter.)

Putting the two halves together (the inverse-mapping loop and the bilinear sample) gives the whole upsampler. Loop over output pixels, send each back through $f^{-1}$ to a fractional input location, split that into its integer cell and fractional offsets, and blend the four corners:

Algorithm 6.5.1
for each output pixel (x, y):
    (u, v) = f⁻¹(x, y)                  # fractional input location
    x0, y0 = floor(u), floor(v)
    s, t = u - x0, v - y0               # fractional offsets in [0, 1)
    top = (1 - s)·I[x0, y0]   + s·I[x0+1, y0]
    bot = (1 - s)·I[x0, y0+1] + s·I[x0+1, y0+1]
    out[x, y] = (1 - t)·top + t·bot     # blend the two rows vertically

Read it back: the two horizontal blends collapse the top and bottom edges of the cell, and the final vertical blend combines them: 1-D linear interpolation done twice. Bilinear is the workhorse default: cheap, no blocky tiles, smooth enough for most resizing, and the reconstruction you get for free from a graphics-hardware texture lookup (Figure 6.4.16, bilinear panel). But it has a ceiling. The reconstructed surface has creases (it is only continuous, not smooth, kinking at every sample) and on a strong upscale it reads a touch soft. To push past it we need to look at more than four neighbors, and the "interpolate between neighbors" story does not generalize gracefully. We need a better lens, but first a word on what any of this costs.

6.5.8 The convolution perspective⧉

Here is the more general formulation. So far we have described reconstruction as interpolating between neighbors. That picture is fine for linear interpolation and for upsampling, but it has two failings: it does not extend cleanly to fancier reconstructions, and it does not extend to downsampling at all. Since this whole chapter is fundamentally about (re)sampling, the Fourier-flavored view from Chapter 6.4, Sampling, downsampling, and aliasing has much better insight to offer. In practice that view means one thing: reconstruction is convolution.

Concretely, the output value at a fractional location is a weighted sum of nearby input pixels, with the weights read off a smooth analytical kernel $k$ centered at the desired (non-integer) location. Slide a continuous kernel so its center sits exactly on the query point $f^{-1}(x,y)$, evaluate the kernel at each surrounding integer input pixel to get that pixel's weight, and sum. The general resampling formula is

$$ \text{out}(x,y) \;=\; \sum_{x'} I(x')\,\cdot\, k\big(f^{-1}(x,y) - x'\big), $$

summing over input pixels $x'$ in the neighborhood, the kernel evaluated at the offset between the query and each pixel. In words: to find an output pixel, inverse-transform it to a fractional input location, drop a smooth kernel centered there, and accumulate kernel-weight times pixel-value over the nearby samples. It is convolution, but where convolution in the earlier chapter lived purely on the integer grid, here the kernel is a continuous function we can evaluate at any fractional offset (ask for $k(-0.7)$ and it answers), because the center is non-integer.

The resampler needs weights for fractional offsets. A formula for a real-valued argument supplies them, as do a lookup table or a bank of filters for quantized fractional phases. A list of integer-offset convolution taps alone does not specify interpolation between taps. Nearest neighbor’s box kernel is discontinuous, so continuity is not a requirement.

fig-resample-kernel
Figure 6.5.6. Why the kernel must be a continuous function. A row of discrete input samples sits on the integer grid; a smooth kernel $k(x)$ is centered at a non-integer query location $x^\*$ between two of them, and the output is the sum of each sample times the kernel's height at that sample's offset. Because the center is fractional, those offsets are fractional too, so the kernel has to be a closed-form curve readable at any real argument, not a fixed stencil of integer taps.

Bilinear is already a special case. The linear-interpolation weights $0.7$ and $0.3$ are exactly the values of a triangular (tent) kernel, one that is $1$ at its center and falls linearly to $0$ one sample away on each side, evaluated at the offsets to the two neighbors. So "interpolate linearly" and "convolve with a tent" are the same operation, seen two ways. Nearest neighbor is the same story with a box kernel (a flat $1$ over a one-pixel window). Once reconstruction is "convolve with a kernel," choosing a better reconstruction is just choosing a better kernel.

Figure 6.5.7 draws these kernels as 2-D kernels $k(x)\,k(y)$, each with its profile, an upsampling made with it, and its Fourier transform, starting with the box and the tent.

fig-resample-kernel-stages
Figure 6.5.7, stage 1. The box kernel of nearest neighbor, $k(x)=1$ for $|x|<\tfrac12$, as the separable 2-D kernel $k(x)\,k(y)$: a column one sample wide and of height 1 over the grid of sample positions. Drag the height field to turn it. On the right, top to bottom: the profile $k(x)$, with a dot wherever the kernel changes sign; eleven pixels of one row of the input (their luma, marked red on the input below) reconstructed as a sum of copies of the kernel, each centered on a sample and scaled by its value (blue); the magnitude of the Fourier transform along $f_y=0$ against the ideal reconstruction filter (green), which passes every frequency up to Nyquist, half a cycle per sample, and none beyond; and the 2-D magnitude $|K(f_x)\,K(f_y)|$ on a log scale, with the Nyquist square outlined. Bottom left, the input, a 16 × 16 window of it, and that window upsampled ×8 with the same kernel, shown at 1:1. The menu picks the input, the same six as Figure 6.5.1, the iguana first. The box's transform is $\operatorname{sinc}(f)$: already down to 0.64 at Nyquist, it leaves 23% of its energy beyond Nyquist in 1-D and 40% outside the square in 2-D, the aliasing that shows as blocks.
fig-resample-kernel-stages
Figure 6.5.7, stage 2. The tent, $k(x)=1-|x|$ for $|x|<1$, in 2-D. The product $k(x)\,k(y)$ is not a pyramid: along a diagonal it falls as $(1-t)^2$, and its four taps are the four bilinear weights. The reconstruction of the row is the polyline through the samples. Its transform is $\operatorname{sinc}^2(f)$, 0.41 at Nyquist, so the passband droops more than the box's, while only 5% of the energy lies beyond Nyquist in 1-D and 10% outside the square in 2-D.

The box and the tent are the crudest kernels. Better ones, the Mitchell–Netravali cubics $k_{BC}$ and the Lanczos windowed sinc, are developed as approximations of the ideal sinc in Chapter 6.4, Sampling, downsampling, and aliasing, where they first serve as prefilters for downsampling. As reconstruction kernels they upsample with less blur than the tent and without the box's blocks, at the price of some overshoot near edges.

The third stage of Figure 6.5.7 draws any member of the family in 2-D, with its Fourier transform.

fig-resample-kernel-stages
Figure 6.5.7, stage 3. The Mitchell–Netravali cubics $k_{BC}$ in 2-D. The menu sets the named members and the sliders any $(B,C)$. Below zero the sheet turns blue: Catmull–Rom's negative lobes are four shallow troughs along the axes, down to −0.07, and on the diagonals, where both factors are negative, the product is positive again. Catmull–Rom's reconstruction passes through the samples and overshoots at sharp transitions; the cubic B-spline has no negative lobe, does not overshoot, and misses the samples. The B-spline is the box convolved with itself four times, so its transform is $\operatorname{sinc}^4(f)$: 0.16 at Nyquist with 0.6% of the energy beyond it, against 0.49 and 4.7% for Catmull–Rom and 0.38 and 3.2% for Mitchell. For every $(B,C)$ the transform is zero at every nonzero integer frequency, because the weights sum to one at every phase.

The sixth stage of Figure 6.5.1 upsamples with any member of this family.

fig-upsampling-stages
Figure 6.5.1, stage 6. The Mitchell–Netravali cubics. The menu sets the named members (Catmull–Rom, which is Keys' cubic convolution with $a=-\tfrac12$; Mitchell; the cubic B-spline; and Keys with $a=-\tfrac34$, the cubic OpenCV uses) and the sliders set any $(B,C)$. The readouts under the kernel give $k(0)=1-B/3$ and $k(\pm1)=B/6$, so only $B=0$ passes through the samples; with the B-spline the row curve misses the dots. Negative lobes take the output outside the input's range, printed above the row plot: on the checkerboard (squares of 0.1 and 0.9, at ×4) Catmull–Rom reaches 1.03, Mitchell 0.96, and the B-spline stays at 0.9.

Its last stage swaps the cubic for the Lanczos kernel, with order $a$ from 2 to 4.

fig-upsampling-stages
Figure 6.5.1, stage 7. Lanczos of order $a$ from 2 to 4: the sinc windowed by a sinc $a$ times wider, $2a$ taps per axis, with the weights renormalized to sum to one. More lobes, more ringing: on the checkerboard at ×4 the output reaches 1.04, 1.11 and 1.14 for $a=2$, 3 and 4.

The last stage of Figure 6.5.7 draws the Lanczos kernel the same way, for orders 2 to 4.

fig-resample-kernel-stages
Figure 6.5.7, stage 4. The Lanczos kernel $\operatorname{sinc}(x)\operatorname{sinc}(x/a)$ in 2-D, order $a$ from 2 to 4: $2a$ taps per axis, $4a^2$ in 2-D, and $a-1$ sign changes on each side. Its transform is the closest of the four kernels to the ideal: about 0.50 at Nyquist for every order, with 4.6%, 2.9% and 2.1% of the energy beyond Nyquist for $a=2$, 3 and 4 (9%, 6% and 4% outside the square in 2-D), and the step from passband to stopband steepens as $a$ grows. The price is the lobes: the 2-D kernel dips to −0.09, −0.15 and −0.18, and the reconstruction rings at sharp transitions. The weights are not renormalized here, and Lanczos weights do not sum exactly to one, so a constant row comes back with a ripple of up to 1.9%, 0.6% and 0.2%; practical resamplers renormalize them.

6.5.9 The full signal-processing view: sample, reconstruct, prefilter, resample⧉

The kernels above are one part of a pipeline that explains every resampler in this chapter. We write it in one dimension, with input samples $I(j)$ at the integers and output samples at a spacing $q$ measured in input pixels, so $q>1$ shrinks the image and $q<1$ enlarges it. Two dimensions repeat it along each axis, and a general warp replaces the output position $mq$ by $f^{-1}$ of the output pixel.

The sampled representation. The input is a list of numbers, not a function. We represent it as a train of impulses, $I_\delta(x)=\sum_j I(j)\,\delta(x-j)$. In the frequency domain, sampling at unit spacing makes the spectrum periodic, $\hat I_\delta(u)=\sum_n \hat I(u-n)$: copies of the spectrum of the underlying continuous image, one every cycle per pixel (Chapter 6.4, Sampling, downsampling, and aliasing). The band $|u|<\tfrac12$ holds the information; the copies are what sampling added.

Figure 6.5.8 follows these steps on one row of an image, one column per step, each signal above its spectrum.

fig-resampling-pipeline-1d
Figure 6.5.8, stage 1. Resampling in one dimension, one step at a time, on a row of the Prudential facade whose windows repeat every 3.7 pixels. This stage shows the input alone: the samples $I(j)$ as impulses, with the row's mean removed (top), and the magnitude of their transform $\hat I_\delta(u)$ (bottom), periodic with a period of one cycle per pixel, with the band $|u|<\tfrac12$ shaded. The facade's repetition is the pair of peaks at ±0.27 cycles per pixel, repeated around every integer frequency. The menus and sliders choose another photograph, row or window, or one of four synthetic inputs sampled at the integers (sine, box, sawtooth, box train); a periodic input has a line spectrum, drawn as arrows.

Reconstruction. Convolving the impulses with a kernel $k$ gives a continuous image, $r(x)=\sum_j I(j)\,k(x-j)$, the formula of the previous section. In frequency the convolution becomes a product, $\hat r(u)=\hat k(u)\,\hat I_\delta(u)$. The ideal $k$ is the sinc, whose transform keeps the central copy and removes the others. A practical kernel does both jobs imperfectly: it attenuates the central copy near $u=\tfrac12$, which blurs, and it keeps part of the neighboring copies, which shows as blocks or staircases. This is imaging, the reconstruction counterpart of aliasing.

The second stage of Figure 6.5.8 adds the reconstruction.

fig-resampling-pipeline-1d
Figure 6.5.8, stage 2. Reconstruction: $r=I_\delta*k$ (top, over the faint samples) and its spectrum $\hat r=\hat k\,\hat I_\delta$ (bottom), with $|\hat k|$ in gray. The buttons choose $k$: nearest neighbor, linear, Keys' cubic or Lanczos-3. Each kernel lowers the central copy toward $u=\tfrac12$, which blurs, and keeps part of the copies around $u=\pm1$, the imaging. At the facade's 0.27 cycles per pixel, $\hat k$ is 0.89 for nearest neighbor, 0.78 for linear, 0.92 for the cubic and 1.00 for Lanczos-3; at its first image, 0.73 cycles per pixel, it is 0.33, 0.11, 0.08 and 0.004.

Prefiltering. If the output grid is coarser than the input, $r$ still holds frequencies up to half a cycle per input pixel, more than the output can represent. A low-pass filter scaled to the output spacing, $h_q(x)=\tfrac1q\,h(x/q)$, removes them: $g=r*h_q$, and $\hat g(u)=\hat h(qu)\,\hat r(u)$, ideally a box that keeps $|u|<1/(2q)$. When enlarging, $q<1$ and that band is wider than what reconstruction already left, so the prefilter has nothing to do.

The third stage adds the prefilter.

fig-resampling-pipeline-1d
Figure 6.5.8, stage 3. The prefilter $h_q$, stretched to the output spacing $q$: $g=r*h_q$ (top, over a faint $r$) and $\hat g=\hat h(qu)\,\hat r$ (bottom), with $\hat h(qu)$ in gray. The slider sets $q$ from 1/4 (enlarging four times) to 4 (shrinking four times); here $q=3$. The shaded band is now the output's, $|u|<1/(2q)$. The facade's peaks lie outside it, so the Gaussian removes them and leaves an almost flat $g$. Below $q=1$ the output band covers the whole frequency window and the prefilter has nothing left to remove.

Sampling on the new grid. The output is $O(m)=g(mq)$. Sampling at spacing $q$ makes the output spectrum periodic again, with copies every $1/q$; whatever $\hat g$ still holds beyond $1/(2q)$ folds back into the band as aliasing.

The fourth stage samples the new grid.

fig-resampling-pipeline-1d
Figure 6.5.8, stage 4. Sampling at spacing $q$: the outputs $O(m)=g(mq)$ (top, over a faint $g$) and their transform (bottom), periodic with period $1/q$, each copy $\tfrac1q\,\hat g(u-k/q)$ in its own color over their faint sum. This stage opens with the prefilter off. At $q=3$ the copies of the facade's peaks then land at ±0.064 cycles per pixel, inside the output band: a period of 16 pixels that the facade does not have. With linear reconstruction, the aliased part is 9% of the output's energy in the band after the Gaussian prefilter, 2% after Lanczos-3, and 88% after a box of width $q$, the area average.

Reconstruction and prefiltering are both convolutions, so they combine into one kernel, $\rho=k*h_q$, and the whole pipeline collapses into a single weighted sum per output pixel:

$$O(m)=\sum_j I(j)\,\rho(mq-j).$$

When enlarging, $\rho$ is essentially the reconstruction kernel $k$, which keeps its width in input pixels. When shrinking, it is essentially the prefilter $h_q$, whose width grows with the reduction. This is the scaled kernel of the next section, whose normalized formula is this sum with $\rho$ replaced by the stretched kernel (Figure 6.5.9). Figure 6.4.9 shows the same steps for the first sampling of a continuous signal.

fig-reconstruction-pipeline
Figure 6.5.9. The full resampling pipeline as one chain: sampled input → reconstruct (continuous signal) → prefilter (low-pass for the output grid) → resample (new grid). The reconstruction and prefilter steps are both convolutions, so they fold into one combined kernel; upsampling is dominated by the reconstruction filter, downsampling by the prefilter.

The last stage of Figure 6.5.8 draws $\rho$ beside the pipeline.

fig-resampling-pipeline-1d
Figure 6.5.8, stage 5. The combined kernel $\rho=k*h_q$ (right, black), with $k$ and $h_q$ in gray at their own heights, and its transform $\hat\rho(u)=\hat k(u)\,\hat h(qu)$; the shaded band is the ideal one, $|u|<\min(\tfrac12,1/(2q))$. Move $q$. With linear $k$ and a Gaussian $h$, $\rho$ is within 3% of $k$ at $q=1/4$, with a support of ±1.5 pixels, and within 2% of $h_q$ at $q=4$, with a support of ±9 pixels (both distances in $L^1$). The outputs drawn in the fourth column are the single sum $O(m)=\sum_j I(j)\,\rho(mq-j)$.

6.5.10 Upsampling vs. downsampling: scale the kernel⧉

We have promised twice that downsampling reuses the same machinery; now we cash it in. Whether you scale up or down, you run the same convolution with the same family of kernels. The only thing that changes is how you scale the kernel.

When upsampling, the kernel keeps a fixed footprint in the input image (about four pixels wide for bicubic, however much you enlarge) because the input samples are already dense enough; you are just reading between them. So for upsampling the kernel width in input space is constant (a useful corollary: clamp the kernel's scale factor at $1$, never narrower).

Downsampling is the harder case, and not by a little. Pick one input pixel per output pixel (nearest neighbor, no blur) and you get aliasing: a fine texture turns into bold, wrong, low-frequency moiré (a brick wall or a striped shirt shot from far away dissolving into wide phantom bands; Figure 6.5.13). This is not a cosmetic glitch but the wagon-wheel / Prudential-tower effect from Chapter 6.4, Sampling, downsampling, and aliasing: detail finer than the output grid can represent does not vanish, it folds down and masquerades as coarse structure.

Widen the filter in proportion to the reduction. In one dimension, let $q=W_{\mathrm{in}}/W_{\mathrm{out}}$ and $a_s=\max(1,q)$. A practical normalized resampler uses

$$w_j=k((u-j)/a_s),\qquad O(x)=\frac{\sum_j w_j I(j)}{\sum_j w_j}.$$

The source coordinate is $u$. The sum covers the stretched kernel’s support. A fivefold reduction stretches a four-sample cubic support to twenty input samples, not five. Normalization preserves a constant image when the denominator is nonzero. Finite support still permits residual aliasing. How much depends on the kernel as well as its width. A box as wide as the reduction, the plain average of each block, is the most common prefilter and one of the leakiest: the sidelobes of its transform pass a fine periodic texture. Shrinking the facade of Figure 6.5.8 threefold, 88% of the energy that ends up in the output band is aliased with a box prefilter, against 9% with a Gaussian and 2% with Lanczos-3 (stage prefilter).

Pixel-center conventions also matter. With centers at integer indices and image edges at half-integers, use $u=(x+1/2)W_{\mathrm{in}}/W_{\mathrm{out}}-1/2$. This differs from aligning the first and last centers. Choose one convention and use it consistently when composing transforms.

fig-bicubic-bc
Figure 6.5.10. The Mitchell–Netravali cubic family. A single pair of parameters $(b, c)$ tunes the cubic resampling kernel across a whole family: from blurry (the cubic B-spline, $b=1$) through the balanced Mitchell ($b=c=\tfrac13$) to sharp-but-ringing (Catmull-Rom, $b=0$, $c=\tfrac12$). Drag the $(b,c)$ point to trade blur against ringing and overshoot, shown on the kernel, its frequency response, and a resampled edge.
fig-kernel-scaling
Figure 6.5.11. The one knob that distinguishes up- from downsampling. Upsampling: the kernel keeps its native width in input space (about four pixels for bicubic), reaching across the same handful of dense input samples no matter how much you enlarge. Downsampling: the kernel is stretched to span the output-pixel spacing, by $s\times$ for an $s\times$ reduction, so it averages over all the input detail crammed into each new pixel, doubling as the anti-aliasing prefilter. Same kernel shape, scaled differently.

Not every resize function widens it. Parmar et al. found that several widely used libraries reduce an image with a kernel of fixed width by default; Figure 6.5.12 puts that behavior beside the alternatives and traces the aliasing to the kernel.

fig-resize-library-behaviors
Figure 6.5.12. Four downsamplers on one input, side by side, each with the kernel it applies drawn under it in source pixels: the dots are the weights one output pixel takes from the source pixels (the ticks), and the shaded band is that output pixel's footprint, $s$ source pixels wide. Decimation reads one source pixel. A fixed-width kernel, bilinear or bicubic at the source's own scale, reads two to four pixels near the middle of the footprint and ignores the rest, the default of several of the libraries Parmar et al. compare on a thin circle. The kernel widened by the factor covers the footprint, and so does the area average. On the thin circle the first two break the line into scattered dots; on the zone plate they bring back rings that are not in the source. The readout under each result measures its aliasing: the part of the source above the output's Nyquist frequency (what an ideal low-pass in the DCT basis removes) is passed through the same downsampler, and whatever comes out is content the output cannot hold. It is given as an RMS in 8-bit code values and as a share of that part's energy. At 1/4 on the circle with bicubic, decimation keeps 98% of it, the fixed kernel 68%, the area average 8.6% and the widened kernel 2.3%. At odd whole factors the fixed tent lands on a single source pixel and becomes decimation; at other factors its taps shift from one output pixel to the next. The readout measures aliasing, not blur: the widened tent scores lowest and is also the softest. The arithmetic is on sRGB code values, as the libraries do, and each result is enlarged by a whole number with nearest neighbor. Peacock photograph by Bernard Spragg (CC0).

The effect is ugliest on real, regular structure: the rows of windows on a glass facade, shrunk for a thumbnail, break into sweeping moiré if you decimate without prefiltering (Figure 6.5.13).

fig-downsample-comparison
Figure 6.5.13. Downsampling a real photograph of a building facade. Left: naive decimation (pick one pixel per output) turns the regular grid of windows and mullions into coarse, wrong moiré bands sweeping across the facade: aliasing on actual content, not a test pattern. Right: the same shrink with a properly scaled prefilter first, a clean, sharp miniature with the window grid faithfully averaged down. The artifact is exactly the aliasing of the previous figure, now on the kind of image you would really downsample.
💡 Big lesson (L6.20) — every resampling pass multiplies the spectrum again, so resample once

Repeated resampling accumulates blur, ringing, and sampling error. For repeated fixed convolution the frequency response is exactly $H^N$; for rotations it is only an illustration because the grid and effective filter change at each pass. Compose geometric transforms and resample the original once.

fig-xkcd-2671
Figure 6.5.14. xkcd #2671 "Rotation", Randall Munroe, xkcd.com/2671 (CC BY-NC 2.5).
fig-rotation-challenge
Figure 6.5.15. Andrew Adams' rotation challenge, and how much the resampling filter matters. Each image is rotated a full $360°$ in $N=10$, $60$, or $360$ equal steps, every step resampling the previous result; it is the same total rotation, just chopped into more or fewer passes, so more (smaller) steps compound more damage. The three rows repeat the experiment with a different resampling filter: nearest-neighbor, bilinear, and bicubic (cubic spline). Nearest jitters and loses registration; bilinear (a mild low-pass) washes the iguana's fine scales into mush by $N=360$; bicubic, whose negative lobes restore some of the contrast each pass would otherwise lose, holds the detail markedly longer. A validity mask blacks out the corners that ever rotate out of frame, so the kept region shrinks toward the inscribed disc. The lesson: the resampling kernel is not a detail: it is the difference between a tool you can trust and one that quietly destroys your image.

Run the same experiment yourself in the interactive version below: one image, your choice of total angle and step size (the total is rounded to a whole number of steps), and any of five filters, and when you pick bicubic, the same $(B, C)$ handle from above to dial the cubic between blur and ring.

For a far deeper practical treatment of resampling filters, their windows, and the subtle quality differences between them, Nicolas Robidoux's analysis for ImageMagick is the reference worth reading (<https://usage.imagemagick.org/filter/nicolas/>); it is the kind of care that separates a library taking resampling seriously from one that ships a default and hopes.

fig-rotation-resample-interactive
Figure 6.5.16. The rotation challenge, interactive. Rotate one image to a chosen total angle in equal steps, each step re-sampling the previous result with the selected filter (nearest · box · bilinear · Lanczos · bicubic). The rotation runs on sRGB code values, and "box" here has two-pixel support, unlike nearest-neighbor. The orientation depends only on the total; the damage depends on the number of passes $N$ and on the filter: drop the step size to compound more passes, switch filters to watch nearest jitter, box and bilinear blur to mush, and Lanczos or a well-tuned bicubic hold the detail. For bicubic, drag the $(B,C)$ point to choose the cubic. Set the total to $360°$ to land back at the original orientation and read off pure resampling loss. A bottom panel plots the kernel's 1-D frequency response $H(f)$ and, as an approximation of the $N$ stacked passes, $H(f)^N$: the curve collapsing toward DC is why detail vanishes: each pass multiplies the spectrum by $H$, so after $N$ passes the surviving spectrum is roughly $H^N$. The plot keeps one scale, 0 to 2, and shades orange where the passes boost a frequency (the sharpening that rings) and blue where they attenuate it. The square panel is the 2-D spectrum of the separable kernel, $|H(f_x)H(f_y)|$, out to one cycle per sample, with the source grid's Nyquist square in black and the turned grid's dashed. Only frequencies inside both pass as themselves; what the kernel lets through elsewhere, shown red, comes from the spectrum's copies or from corners the turned grid cannot hold, and folds back as aliasing. The readout gives its share of the kernel's energy: about 29% for nearest and 8% for Lanczos at a $10°$ step. Switched to $N$ passes, the panel multiplies the kernel turned by every step, and only what fits inside every grid along the way passes as itself, outlined dashed: the square rounds into a near disc of radius one half, and the readout gives both the share left outside after all the passes and the share a single pass aliases.

6.5.11 Prefiltering when the transform isn't a clean scale⧉

An affine map sends a square output footprint to a parallelogram in the input. Elliptical weighted average (EWA) filtering approximates the footprint with an ellipse. Its shape can be obtained from the Jacobian $J$ of the inverse map: a circular footprint with covariance $C$ maps locally to covariance $JCJ^\top$. An affine map has the same Jacobian everywhere; a perspective or nonlinear map generally does not. A pure rotation changes orientation without changing area or introducing unequal axis scales.

Mipmaps choose among prefiltered resolutions cheaply; EWA accounts more closely for directional compression. Both approximate the combined reconstruction and output filtering. Chapter 11.2, Antialiasing for complex transforms develops the footprint model.

Margin note: resample in linear light, not in gamma

A resample is a weighted average of neighboring pixels, and an average is only physically meaningful on values proportional to light, i.e. linear radiance. Most image files are gamma-encoded (sRGB), and averaging the encoded values is wrong: the encoding is concave, so the mean of two encoded values is darker than the encoding of their mean. The visible symptom is the classic "gamma-incorrect" darkening on downsampling and on antialiased edges: thin bright features on dark backgrounds lose energy, and mip-mapped textures go muddy. The fix is the standard ISP discipline (Chapter 45.15, Rendering a raw DNG): decode to linear, resample, re-encode, once, around the whole pass. Log space is a third option, but a perceptual / dynamic-range one rather than a physically correct one: averaging in log is a geometric mean that down-weights a few very bright pixels and can look better on extreme-HDR content, but for ordinary resampling linear is the correct default and gamma is the trap.

This matters on a receding surface. One output pixel may cover little source extent in one direction and a long strip in the other. Choosing one isotropic mipmap level large enough for the strip can erase useful detail across it; choosing a finer level can leave aliasing along it. An anisotropic footprint spends its averaging in the compressed direction. Ripmaps, multiple texture probes and EWA make different approximations to this same coverage problem, developed in Chapter 11.2, Antialiasing for complex transforms.

6.5.12 Resampling and encoding⧉

To approximate spatial integration of light, decode the transfer curve, filter in linear light, then encode again (Chapter 4.3, Encoding: linear, gamma or log). With nonnegative normalized weights, sRGB decoding is convex, so decoding the average code value gives no more light than averaging decoded values. Signed kernels such as Lanczos do not obey this inequality.

An equal mixture of black and white illustrates the size of the error: their sRGB codes average to 0.5, which decodes to about 0.214 linear light; their light average is 0.5, encoded as about 0.735. Low-contrast regions have smaller discrepancies. A deliberately chosen perceptual-space operation is possible, but it is a different operation from integrating light.

fig-resample-encoding-2d
Figure 6.5.17. The same downsample, twice: on the code values, and on the light (interactive). The large panel shows the original at 1:1 beside the reductions, enlarged, the one averaged in linear light in the top half and the one averaged in code values in the bottom half (a toggle shows either alone): drag the line between original and reduction, or switch the panel to either one whole, and the readout gives the mean light of the original and of each reduction as a code value. A pattern with fine high-contrast detail is shrunk by a factor you choose (2 to 16, 8 at first) with the filter you choose, once by averaging the sRGB code values directly and once by decoding to light, averaging, and re-encoding. The gamma-space result is darker, and it is darker exactly where the picture has detail: the zone plate develops a dark bruise that spreads outward as its rings tighten, while a flat region does not move at all. The reason is Jensen's inequality. sRGB's decoding is convex, so for weights that are nonnegative and sum to one, decoding the average of the codes lands at or below the average of the light, with equality where the neighborhood is flat or stays within the curve's linear toe. Change the filter and watch how little changes: box, bilinear, bicubic and Lanczos differ in sharpness and ringing, while the mean darkening moves only modestly, because most of the error comes from averaging codes instead of light, whatever the weights. The exact numbers depend on the pattern, the contrast and the filter. One caveat the figure keeps rather than hides: the pixel-by-pixel guarantee holds only for filters whose weights are all non-negative, box and bilinear here. Bicubic and Lanczos have negative lobes, so their output is not an average, and some pixels do come out brighter; the signed difference map shows those in blue. Raise the contrast and the gap grows, because contrast and frequency both widen the spread of values inside one output pixel. The patterns are generated in linear light on purpose: a zone plate authored straight into sRGB would already contain the mistake before any resampling happened.
fig-resample-encoding-1d
Figure 6.5.18. The same effect on one signal, where the mechanism is visible. A square-wave chirp in linear light is box-averaged into a quarter as many samples, two ways: averaging the light (correct) and averaging the sRGB code values and then decoding. The two traces are overlaid in light, and the gamma one sits systematically below, with the gap opening as the oscillation speeds up and closing to nothing on the one bin where the signal happens to be flat. That is Jensen's inequality exactly: sRGB's decoding is convex, so decoding the average of the codes is at most the average of the light, with equality only when there is nothing to average over. The bottom strip renders both as you would see them: at this contrast the gamma row is darker by up to 45 code values out of 255, a third of the way to black, on a difference nobody chose.

Do not infer the processing domain from a library or filter name. Set it explicitly and check with a black-and-white pattern. For transparent images, filter premultiplied linear RGB and alpha together, then unpremultiply where alpha is nonzero. Otherwise arbitrary colors in transparent pixels can bleed into visible edges.

6.5.13 Beyond Nyquist-limited reconstruction⧉

The fixed kernels above are linear for fixed geometry and weights, but separable kernels are generally not isotropic. They cannot determine arbitrary detail missing from the samples. Content-dependent reconstruction uses extra assumptions, such as locally smooth regions separated by edges.

But, exactly as with deblurring, we are rarely without further assumptions. Natural images are not arbitrary signals: they occupy a thin, structured manifold (the natural-image prior of Chapter 7.1, Imaging as Bayesian inference with probabilistic priors), and we can exploit those regularities to beat the worst case. One regularity matters here especially: images are locally anisotropic. Near an edge there is content in essentially one direction (across the edge) and almost none along it. That has a precise Fourier consequence: the local spectrum is concentrated on a line (or a narrow wedge), so when sampling replicates the spectrum, the alias copies land in only part of the Fourier plane, leaving large regions where we know no alias can be.

If the local edge orientation is known, an elongated filter can average along the edge while preserving variation across it. This can improve reconstruction under that scene model. It does not separate every overlap of spectral replicas: some edge orientations and textures remain ambiguous. Estimating orientation from the same aliased data introduces another possible error. The complete adaptive procedure is generally nonlinear even though filtering with frozen weights is linear.

This steered, content-adaptive strategy is the workhorse of demosaicking and of older super-resolution. Many demosaicking methods follow exactly this two-step recipe: first analyze the local orientation and pattern (is there an edge here, and which way does it run?), then steer a reconstruction filter to interpolate along the edge rather than across it, which is what keeps demosaicked edges from breaking into the colored zipper artifacts a naive isotropic interpolation produces.

The modern version of the same idea drops the hand-built orientation analysis and instead learns the regularities of natural images directly (and learns how to reconstruct them) from data. That is where super-resolution and demosaicking have largely gone (and where this book picks them up again in the learning chapters): the prior is no longer "images are locally oriented," written by hand, but a model trained to know what the high-resolution image probably was.

fig-sampling-steered-recon
Figure 6.5.19. Steering the reconstruction beyond Nyquist. The same 2-D sampling pipeline as the sampling chapter's figure, but now the reconstruction filter is an oriented, anisotropic kernel you aim by hand: set its orientation and its scale along each of its two axes to form an elliptical (sharp-bicubic) footprint. Line it up along an edge (wide along, narrow across) and the undersampled edge reconstructs with far less staircasing, because the filter averages where the signal is smooth and stays sharp where the detail is; misalign it and the aliasing returns. This relies on the edge being straight and on the direction being right, which you supply here and a practical method must estimate; it does not recover arbitrary aliased detail. It is the hand-driven version of what demosaicking and classical super-resolution do automatically: estimate the local orientation, then steer. Flip on the isotropic (round) filter toggle to see the standard, non-steered reconstruction for contrast: round blurs the staircase without removing it; steering along the edge does much better. The spectral panels below are schematic envelopes, not transforms of the displayed samples.

Finally, there is a second kind of luck, and it changes things entirely. Sometimes we are not reconstructing blind: we already hold a high-resolution image that we know is correlated with the very quantity we are trying to recover. This is the common situation in colorization and in local image adjustments, where the thing being upsampled (the chrominance, or a per-pixel correction) is often correlated with the high-resolution guide’s edges. Two families of methods turn that correlation into resolution. Cross / joint bilateral filtering reads its edge-stopping weights from the high-resolution guide while it filters the low-resolution quantity, so the upsampled result inherits the guide's edges (forward ref: Chapter 10.1, Bilateral filtering §cross/joint, and Chapter 10.6, Guided image filtering). Or we can pose an optimization that asks the output to differ between two neighboring pixels only as much as the input does, encouraging the correlations among output pixels to match those among input pixels (forward ref: Chapter 10.3, Edge-preserving optimization — colorization and Chapter 10.5, Matting Laplacian). Either way the high-resolution guide supplies candidate structure; mismatched edges can transfer texture or create false boundaries, a recurring theme of the chapters in Edges matter.

6.5.14 Practical implementation choices⧉

For a reproducible resize, specify the coordinate convention, kernel and scale, boundary extension, processing domain, and alpha treatment. Check that a constant image stays constant, an identity transform reproduces samples for an interpolating kernel, and a fine pattern is adequately filtered during reduction. Compose successive transforms before sampling.

Video 6.5.1. Video: Interpolation in 5 minutes, a short visual tour of nearest-neighbor, linear, and cubic interpolation, by Graphics in 5 Minutes (YouTube).
fig-xkcd-1416
Figure 6.5.20. xkcd #1416 "Pixels", Randall Munroe, xkcd.com/1416 (CC BY-NC 2.5). A recursive zoom into ever-finer pixels, an apt closing image for a chapter about reconstructing values between samples.

Recap: big lessons of this chapter

💡 Big lesson (L6.19) — resample by gathering, rarely by scattering

Loop over output pixels and pull back through $f^{-1}$ to find where each one reads from. Gathering fills every output exactly once; scattering pushes input pixels forward and leaves holes whose shape depends on the transform, a lattice of gaps under a scale and tilted slivers under a rotation. The same reflex governs warping, morphing, texture mapping, and flow-based reprojection later. Scattering earns its place only when the forward map is the thing you actually have, as in splatting a point cloud, and then it needs an explicit accumulation and normalization pass.

fig-resample-forward-inverse
Figure 6.5.2 (reproduced). Forward vs. inverse mapping on a 2× upscale of a small rainbow grid. Left (forward): loop over input pixels and push each to $(2i,2j)$; only one cell of every 2×2 output block gets written, leaving the other three as black gaps (a regular lattice of holes; worse under a rotation, where the holes become tilted slivers). Right (inverse): loop over output pixels and, for each, apply $f^{-1}$ to reach back into the input and sample, so every pixel is filled, no gaps. The arrows trace the two routes for a few pixels: forward arrows start at an input pixel and end where it is written, inverse arrows start at an output pixel and end at the input pixel it reads. The whole chapter loops over the output and reaches back into the input.
💡 Big lesson (L6.20) — every resampling pass multiplies the spectrum again, so resample once

Repeated resampling accumulates blur, ringing, and sampling error. For repeated fixed convolution the frequency response is exactly $H^N$; for rotations it is only an illustration because the grid and effective filter change at each pass. Compose geometric transforms and resample the original once.