@article{ipol.2011.g_lmii, title = {{Linear Methods for Image Interpolation}}, author = {Getreuer, Pascal}, journal = {{Image Processing On Line}}, volume = {1}, year = {2011}, doi = {10.5201/ipol.2011.g_lmii}, } % if your bibliography style doesn't support doi fields: note = {\url{http://dx.doi.org/10.5201/ipol.2011.g_lmii}}
 published
 20110927
 reference
 Pascal Getreuer, Linear Methods for Image Interpolation, Image Processing On Line, 1 (2011). http://dx.doi.org/10.5201/ipol.2011.g_lmii
Communicated by Gabriele Facciolo
Demo edited by Pascal Getreuer
 Pascal Getreuer
pascal.getreuer@cmla.enscachan.fr
, CMLA, ENS Cachan
Overview
Given input image v with uniformlysampled pixels v_{m,n}, the goal of interpolation is to find a function u(x,y) satisfying
such that u approximates the underlying function from which v was sampled. Another way to interpret this is v was created by subsampling, and interpolation attempts to invert this process.
We discuss linear methods for interpolation, including nearest neighbor, bilinear, bicubic, splines, and sinc interpolation. We focus on separable interpolation, so most of what is said applies to onedimensional interpolation as well as Ndimensional separable interpolation.
References
 C.E. Shannon, “A Mathematical Theory of Communication,“ Bell System Technical Journal, vol. 27, pp. 379–423, 623–656, 1948.
 G. Strang and G. Fix, “A Fourier analysis of the finite element variational method,” Constructive Aspects of Functional Analysis, Edizioni Cremonese, Rome, pp. 795–840, 1973.
 R. Keys, “Cubic Convolution Interpolation for Digital Image Processing,” IEEE Trans. Acoustics, Speech, and Signal Processing, 29(6), 1981.
 A. Aldroubi, M. Unser, and M. Eden, “Cardinal spline filters: stability and convergence to the ideal sinc interpolator,” Signal Processing, vol. 28, no. 2, pp. 127–138, 1992.
 A. Schaum, “Theory and design of local interpolators,” CVGIP: Graph. Models Image Processing, vol. 55, no. 6, pp. 464–481, 1993.
 M. Unser, “Splines: A Perfect Fit for Signal/Image Processing,” IEEE Signal Processing Magazine, 16(6), pp. 22–38, 1999.
 P. Thévenaz, T. Blu, and M. Unser, “Interpolation Revisited,” IEEE Trans. Medical Imaging, vol. 19, no. 7, pp. 739–758, 2000.
 T. Blu, P. Thévenaz, and M. Unser, “MOMS: MaximalOrder Interpolation of Minimal Support,” IEEE Trans. Image Processing, vol. 10, no. 7, pp. 1069–1080, 2001.
 E. Meijering, “A Chronology of Interpolation: From Ancient Astronomy to Modern Signal and Image Processing,” In Proceedings of the IEEE, 90, pp. 319–342, 2002.
 F. Malgouyres and F. Guichard, “Edge direction preserving image zooming: A mathematical and numerical analysis,” SIAM J. Numerical Analysis, 39(1), pp. 1–37, 2002.
 B. Fornberg, E. Larsson, N. Flyer, “Stable computations with Gaussian radial basis functions,” SIAM J. Sci. Comput. 33, pp. 869–892, 2011.
 FFmpeg, http://www.ffmpeg.org.
