Image Interpolation

Image interpolation estimates an image’s values between its pixels. Every operation that reads an image at a position that is not a pixel centre needs it: resizing, rotating, rectifying and registering images, correcting lens distortion, mapping textures, and sampling features at sub-pixel positions inside a network. The choice of method decides how much an image blurs, whether edges ring or turn into staircases, whether a flat region keeps its brightness, and how fast errors pile up when an image is resampled again and again. Image warping introduces nearest-neighbour, bilinear and bicubic interpolation as the step after the coordinate mapping, and sampling and aliasing treats interpolation as reconstruction from samples. This article is about the kernels themselves: how they are built, what their frequency responses say, and how they compare when measured on the same footing.

Definition

Given samples $v_{m,n}$ on the integer grid, interpolation constructs a function $u(x, y)$ defined at every real position that passes through the samples, $u(m, n) = v_{m,n}$, and approximates the continuous image from which they were taken (Getreuer, 2011, §1). A method that is linear, shift-invariant and local can always be written as a convolution of the samples with an interpolation kernel $K$ (Getreuer, 2011, §3):

$$u(x, y) = \sum_{m, n} v_{m,n}\, K(x - m,\, y - n).$$

The kernel is usually designed separable, $K(x, y) = K_1(x)\,K_1(y)$ (Getreuer, 2011, §3), so a 2D interpolation is a 1D interpolation along the rows followed by one along the columns, and everything below is stated for the 1D kernel $K_1$, written $K$. For $u$ to pass through the samples, $K$ must be interpolating: $K(0) = 1$ and $K(n) = 0$ at every other integer.

Intuition

Each sample spreads its value over its neighbourhood in the shape of the kernel, and the interpolated image is the sum of these spread-out copies. A box one pixel wide gives a piecewise-constant image of flat tiles: nearest-neighbour interpolation. A triangle two pixels wide joins neighbouring samples with straight lines: linear interpolation, and bilinear in 2D. Wider, smoother kernels with small negative lobes, such as cubic convolution or a windowed sinc, produce smoother images and sharper transitions, at the price of a slight dip and overshoot on either side of an edge.

The same picture has a frequency-domain reading. Sampling makes the spectrum of the image periodic, with copies of the original spectrum centred on every multiple of the sampling rate (sampling and aliasing). Interpolation multiplies that periodic spectrum by the kernel’s Fourier transform $\hat K(f)$. The ideal kernel keeps everything below the Nyquist frequency of half a cycle per pixel unchanged and removes everything above it; that kernel is the sinc. A practical kernel does two things wrong at once: it attenuates detail just below Nyquist, which blurs, and it lets through part of the spectral copies above Nyquist, which Mitchell and Netravali (1988) call postaliasing, as described in sampling and aliasing. Briand (2019, §1) describes the Fourier transform of a typical kernel as a smooth approximation of the indicator of the Nyquist band.

Formal Definition

Interpolating kernels

Nearest neighbour uses the box $K(t) = 1$ for $-\tfrac12 \le t < \tfrac12$ and $0$ otherwise; the half-open interval decides which sample wins at an exact tie (Getreuer, 2011, §7). Linear interpolation uses the triangle $K(t) = \max(0, 1 - |t|)$; its weights are non-negative and sum to one, so the result stays between the smallest and largest sample and never overshoots (Getreuer, 2011, §8). Their Fourier transforms are $\operatorname{sinc}(f)$ and $\operatorname{sinc}^2(f)$, with $\operatorname{sinc}(f) = \sin(\pi f)/(\pi f)$.

Cubic convolution (Keys, 1981) is the piecewise cubic that is symmetric, continuously differentiable and supported on $(-2, 2)$; these conditions leave one free parameter $a$ (Getreuer, 2011, §9; Thévenaz et al., 2000, §V-F):

$$K_a(t) = \begin{cases} (a + 2)|t|^3 - (a + 3)|t|^2 + 1, & |t| < 1,\\ a|t|^3 - 5a|t|^2 + 8a|t| - 4a, & 1 \le |t| < 2,\\ 0, & \text{otherwise}. \end{cases}$$

Getreuer notes that the values $-1$, $-0.75$ and $-0.5$ have been proposed for $a$, and that $a = -0.5$ makes the method third-order accurate while every other value gives first order; Thévenaz et al. likewise give approximation order 3 as the best a Keys kernel reaches, at $a = -1/2$. Getreuer calls interpolation with this kernel bicubic. Pillow 12.0.0 hard-codes $a = -0.5$ in its BICUBIC filter, and ImageMagick’s source lists it as Catmull–Rom, the member $B = 0$, $C = \tfrac12$ of the two-parameter family below.