Online Demo
An online demo of this algorithm is available.
Notations
We use the following normalization of the Fourier transform,
the transform of a function f is denoted by hat. We will also use the bilateral Ztransform. The Ztransform of a sequence c_{n} is
We denote the sequence with a subscript c_{n} and its Ztransform as c(z).
Interpolation Kernels
Linear interpolation can be abstractly described as a linear operator Z mapping samples v to a function u := Z(v). We define two properties on Z.
Let S_{k,l} denote shifting, S_{k,l}(v)(m,n) := v_{m−k,n−l}. Z is said to be shiftinvariant if for every v and (k, l),
Let v^{N} denote the restriction of v to the set {−N,…,+N}^{2}, then Z is said to be local if for every v and (x, y),
If Z is a linear, shiftinvariant, local operator from to , then there exists a unique function K ∈ L^{1} called the interpolation kernel such that Z(v) can be written as a convolution [10]:
The conditions for this result are mild and are satisfied by most linear interpolation methods of practical interest.
The properties of the kernel K, like symmetry, regularity, and moments, imply properties of the interpolation. For fourfold rotational symmetry, K is usually designed as a tensor product K(x,y) = K_{1}(x)K_{1}(y), where K_{1} is a symmetric function. Such a kernel also has the benefit that the interpolation is separable, the interpolation can be decomposed into onedimensional interpolations along each dimension.
Similarly in three (or higher) dimensions, it is convenient to use a separable kernel K(x,y,z) = K_{1}(x)K_{1}(y)K_{1}(z) so that the interpolation decomposes into onedimensional interpolations. More generally, one can use different onedimensional kernels along different dimensions, K(x,y,z) = K_{1}(x)K_{2}(y)K_{3}(z). For instance, it may be more appropriate to use a different kernel for the temporal dimension than for the spatial dimensions in video interpolation.
In order to have the interpolant agree with the known samples Z(v)(m,n) = v_{m,n}, it is easy to see that the requirements are K(m,n) = 0 for all integer m and n, except at the origin where K(0,0) = 1. K is said to be an interpolating function if it satisfies these requirements.
TwoStep Interpolation
That K must be interpolating can be an inconvenient restriction. It complicates the problem of finding a kernel K that also satisfies other design objectives. It can also be an obstacle computationally. In some cases such as spline interpolation, K has infinite support, yet the interpolant can be expressed as a linear combination of compact support functions (the Bsplines).
Suppose that one dimensional data f_{m} is to be interpolated. An alternative approach is to begin with a basis function φ(t) and then express the interpolant as a linear combination of {φ(t − n)},
where the coefficients c_{n} are selected such that
That is, interpolation is a twostep procedure: first, we solve for the coefficients c_{n} and second, the interpolant is constructed as ∑c_{n}φ(t − n).
To solve for the coefficients, notice that f_{m} = ∑c_{n}φ(m − n) is equivalent to a discrete convolution f_{m} = (c ∗ p)_{m} where p_{m} = φ(m). This suggests solving for c_{n} under the Ztransform as
Let (p)^{−1} deonte the convolution inverse of p, which is the sequence such that p ∗ (p)^{−1} = δ where δ is the unit impulse (δ_{0} := 1 and δ_{m} := 0 otherwise). Using (p)^{−1} we obtain c_{n} by “prefiltering” as
A unique convolution inverse exists in many practical cases of interest [7]. If it exists and φ is real and symmetric, then (p)^{−1} can be factored into pairs of recursive (infinite impulse response) filters, allowing for an efficient inplace calculation. Finally, the interpolant is obtained as ∑c_{n}φ(t − n).
It is possible to express this twostep interpolation procedure in terms of an interpolation kernel,
For example, the following figures illustrate this relationship for cubic and septic (7th) Bspline interpolation where φ are the cubic and septic Bsplines. In both cases, φ is not an interpolating function, but prefiltering with (p)^{−1} yields a kernel K_{1} that is interpolating.
Twostep interpolation where φ is the cubic Bspline. 
Twostep interpolation where φ is the septic Bspline. 
Interpolation Properties
An interpolation method Z has approximation order J if it reproduces polynomials up to degree (J − 1).
Let Z and K_{1} be defined from a function φ as in the previous section. The Strang–Fix conditions [2] provide the following equivalent conditions to determine the approximation order.
 Z has approximation order J,
where denotes the j_th derivative of and . Aside from approximation order, another property of interest is the L^{2} interpolation error. Let and define f_{}δ(t) = ∑f(δn)K_{1}(t − n) for δ_ > 0. The interpolation error is approximately
This approximation is exact for bandlimited functions. Otherwise, it is equal to the average L^{2} error over all translations of f. The function E_{int} is called the interpolation error kernel, and it can be interpreted as predicting the error made in interpolating sin(2πξt). Asymptotically as δ → 0,
The decay of the error is dominated by the approximation order J. When comparing methods of the same approximation order, the constant C_{int} is also of interest [7].
Kernel Normalization
A basic and desirable quality of an interpolation method is that it reproduces constants, which is equivalent to
If the kernel does not reproduce constants, it can be the normalized to fix this as
provided that the denominator does not vanish. Normalization ensures that the interpolation weights sum to 1. Interpolation with the normalized kernel reproduces constants: suppose f_{n} = c, then
Nearest Neighbor
The most widely used methods for image interpolation are nearest neighbor, bilinear, and bicubic interpolation.
Nearest Neighbor  Bilinear  Bicubic 
Nearest neighbor, bilinear, and bicubic applied to the same uniformlyspaced input data. 
The nearest neighbor interpolation of v is the piecewise constant function
where [⋅] denotes rounding to the nearest integer. That is, u(x,y) is defined as the value of the input sample that is closest to (x,y). For this reason, nearest neighbor interpolation is sometimes called “pixel duplication.” The interpolation kernel for nearest neighbor is
The nearest neighbor idea can be applied to very general kinds of data, to any set of samples that has a notion of distance between sample locations. The interpolant is constant within the Voronoi cell around each sample location.
Nearest neighbor scattered data interpolation. 
Bilinear
The bilinear interpolation of v is the continuous function
where ⌊⋅⌋ denotes the floor function and is the fractional part. The interpolation kernel is
where (⋅)^{+} denotes positive part.
Within each cell [m,m+1]×[n,n+1], the interpolation is a convex combination of the samples located at the cell corners v_{m,n}, v_{m+1,n}, v_{m,n+1}, v_{m+1,n+1}. Because this combination is convex, the interpolation is bounded between min v and max v and does not produce overshoot artifacts. Bilinear interpolation reproduces affine functions: if v_{m,n} = am + bn + c then u(x,y) = ax + by + c.
Bilinear interpolation is arguably the simplest possible separable method that produces a continuous function. It is extremely efficient and on many platforms available in hardware, making it practical for realtime applications.
A closely related method to bilinear interpolation is linear interpolation. Given a triangulation of the sample locations, each triangle is interpolated by the affine function that satisfies the data at the triangle vertices.
Bicubic
Cubic interpolation kernel. 
Surface plot of K_{1}(x)K_{1}(y). Yellow circles indicate local maxima. 
Bicubic interpolation [3] uses the interpolation kernel
where α is a free parameter. This function is derived by finding a piecewise cubic polynomial with knots at the integers that is required to be symmetric, C^{1} continuous, and have support in −2<t<2. These conditions leave one remaining degree of freedom represented by α. For any value of α, K_{1} has extremal points at t = 0 and ±4/3.
The values −1, −0.75, and −0.5 have been proposed for α, motivated by various notions of optimality [9]. The choice α = −0.5 is particularly compelling. Bicubic interpolation is thirdorder accurate with α = −0.5 and only firstorder accurate for any other α. Furthermore, α = −0.5 is optimal in a sense of matching the sinc kernel and is also optimal in terms of E_{int}. All examples and comparisons with bicubic here use α = −0.5.
Comparison of the nearest neighbor, bilinear, and bicubic kernels and the sinc. 
Sinc
The sinc function. 
Surface plot of sinc(x)sinc(y). Yellow circles indicate local maxima. 
The Whittaker–Shannon interpolation [1] of v is
where sinc(t) := sin(πt)/(πt) for t ≠ 0 and sinc(0) := 1 is the interpolation kernel. The interpolation is such that u is bandlimited, that is, its Fourier transform is zero for frequencies outside of [−½,½]. This interpolation is also often called “sinc interpolation” or “Fourier zeropadding.” The extremal points of the sinc kernel are solutions of the equation πt = tan(πt), which include t ≈ 0, ±1.4303, ±2.4590, ±3.4709, ±4.4774.
The powerful property of sinc interpolation is that it is exact for bandlimited functions, that is, if v_{m,n} = f(m,n) with
then the interpolation reproduces f. Proof
In some ways, sinc interpolation is the ultimate interpolation. It is exact for bandlimited functions, so the method is very accurate on smooth data. Additionally, Fourier zeropadding avoids staircase artifacts, it is effective in reconstructing features at different orientations. Under an aliasing condition [10], Fourier interpolation reproduces cylindrical functions f(x,y) = h(αx + βy).
The disadvantage of sinc interpolation is that in aliased images it produces significant ripple artifacts (Gibbs phenomenon) in the vicinity of image edges. This is because sinc(x) decays slowly, at a rate of 1/x, so the damage from meeting an edge is spread throughout the image. Moreover, bandlimitedness can be a distorted view of reality. Thévenaz et al. [7] give the following amusing example: consider the air/matter interface of a patient in a CT scan, then according to classical physics, this is an abrupt discontinuity and cannot be expressed as a bandlimited function. Antialias filtering is not possible on physical matter, and any attempt to do so would probably be harmful to the patient.
Nearest Neighbor  Bilinear  Bicubic  Lanczos3  Cubic BSpline  Sinc 
Linear interpolation of a step edge: a balance between staircase artifacts and ripples. 
Windowed Sinc Approximations
The Lanczos3 kernel. 
Surface plot of the Lanczos3 kernel. Yellow circles indicate local maxima. 
A solution to limiting the ripple artifacts of the sinc kernel is to approximate it with a compactlysupported function
where w is a window function. A popular choice in image processing is the Lanczos window,
where n is a positive integer usually set to 2 or 3. The Lanczos kernels do not reproduce constants exactly, but can be normalized to fix this as described in the section on Kernel Normalization. There are many other possibilities for the window, for example Hamming, Kaiser, and Dolph–Chebyshev windows to name a few, each making different tradeoffs in frequency characteristics.
Comparison of the normalized Lanczos kernels and the sinc. 
Splines
Define the Bsplines,
The Bspline functions are optimal in the sense that, among all piecewise polynomials with uniformly spaced knots, the Bsplines have the maximal approximation order and are maximally continuous for a given support.
The Bspline β^{J−1} has approximation order J. For J > 2, β^{J−1} is not interpolating, so prefilting must be applied as described in the section on TwoStep Interpolation where β^{J−1} takes the role of φ.
As J→∞, Bspline interpolation converges to Whittaker–Shannon interpolation in a strong sense: the associated interpolation kernel K_{1} converges to the sinc both in the spatial and Fourier domains in L^{p} for any 1 ≤ p < ∞, see [4], [6]. This convergence is illustrated in the figure below with K_{1} and its Fourier transform for degrees 1, 3, 5, and 7.
Comparison of the Bspline interpolation kernels and the sinc. 
Aside from the Bsplines, there are numerous other piecewise polynomial methods for interpolation in the literature. One example is the class of Maximum Order and Minimal Support (Moms) functions introduced by Blu, Thévenaz, and Unser [8], which are linear combinations of β^{J−1} and its derivatives,
The Bsplines are a special case within the class of Moms. The Bsplines are the Moms with maximum regularity, they are (J − 2) times continuously differentiable.
Within the Moms, the “oMoms functions” are the Moms defined by minimizing the quantity
The oMoms functions have the same approximation order and support as the Bsplines, but the C_{int} constant is smaller. So compared to the Bsplines, the advantage of oMoms is that they have lower asymptotic L^{2} interpolation error. On the other hand, oMoms are less regular than the Bsplines.
The splines proposed by Schaum [5] are also within the class of Moms. Schaum splines have the property that they are interpolating, so prefiltering is not needed. However, for a given support size, the approximation constant C_{int} is worse than with the oMoms of the same support.
The first few evenorder oMoms are
Order J  φ(t)  
2  
4  
6  
8  
10  
12 
A note on numbering: Usually, with a piecewise polynomial method, its approximation order is L but its highest degree is (L − 1). We refer to methods by degree, for instance “oMoms 3” and “β^{3}” are methods that are locally cubic polynomial and have approximation order 4.
Radial Basis Functions
A radial basis function (RBF) is a function φ such that
RBF interpolation is useful for scattered data interpolation. Let (x_{n},y_{n}) denote the location of the nth sample, then the interpolant is
where the coefficients c_{n} are solved such that u(x_{n},y_{n}) = v_{n}. If the samples are uniformly spaced, then the c_{n} can be found as described in the section on TwoStep Interpolation.
Popular choices for φ include
•  Thin plate splines  φ(r) = εr^{n}logεr, n even 
•  Bessel  φ(r) = J_{0}(εr) 
•  Gaussian  φ(r) = exp(−(εr)^{2}) 
•  Multiquadric  φ(r) = (1 + (εr)^{2})^{½} 
•  Inverse multiquadric  φ(r) = (1 + (εr)^{2})^{−½} 
•  Inverse quadratic  φ(r) = (1 + (εr)^{2})^{−1} 
where ε is a positive parameter controlling the shape. A challenging problem is that smaller ε improves interpolation quality but worsens the numerical conditioning of solving the c_{n}. A topic of interest beyond the scope of this article is RBF interpolation in the limit ε → 0. For example, [11] investigates numerically stable interpolation for small ε with Gaussian φ.
Methodology
Boundary handling
A technicality is how to handle the image boundaries. The interpolation formulas may refer to samples that are outside of the domain of definition, especially when interpolating at points near or on the image boundaries.
The usual approach is to extrapolate (pad) the input image. Several methods for doing this are
 Constant extension, …aaaabcdeeee…
 Halfsample symmetric extension, …cbaabcdeedc…
 Wholesample symmetric extension, …dcbabcdedcb…