Windowed sinc. The ideal kernel $\operatorname{sinc}(t)$ has infinite support and decays only like $1/t$. Truncating it with a window gives a compact kernel; the Lanczos window is itself a stretched sinc (Getreuer, 2011, §11):

$$L_n(t) = \operatorname{sinc}(t)\,\operatorname{sinc}(t/n), \qquad |t| < n,$$

with $n = 2$ or $3$ in practice. Lanczos kernels are interpolating, but their shifted copies do not sum exactly to one, so Getreuer recommends normalizing the weights at each position by their sum. Pillow’s LANCZOS filter is $L_3$, and its resampler divides every set of weights by their sum for all filters (Pillow 12.0.0, Resample.c).

Non-interpolating kernels

Mitchell–Netravali. Mitchell and Netravali (1988) studied a family of piecewise cubic filters to balance blurring against ringing, as sampling and aliasing describes. ImageMagick’s source writes the family with two parameters $B$ and $C$:

$$K_{B,C}(t) = \frac{1}{6} \begin{cases} (12 - 9B - 6C)|t|^3 + (-18 + 12B + 6C)|t|^2 + (6 - 2B), & |t| < 1,\\ (-B - 6C)|t|^3 + (6B + 30C)|t|^2 + (-12B - 48C)|t| + (8B + 24C), & 1 \le |t| < 2, \end{cases}$$

and labels $B = C = \tfrac13$ “Mitchell” (ImageMagick 7.1.2-35, resize.c). With $B > 0$ the kernel is not interpolating: $K(0) = 1 - B/3$ and $K(\pm1) = B/6$, so it smooths the samples it passes through. $B = 1$, $C = 0$ is the cubic B-spline used without a prefilter.

B-splines and generalized interpolation. The B-spline of degree $n$, $\beta^n$, is the box convolved with itself $n$ times; the cubic is

$$\beta^3(t) = \begin{cases} \tfrac23 - |t|^2 + \tfrac12|t|^3, & |t| < 1,\\ \tfrac16 (2 - |t|)^3, & 1 \le |t| < 2. \end{cases}$$

For $n \ge 2$ it is not interpolating, but it can still interpolate if the samples are first converted to coefficients $c_k$ with $\sum_k c_k \beta^n(m - k) = v_m$ at every sample (Thévenaz et al., 2000, §II; Getreuer, 2011, §4):

$$u(t) = \sum_k c_k\, \beta^n(t - k), \qquad c = (b^n)^{-1} * v,$$

where $b^n_m = \beta^n(m)$. Thévenaz, Blu and Unser call this generalized interpolation. The inverse filter $(b^n)^{-1}$ factors into pairs of first-order causal and anti-causal recursive filters, one pair per pole; for the cubic spline there is a single pole $r = \sqrt3 - 2$, the prefilter costs two additions and three multiplications per coefficient (Thévenaz et al., 2000, §II-D), and Briand and Monasse (2018) credit the recursive algorithm to Unser et al. (1991) and treat its boundary handling in detail. The two steps together are equivalent to a single interpolating kernel of infinite support, the cardinal spline; for the cubic,

$$\eta^3(t) = \sum_k \sqrt3\, r^{|k|}\, \beta^3(t - k), \qquad \hat\eta^3(f) = \frac{3\operatorname{sinc}^4(f)}{2 + \cos 2\pi f},$$

which oscillates like a sinc and converges to it as the degree grows (Getreuer, 2011, §12; Briand & Monasse, 2018, §2.1). The coefficients follow from Briand and Monasse’s exponential filters (their eq. 21, with $\gamma = 6$), and SciPy’s ni_splines.c stores the same pole.

Approximation order

Thévenaz et al. (2000, §IV) rank kernels by the rate at which interpolation error falls as the sampling step $h$ shrinks: a kernel of approximation order $L$ gives an error proportional to $h^L$ for smooth images. By the Strang–Fix conditions this is equivalent to reproducing every polynomial of degree $L - 1$ exactly; order 1 is the partition of unity, $\sum_n K(t - n) = 1$ for every $t$, which means a constant image stays constant. Nearest neighbour has order 1, linear order 2, Keys with $a = -1/2$ order 3, and the cubic B-spline order 4, which Thévenaz et al. contrast with Keys: both have support 4, but the spline has the higher order and is twice continuously differentiable where Keys is once. For sincs truncated by a rectangular (Dirichlet) or Hanning window, which do not reproduce constants, they find approximation order $L = 0$, so a finer sampling step need not reduce the error (§VI); the Lanczos window is not among the windows they analyze, but the example below shows that Lanczos kernels miss the partition of unity too. Getreuer (2011, §5) summarizes the same conditions and the error kernel that predicts interpolation error frequency by frequency.