Image scaling
Suppose v is an M×N image with samples v_{m,n}, m = 0,…M − 1, n = 0,…N − 1, which are located on the integer grid .
For scaling v by a (possibly noninteger) factor d, a choice for where to resample the function u is on a topleftanchored grid of points
with M ′ = [dM], N ′ = [dN], where [⋅] is a rounding convention.
Example topleftanchored grids. Open circles denote input samples and filled circles denote interpolation samples. 
Another choice is the centered grid of points
where s_{r,M,M ′} = (1/d − 1 + M − M ′ /d)/2. If d is integer, the offsets simplify to s = (1/d − 1)/2. Although more complicated, the centered grid has the advantage that interpolation with a symmetric kernel on this grid commutes with flipping the image.
Example centered grids. Open circles denote input samples and filled circles denote interpolation samples. 
For an odd integer scale factor and interpolation with a symmetric kernel, the two grids are related in a simple way through adding and removing border samples. Interpolation performed on the topleftanchored grid is converted to an interpolation on the centered grid by removing (d − 1)/2 rows and columns from the bottom and right borders and extrapolating (according to the boundary extension) the same number of rows and columns on the top and left.
Domain mapping
Aside from changing image resolution, another application is to deform or map the domain of the image. That is, the mapped coordinates (x',y') are related to the original coordinates (x,y) by
Let u(x,y) be an interpolation of v, then the deformed image is obtained by resampling as
For example, two simple mappings are
 Subpixel translation
 Rotation
Other applications include parallax correction, lens distortion correction, image registration, and texture mapping.
Algorithm
Linear interpolation amounts to evaluating the sum
Choices of algorithms to do this efficiently depend on K and the sampling locations. Keep in mind that the sum is conceptually over the infinite integer grid, so v_{m,n} should be replaced with an appropriate boundary extension when m,n is beyond the domain of the image.
Nearest neighbor, bilinear, and bicubic
Suppose the interpolation method has a separable kernel K(x,y) = K_{1}(x)K_{1}(y) with compact support, K_{1}(x) = 0 for x ≥ R. Then interpolation at point (x,y) can be implemented as
where there are L×L nonzero terms in the sum. The values of K_{1} can be reused so that K_{1} only needs to be evaluated for only 2_L_ different arguments.
For image scaling, it is possible to interpolate more efficiently by exploiting the separable structure of the sampling points. Let v_{m,n}, m = 0,…, M − 1, n = 0,…, N − 1 be the given image, which is to be interpolated on the grid of points (x_{m′},y_{n′}) for m′ = 0,…, M ′ − 1, n′ = 0,…, N ′ − 1. Interpolation is performed by first interpolating along the y direction,
for m = 0,…, M − 1, n′ = 0,…, N ′ − 1, and then along the x direction,
for m′ = 0,…, M ′ − 1, n′ = 0,…, N ′ − 1.
An efficient way to implement the one dimensional interpolations is to express them as multiplication with a sparse matrix:
A is a sparse matrix of size N ′ × N and B is a sparse matrix of size M ′ × M. Away from the boundaries, the matrix entries are
Near the boundaries, the matrix entries need to be adjusted according to the boundary extension. Define
Halfsample symmetric  
Wholesample symmetric  
Constant extension 
then the matrix entries near the boundaries are computed as
and similarly for B. This sparse matrix approach is used by the libswscale library in FFmpeg [12] for scaling video frames. It is especially wellsuited for video, as A and B only need to be constructed once to interpolate any number of frames.
If the interpolation locations are uniformly spaced with a rational period, t_{n′} = t_{0} + n′⋅a/b, then another approach is to compute onedimensional interpolation through convolution. Given n′, let s = b⌊n′/b⌋ and r be such that n′ = bs + r, then
Lanczos
Interpolation with the Lanczos kernels does not exactly reproduce constants. To fix this, the Lanczos kernels should be normalized as described in the section on Kernel Normalization.
Kernel normalization can be incorporated efficiently into the algorithms discussed in the previous section. For the interpolation of a single point (x,y), the normalized interpolation is
For the sparse matrix approach, kernel normalization is achieved by scaling each matrix row so that it sums to one. Similarly for the convolution approach, h^{r} should be scaled so that it sums to one.
Splines
For spline methods, implementation is as described in the section on TwoStep Interpolation. The method is determined by the choice of basis function φ, for example, φ may be a Bspline or an oMoms function.
Given a basis function φ, define p_{m} = φ(m). For the prefiltering step, the convolution inverse (p)^{−1} is needed. We look first at the cubic Bspline as an example. For the cubic Bspline, p(z) = (z + 4 + z^{−1})/6 and
where r = 3^{½} − 2. The first factor corresponds to the firstorder causal recursive filter
The left endpoint c'_{0} can be computed depending on the boundary extension,
Halfsample symmetric  
Wholesample symmetric 
Since r < 1, the terms of the sum are decaying, so in practice we only evaluate as many terms as are needed for the desired accuracy.
The second factor corresponds to the firstorder anticausal recursive filter
The right endpoint c''_{N−1} can be computed according to the boundary extension as
Halfsample symmetric  
Wholesample symmetric 
Finally, the constant scale factor is applied c_{n} = 6_c''_{}n_. The cost of prefiltering is linear in the number of pixels. Prefiltering can be computed inplace so that the prefiltered values overwrite the memory used by the input. For prefiltering in twodimensions, the image is first prefiltered along each column, then the columnprefiltered image is prefiltered along each row.
More generally, if φ is symmetric and compactly supported, then p(z) has the form
where the r_{j} and 1/r_{j} are the roots of the polynomial z^{−J}p(z). We enumerate the roots so that r_{j} < 1. The Ztransform of (p)^{−1} is
Therefore, prefiltering can be performed as a cascade of firstorder recursive filters using the same algorithm as for the cubic Bspline. The table below lists the roots and the constant scale factor for several Bsplines and oMoms.