Properties

• Interpolating is not the same as accurate. The central argument of Thévenaz et al. (2000) is that the kernels that perform best are not interpolating; with a prefilter they still pass through the samples. They attribute the belief that B-splines blur to implementations that skip the prefilter.
• Support sets the cost. A separable kernel of support $W$ needs $W^2$ weights per output pixel in 2D: 1 for nearest, 4 for bilinear, 16 for the cubics and Lanczos-2, 36 for Lanczos-3. The B-spline’s prefilter adds a fixed cost per pixel, not per output sample, and Briand and Monasse (2018, §6.1) measured it as negligible next to the interpolation step.
• Negative lobes ring. A kernel with negative values overshoots at edges; non-negative kernels (box, triangle, unprefiltered B-splines) cannot, but they blur. Getreuer (2011, Fig. 12) shows a step edge interpolated by six methods as a balance between staircases and ripples.
• Boundaries need a rule. Near the image border the sum reaches samples that do not exist. Getreuer (2011, §14.1) lists constant, half-sample symmetric and whole-sample symmetric extensions; Thévenaz et al. (2000, §II-C) advocate mirroring. For the B-spline the boundary rule also enters the prefilter, and SciPy’s documentation warns that results are only correct if the mode used for interpolation matches the one used for filtering.
• Grids are a convention. When scaling by a factor $d$, output samples can sit on a grid anchored at the top-left sample or on a centred grid; the centred grid makes interpolation commute with flipping the image (Getreuer, 2011, §14.2). Image warping lists the half-pixel offsets that follow from such conventions as a common bug.
• Shrinking needs a prefilter. Interpolation reconstructs; it does not remove frequencies that the smaller output grid cannot hold. When an image is reduced, the kernel must be widened by the reduction factor or the image smoothed first (sampling and aliasing, image pyramid). Pillow does the former: its resampler stretches the filter’s support by the scale factor whenever it downscales (Pillow 12.0.0, Resample.c).

Examples

Kernels compared on one footing

The script below computes, for each kernel, four properties from its formula and one from an image:

• frequency response $\hat K(f)$ at $f = 0.25$ and $0.5$ cycles per pixel, by Gauss–Legendre quadrature, checked against the closed forms for the box, the triangle and the cardinal cubic spline;
• energy above Nyquist, the fraction of $\int K^2 = \int |\hat K|^2$ at $|f| > \tfrac12$, the part of the kernel that passes spectral copies (postaliasing);
• partition of unity, $\max_t |\sum_n K(t - n) - 1|$, together with $|K(0) - 1|$ and $K(1)$, which show whether the kernel is interpolating;
• step overshoot, the maximum of the interpolated unit step (samples 0 to the left, 1 from the origin on) minus 1;
• repeated half-pixel shifts: a grey crop of scikit-image’s bundled rocket photograph (rows 160–415, columns 256–511), mirrored along $x$ so that its rows are periodic with an odd width of 511, is shifted right by half a pixel $N = 2, 8, 32, 128$ times. After $N$ shifts it should equal the input moved by exactly $N/2$ pixels, which np.roll gives without error, so the RMS difference is the error the kernel has accumulated. The script also prints the change of the image mean.

Lanczos kernels are printed raw and with their weights normalized, because implementations differ (Getreuer recommends normalizing; Pillow does it). The cubic B-spline is used with its prefilter, as the cardinal kernel $\eta^3$ truncated at $|k| \le 40$. An ideal shift by a phase factor in the discrete Fourier transform is the reference; the odd width means there is no Nyquist bin whose half-pixel shift would be ambiguous. Conventions: arrays are indexed (row, column), shifts act along rows, one half-pixel shift computes $u(x - \tfrac12)$ on the circular grid, and frequencies are in cycles per pixel. Three checks compare the hand-written shifts with their prediction in the Fourier domain and with scipy.ndimage.shift for the kernels SciPy offers (orders 0, 1 and 3, with mode="grid-wrap").

Before the script was first run we wrote down six hypotheses from the sources, and computed the values a correct implementation reaches with an independent script (closed-form spectra, dense grids and Fourier-domain prediction of the shifts); the example therefore illustrates the analysis rather than testing it. The hypotheses: (H1) at $f = 0.25$ the cubic spline passes at least as much as Keys, Keys more than linear, Mitchell–Netravali less than Keys, and nearest more than linear; (H2) nearest has more energy above Nyquist than linear, linear more than every cubic and Lanczos kernel, Lanczos-3 less than Lanczos-2 and the cubic spline less than Keys; (H3) nearest and linear do not overshoot, every other kernel does, Lanczos-3 more than Lanczos-2, and Mitchell–Netravali least; (H4) the partition of unity holds to $10^{-12}$ except for raw Lanczos, and holds for normalized Lanczos; (H5) the shift error grows with $N$ for every kernel but the ideal shift, the cubic spline stays below Keys and Keys below linear at every $N$, Mitchell–Netravali stays above Keys, and nearest neighbour does not move the image at all, because at a tie of exactly half a pixel the half-open box picks the same sample; (H6) raw Lanczos changes the image mean and no other kernel does. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 on one BLAS thread, runs in about 7 seconds, and two runs printed identical output. After a pilot run and before the first run, the editor approved three changes: the partition-of-unity cell of normalized Lanczos, which the pilot printed as a constant instead of computing it, is now computed; $K(0)$, which the pilot printed to four decimals, is printed as $|K(0) - 1|$ in scientific notation so that the $10^{-12}$ check can be read, with $K(0)$ and $K(1)$ of Mitchell–Netravali on their own line; and a line printing the change made by one nearest-neighbour half-pixel shift was added. No other printed value changed.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                      # one BLAS thread, set before NumPy is imported
import numpy as np
from scipy import ndimage
from skimage import data, color

# ---- Kernels (1-D; images are interpolated separably or, here, along rows) ----
def nearest(t):
    t = np.asarray(t, float)
    return ((t >= -0.5) & (t < 0.5)).astype(float)      # box on [-1/2, 1/2)

def linear(t):
    return np.clip(1 - np.abs(t), 0, None)

def keys(t, a=-0.5):                                    # Keys cubic convolution
    x = np.abs(t)
    return np.where(x < 1, (a + 2) * x**3 - (a + 3) * x**2 + 1,
                    np.where(x < 2, a * x**3 - 5 * a * x**2 + 8 * a * x - 4 * a, 0.0))

def mitchell(t, B=1/3, C=1/3):                          # Mitchell-Netravali (B, C) cubic
    x = np.abs(t)
    p = ((12 - 9*B - 6*C) * x**3 + (-18 + 12*B + 6*C) * x**2 + (6 - 2*B)) / 6
    q = ((-B - 6*C) * x**3 + (6*B + 30*C) * x**2 + (-12*B - 48*C) * x + (8*B + 24*C)) / 6
    return np.where(x < 1, p, np.where(x < 2, q, 0.0))

def lanczos(n):
    return lambda t: np.where(np.abs(t) < n, np.sinc(t) * np.sinc(np.asarray(t) / n), 0.0)

def bspline3(t):
    x = np.abs(t)
    return np.where(x < 1, 2/3 - x**2 + x**3 / 2, np.where(x < 2, (2 - x)**3 / 6, 0.0))

R3 = np.sqrt(3) - 2                                     # pole of the cubic B-spline prefilter

def cubic_spline(t, M=40):
    """Cardinal cubic spline: B-spline interpolation with its prefilter folded into one kernel,
    sum_k sqrt(3) r^|k| beta3(t - k), truncated at |k| <= M (r^40 is about 1e-23)."""
    t = np.asarray(t, float)
    return sum(np.sqrt(3) * R3**abs(k) * bspline3(t - k) for k in range(-M, M + 1))

KERNELS = {"nearest": (nearest, 0.5), "linear": (linear, 1), "Keys a=-0.5": (keys, 2),
           "Mitchell-Netravali": (mitchell, 2), "Lanczos-2": (lanczos(2), 2),
           "Lanczos-3": (lanczos(3), 3), "cubic B-spline": (cubic_spline, 42)}

# ---- 1. Frequency response by piecewise Gauss-Legendre quadrature ----
xg, wg = np.polynomial.legendre.leggauss(16)

def nodes(a, b, piece):
    edges = np.arange(a, b + piece / 2, piece)
    lo, hi = edges[:-1, None], edges[1:, None]
    return ((hi - lo) / 2 * xg + (hi + lo) / 2).ravel(), ((hi - lo) / 2 * wg).ravel()