Once the image has been prefiltered c = (p)^{−1} ∗ v, interpolation at a point (x,y) is computed as
This second step is essentially the same formula as in the section on nearest neighbor, bilinear, and bicubic, and the same algorithms can be used to perform this step. The only differences are that the prefiltered image c is used in place of v and φ is used in place of K_{1}.
Sinc
For sinc interpolation with an integer scale factor, interpolation can be implemented using the fast Fourier transform (FFT). This approach follows from observing that the sinc interpolant is the unique image that is bandlimited and agrees with the input data.
The input image is first padded to twice its size in each dimension with halfsymmetric extension and transformed with the FFT. The interpolation is then constructed in the Fourier domain by copying the transform coefficients of the padded input image for the lower frequencies and filling zeros for the higher frequencies. The final interpolation is obtained by inverse FFT and removal of the symmetric padding. For realvalued input data, some computational savings can be made by using a realtocomplex FFT and complextoreal inverse FFT to exploit the complex conjugate redundancy in the transform. The complexity of the interpolation is O(N log N) for N output pixels.
Implementation
This software is distributed under the terms of the simplified BSD license.
 source code zip tar.gz
 online documentation
Fourier transforms are implemented using the FFTW library.
Please see the readme.html
file or the online documentation for details.
Examples
Scaling smooth data
The smooth function cos((x^{2} + y^{2})/10) is sampled on the grid {0.5,1.5,…,15.5}×{15.5,14.5,…,15.5} and interpolated on the centered grid with scale factor 4. The figure below compares the interpolation result using various linear methods, along with the root mean squared errors (RMSE) between the interpolation and the exact image.
For this example, higherorder methods tend to perform better. The oMoms perform slightly better than the Bsplines of the same order. Sinc interpolation yields the best result.
Scaling piecewise smooth data
We repeat the previous experiment with piecewise smooth data. The image is discontinuous at the edges, so the bandlimited hypothesis of sinc interpolation is severely violated. Sinc interpolation and other higherorder methods produce strong artifacts near edges. The bicubic and Lanczos2 interpolations have the lowest RMSE in this case.
Texture mapping
A 3D scene of a plane and a sphere is texture mapped with an image. Texture mapping is an example of domain mapping, where in this case the image domain is mapped to the projected plane and sphere surfaces.
Texture  Nearest  Bicubic 
This mapping is derived through reverse ray tracing. For each pixel of the rendering of the 3D scene, a ray of light is virtually traced in reverse. The light ray begins from the eye of the observer, passes through a pixel of the screen, and continues into the 3D scene. The ray may then intersect the plane or the sphere (or nothing, if the ray continues into the background).
Reverse ray tracing. A ray of light is traced from the observer into the 3D scene. 
The ray's intersection point (X,Y,Z) with the 3D scene determines the color of its pixel. The plane is texture mapped as x = X mod M, y = Z mod N. The sphere is texture mapped as x = θ, y = φ, where θ and φ are the spherical angles of the intersection.
Image rotation
Image rotation is another application that for a general rotation angle requires interpolation. The choice of interpolation method significantly affects the quality of the output image. In this experiment, an image is successively rotated by 5 degrees increments using several linear methods. The initial image is rotated by 5 degrees, then the rotated image is rotated by another 5 degrees, and so on.
Nearest  Bilinear  Bicubic  β^{11} 
Rotation using nearest neighbor produces jagged artifacts along edges. Rotation with bilinear is much better but blurs the image. Rotation with bicubic is sharper, but blurring is still noticeable. Finally, the sharpest rotation among these results is with the 11th order Bspline β^{11}. By interpolating with a highorder Bspline, it closely approximates rotation via sinc interpolation.
This material is based upon work supported by the National Science Foundation under Award No. DMS1004694. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Work partially supported by the MISS project of Centre National d'Etudes Spatiales, the Office of Naval Research under grant N000149710839 and by the European Research Council, advanced grant “Twelve labours.”