def response(K, R, f):                                  # K even: H(f) = 2 int_0^R K cos
    t, w = nodes(0, R, 0.5)
    return 2 * (w * K(t)) @ np.cos(2 * np.pi * np.outer(t, f))

def energy_above_nyquist(K, R):
    t, w = nodes(0, R, 0.5)
    total = 2 * np.sum(w * K(t)**2)                     # Parseval: int K^2 = int |H|^2
    fq, wf = nodes(0, 0.5, 0.05)
    inband = 2 * np.sum(wf * response(K, R, fq)**2)
    return (total - inband) / total

analytic = {"nearest": lambda f: np.sinc(f), "linear": lambda f: np.sinc(f)**2,
            "cubic B-spline": lambda f: 3 * np.sinc(f)**4 / (2 + np.cos(2 * np.pi * f))}
dev = max(np.abs(response(*KERNELS[k], np.array([0.25, 0.5])) - h(np.array([0.25, 0.5]))).max()
          for k, h in analytic.items())
print(f"check: quadrature vs closed-form spectra (nearest, linear, cubic B-spline): max dev {dev:.1e}")

print(f"\n{'kernel':19s} {'H(0.25)':>8s} {'H(0.5)':>7s} {'E>Nyq':>7s} {'sum-1':>8s} "
      f"{'|K(0)-1|':>8s} {'K(1)':>8s} {'overshoot':>9s}")
tt = np.arange(0, 1, 1e-3)
ts = np.arange(-4, 4, 1e-4)
for name, (K, R) in KERNELS.items():
    H = response(K, R, np.array([0.25, 0.5]))
    rng = range(-int(R) - 2, int(R) + 3)
    wsum = sum(K(tt - n) for n in rng)
    pou = np.abs(wsum - 1).max()                          # partition of unity
    step = sum(K(ts - n) for n in range(0, int(R) + 8))   # samples 0 for n < 0, 1 for n >= 0
    print(f"{name:19s} {H[0]:8.4f} {H[1]:7.4f} {energy_above_nyquist(K, R):7.4f} {pou:8.1e} "
          f"{abs(float(K(0.0)) - 1):8.1e} {float(K(1.0)):8.1e} {step.max() - 1:9.4f}")
    if name.startswith("Lanczos"):                       # normalized: weights divided by their sum
        den = sum(K(ts - n) for n in range(-int(R) - 8, int(R) + 8))
        pou_n = np.abs(sum(K(tt - n) / wsum for n in rng) - 1).max()
        print(f"{'  normalized':19s} {'':8s} {'':7s} {'':7s} {pou_n:8.1e} {'':8s} {'':8s} "
              f"{(step / den).max() - 1:9.4f}")

print(f"Mitchell-Netravali: K(0) = {float(mitchell(0.0)):.4f}, K(1) = {float(mitchell(1.0)):.4f}")

# ---- 2. Repeated half-pixel shifts of a periodic image ----
g = color.rgb2gray(data.rocket())[160:416, 256:512]     # 256 x 256 crop: rocket, tower, lights
img = np.hstack([g, g[:, :0:-1]])                       # mirrored: periodic rows of odd width 511
W = img.shape[1]
print(f"\ndata: image {img.shape}, values {img.min():.3f}..{img.max():.3f}, "
      f"wrap-around step {np.abs(img[:, 0] - img[:, -1]).max():.4f}")

def half_shift_taps(K, R, normalize=False):
    j = np.arange(-int(R) - 2, int(R) + 3)
    h = K(j - 0.5)                                      # out[x] = sum_j h_j img[x - j] = u(x - 1/2)
    return j, (h / h.sum() if normalize else h)

def shift_spatial(a, j, h):
    return sum(hj * np.roll(a, jj, axis=1) for jj, hj in zip(j, h) if hj != 0)

j0, h0 = half_shift_taps(*KERNELS["nearest"])
print(f"nearest, one half-pixel shift: max |output - input| = {np.abs(shift_spatial(img, j0, h0) - img).max():.1e}")

xi = np.fft.fftfreq(W)
F = np.fft.fft(img, axis=1)
Ns = [2, 8, 32, 128]
methods = [(n, *KERNELS[n], False) for n in KERNELS] + \
          [("Lanczos-2 normalized", *KERNELS["Lanczos-2"], True),
           ("Lanczos-3 normalized", *KERNELS["Lanczos-3"], True)]
print(f"{'RMS error after N half-pixel shifts':37s}" + "".join(f"{'N=' + str(N):>9s}" for N in Ns)
      + f"{'mean drift N=128':>18s}")
worst_fft = 0.0
for name, K, R, norm in methods:
    j, h = half_shift_taps(K, R, norm)
    G = (h[:, None] * np.exp(-2j * np.pi * j[:, None] * xi[None, :])).sum(0)  # per-shift response
    a, row = img.copy(), []
    for n in range(1, Ns[-1] + 1):
        a = shift_spatial(a, j, h)
        if n in Ns:
            truth = np.roll(img, n // 2, axis=1)        # exactly shifted by n/2 pixels
            row.append(np.sqrt(np.mean((a - truth)**2)))
            pred = np.real(np.fft.ifft(F * G**n, axis=1))
            worst_fft = max(worst_fft, np.abs(a - pred).max())
    print(f"{name:37s}" + "".join(f"{e:9.5f}" for e in row) + f"{a.mean() - img.mean():+18.2e}")

a, row = img.astype(complex), []                        # ideal: DFT phase shift (no Nyquist bin)
for n in range(1, Ns[-1] + 1):
    a = np.fft.ifft(np.fft.fft(a, axis=1) * np.exp(-1j * np.pi * xi), axis=1)
    if n in Ns:
        row.append(np.sqrt(np.mean((a.real - np.roll(img, n // 2, axis=1))**2)))
print(f"{'ideal (DFT phase shift)':37s}" + "".join(f"{e:9.1e}" for e in row)
      + f"{a.real.mean() - img.mean():+18.2e}")
print(f"\ncheck: spatial shifts vs DFT prediction, max abs diff {worst_fft:.1e}")

for order, name in [(0, "nearest"), (1, "linear"), (3, "cubic B-spline")]:
    j, h = half_shift_taps(*KERNELS[name])
    a, b = img.copy(), img.copy()
    for n in range(Ns[-1]):
        a = shift_spatial(a, j, h)
        b = ndimage.shift(b, (0, 0.5), order=order, mode="grid-wrap")
    print(f"check: {name} vs scipy.ndimage.shift(order={order}, mode='grid-wrap'), "
          f"128 shifts: max abs diff {np.abs(a - b).max():.1e}")

Output:

check: quadrature vs closed-form spectra (nearest, linear, cubic B-spline): max dev 1.0e-15

kernel               H(0.25)  H(0.5)   E>Nyq    sum-1 |K(0)-1|     K(1) overshoot
nearest               0.9003  0.6366  0.2263  0.0e+00  0.0e+00  0.0e+00    0.0000
linear                0.8106  0.4053  0.0502  0.0e+00  0.0e+00  0.0e+00    0.0000
Keys a=-0.5           0.9390  0.4928  0.0470  2.2e-15  0.0e+00  0.0e+00    0.0741
Mitchell-Netravali    0.8450  0.3833  0.0322  2.2e-15  1.1e-01  5.6e-02    0.0363
Lanczos-2             0.9499  0.5006  0.0459  1.9e-02  0.0e+00  2.5e-17    0.1001
  normalized                                  4.4e-16                      0.0854
Lanczos-3             1.0084  0.5002  0.0287  5.7e-03  0.0e+00  3.2e-17    0.1121
  normalized                                  4.4e-16                      0.1178
cubic B-spline        0.9855  0.4928  0.0266  1.3e-15  0.0e+00 -9.0e-17    0.1078
Mitchell-Netravali: K(0) = 0.8889, K(1) = 0.0556

data: image (256, 511), values 0.058..1.000, wrap-around step 0.0157
nearest, one half-pixel shift: max |output - input| = 0.0e+00
RMS error after N half-pixel shifts        N=2      N=8     N=32    N=128  mean drift N=128
nearest                                0.05471  0.08421  0.15011  0.17582         +0.00e+00
linear                                 0.02055  0.03212  0.04128  0.05203         +0.00e+00
Keys a=-0.5                            0.01520  0.02339  0.03077  0.03597         +0.00e+00
Mitchell-Netravali                     0.01747  0.02784  0.03691  0.04542         -8.66e-15
Lanczos-2                              0.01952  0.05824  0.27476  3.33515         +3.07e+00
Lanczos-3                              0.01189  0.02199  0.06126  0.27485         -1.59e-01
cubic B-spline                         0.01125  0.01710  0.02316  0.02922         -1.52e-14
Lanczos-2 normalized                   0.01520  0.02339  0.03077  0.03597         -1.05e-15
Lanczos-3 normalized                   0.01125  0.01649  0.03276  0.46948         -2.22e-16
ideal (DFT phase shift)                1.8e-16  4.1e-16  1.0e-15  3.0e-15         -2.00e-15

check: spatial shifts vs DFT prediction, max abs diff 1.2e-13
check: nearest vs scipy.ndimage.shift(order=0, mode='grid-wrap'), 128 shifts: max abs diff 0.0e+00
check: linear vs scipy.ndimage.shift(order=1, mode='grid-wrap'), 128 shifts: max abs diff 0.0e+00
check: cubic B-spline vs scipy.ndimage.shift(order=3, mode='grid-wrap'), 128 shifts: max abs diff 8.2e-14

The checks pass: the quadrature matches the closed-form spectra, the hand-written shifts match their Fourier-domain prediction, and SciPy’s nearest, linear and cubic-spline shifts agree with ours over all 128 steps. SciPy’s order=3 is exactly the prefiltered cubic B-spline, and its order=0 breaks the half-pixel tie the same way as the half-open box. The ideal shift returns the exactly moved image to rounding error at every $N$.

Frequency response (H1, H2). H1 and H2 held: every planned ordering appears in the table. At a quarter cycle per pixel, halfway to Nyquist, linear interpolation keeps 81% of the amplitude, Mitchell–Netravali 85%, Keys 94% and the cubic spline 99%: this is the blur that each kernel adds to detail it could represent. Nearest neighbour keeps more of that detail than linear, but it also leaks the most above Nyquist, 23% of its energy against 5% for linear, because the box’s sinc spectrum decays slowly. The kernels of support 4 and more all leak less than linear, with Lanczos-3 and the cubic spline lowest. Lanczos-3 is the one kernel whose response exceeds 1 at $f = 0.25$: windowing a sinc leaves ripples in the passband, so it slightly amplifies some frequencies.

Interpolation and the partition of unity (H4). H4 held. The table confirms which kernels pass through the samples: all but Mitchell–Netravali, whose $K(0) = 8/9$ and $K(\pm1) = 1/18$ make it a smoothing filter. Nearest, linear, Keys, Mitchell–Netravali and the cubic spline reproduce constants to rounding error. Raw Lanczos-2 and Lanczos-3 do not: their shifted copies sum to values that deviate from 1 by up to 1.9% and 0.57%, the defect Getreuer’s normalization repairs, and the normalized rows meet the partition of unity to rounding error.

Ringing (H3). H3 held. Nearest and linear do not overshoot, as their non-negative weights guarantee; every kernel with negative lobes does. Mitchell–Netravali overshoots least, about 3.6% of the step, Keys about 7.4%, and the cubic spline about 10.8%. Lanczos-3 rings more than Lanczos-2, raw or normalized; its extra lobes bring it closer to the sinc, whose truncation is the source of the ringing.

Repeated shifts (H5, H6). H5 held: every kernel’s error grows with $N$, and the orderings held at every $N$: the cubic spline below Keys, Keys below linear, Mitchell–Netravali above Keys. At a half-pixel shift every output value lies midway between samples, and the shift’s frequency response is applied once per shift, so in our reading the attenuation compounds: after 128 shifts linear interpolation reaches an RMS error of 0.052 on a $[0, 1]$ scale, against 0.029 for the cubic spline. Nearest neighbour never moves the image, as the one-shift line shows, so after $N$ shifts the image is simply $N/2$ pixels behind where it should be. Mitchell–Netravali loses to Keys although it rings less; in our reading, because it is not interpolating, it smooths even where Keys reproduces the samples. Raw Lanczos fails differently: its weights at the half-pixel position do not sum to one, so each shift scales the mean, which after 128 shifts has risen by 3.07 for Lanczos-2 and fallen by 0.16 for Lanczos-3, and their errors grow far beyond those of the other kernels. Thévenaz et al. (2000, §VIII-D) report the same drift of the mean for apodized sincs in their repeated-rotation experiment. H6 held: normalizing removes the drift, and no other kernel changes the mean beyond rounding error.

Two results were not among the hypotheses, although the independent script had already shown them before the pilot. First, normalized Lanczos-2 gives exactly the same shifts as Keys. At the half-pixel position the weights are $L_2(\tfrac12) = 4\sqrt2/\pi^2$ and $L_2(\tfrac32) = -4\sqrt2/(9\pi^2)$, in the ratio $-9$, and Keys’ weights $K(\tfrac12) = 9/16$ and $K(\tfrac32) = -1/16$ are in the same ratio, so after normalization the two filters coincide (our derivation; it holds at this position only, not for the kernels in general, as their different step overshoots in the first table show). Second, normalized Lanczos-3 equals the cubic spline at $N = 2$ to the printed precision and is below it at $N = 8$, but its error jumps to 0.47 at $N = 128$, worse than linear. In our reading this follows from the passband ripple that the first table shows for the raw kernel: a filter that amplifies some frequencies by a little more than 1 amplifies them geometrically when applied 128 times. A kernel that slightly sharpens is harmless in one resampling and harmful in a long chain.

Thévenaz et al. (2000, §VIII) ran the 2D version of this experiment, 15 rotations of 24° each in floating point, and compared the result with the original. On their CT image (Table IV) the signal-to-noise ratio was 19.44 dB for nearest neighbour, 24.21 dB for linear, 38.74 dB for Keys and 50.77 dB for the cubic spline, the same ranking as here; the cubic o-Moms kernel, which has the cubic spline’s support and order and a smaller error constant, reached 62.05 dB. Briand (2019) turned this idea into a measure of interpolation quality, the reversibility error: the error after a random homography close to the identity, a crop and the inverse homography, averaged over many homographies.

Common Misconceptions

• “B-spline interpolation blurs.” Only if the B-spline is used directly as a kernel ($B = 1$, $C = 0$ in the Mitchell–Netravali family), which smooths the samples. With its prefilter it passes through every sample and has the sharpest passband of the cubic kernels compared above. Thévenaz et al. (2000, §V-B) single out this mistake.
• “Bicubic means one method.” Bicubic in Getreuer’s usage and in Pillow is Keys’ cubic convolution, for which values of $a$ from $-1$ to $-0.5$ have been proposed (Getreuer, 2011, §9) and Pillow uses $-0.5$; cubic interpolation in SciPy is the prefiltered cubic B-spline, whose order=3 is the default of map_coordinates, shift and zoom. The results differ, so a pipeline that must match another needs the same kernel, not just the same name (model deployment gives library defaults and the effect of a mismatch).
• “The sinc is the best interpolator.” It is exact only for band-limited images, and images with edges are not band-limited; Getreuer (2011, §10) shows sinc interpolation ringing around edges, and on his piecewise-smooth test image bicubic and Lanczos-2 had the lowest error of the methods compared (§16.2), while on smooth data the higher-order methods and sinc did best (§16.1).
• “Windowing a sinc keeps it exact on constants.” It does not: raw Lanczos weights do not sum to one, which shows as a drift of brightness over repeated resampling unless the weights are normalized.
• “A higher-order kernel always helps.” Order describes behaviour on smooth content. On edges, kernels with larger negative lobes ring more, as the step overshoots above show. Briand and Monasse (2018, §6.4), who recommend B-splines of higher order than cubic, measured the gain against Shannon–Whittaker interpolation of a photograph resampled by a homography and in a consistency test of ten shifts of 0.1 pixel followed by one of −1 pixel, an experiment close to the one above.
• “Interpolation adds detail.” It estimates values between samples under a smoothness assumption and adds no information; recovering detail beyond the grid is super-resolution, whose interpolation baselines are these kernels.

Where It Is Used

• Resizing. Enlarging an image is pure interpolation; shrinking it combines interpolation with a prefilter whose width follows the scale factor (sampling and aliasing, image pyramid). Dataset preparation depends on these choices: FFHQ was reduced with Pillow’s Lanczos filter, and the resizing inside the Fréchet Inception Distance changes its scores.
• Geometric transformations. Rotation, affine and projective warps, rectification and lens distortion correction evaluate the input at non-integer positions by backward mapping (image warping), and Getreuer (2011, §14.3) adds texture mapping and parallax correction.
• Registration and motion. Image registration and optical flow warp one image towards another at every iteration and differentiate the interpolated image; Thévenaz et al. (2000, §IX) prefer B-splines to o-Moms for such tasks because B-splines are smoother.
• Medical volumes. Resampling scans to isotropic voxels is interpolation along the slice axis, where thick slices limit what any kernel can recover (medical image preprocessing); Thévenaz et al. (2000, §I) motivate their analysis with such medical uses, where interpolation artifacts can affect a clinician’s reading.
• Networks. Bilinear upsampling by an integer factor is a transposed convolution with a fixed kernel, which can instead be learned; most super-resolution benchmarks synthesize their low-resolution inputs by bicubic downscaling, and interpolation is the baseline in their tables.