Editing Bicubic interpolation

{{Short description|Extension of cubic spline interpolation}}
{{comparison of 1D and 2D interpolation.svg}}
In [[mathematics]], '''bicubic interpolation''' is an extension of [[Cubic interpolation#Interpolation on a single interval|cubic spline interpolation]] (a method of applying cubic interpolation to a data set) for [[interpolation|interpolating]] data points on a [[two-dimensional]] [[regular grid]]. The interpolated surface (meaning the kernel shape, not the image) is [[Smooth function|smoother]] than corresponding surfaces obtained by [[bilinear interpolation]] or [[nearest-neighbor interpolation]]. Bicubic interpolation can be accomplished using either [[Lagrange polynomial]]s, [[cubic spline]]s, or [[#Bicubic convolution algorithm|cubic convolution]] algorithm.

In [[image processing]], bicubic interpolation is often chosen over bilinear or nearest-neighbor interpolation in [[resampling (bitmap)|image resampling]], when speed is not an issue. In contrast to bilinear interpolation, which only takes 4 [[pixel]]s (2×2) into account, bicubic interpolation considers 16 pixels (4×4). Images resampled with bicubic interpolation can have different interpolation [[Spatial anti-aliasing|artifacts]], depending on the b and c values chosen.

==Computation==

[[Image:Interpolation-bicubic.svg|thumb|right|Bicubic interpolation on the square <math>[0,4] \times [0,4]</math> consisting of 25 unit squares patched together. Bicubic interpolation as per [[Matplotlib]]'s implementation. Colour indicates function value. The black dots are the locations of the prescribed data being interpolated. Note how the color samples are not radially symmetric.]]
[[Image:Interpolation-bilinear.svg|thumb|right|[[Bilinear interpolation]] on the same dataset as above. Derivatives of the surface are not continuous over the square boundaries.]]
[[Image:Interpolation-nearest.svg|thumb|right|[[Nearest-neighbor interpolation]] on the same dataset as above.]]

Suppose the function values <math>f</math> and the derivatives <math>f_x</math>, <math>f_y</math> and <math>f_{xy}</math> are known at the four corners <math>(0,0)</math>, <math>(1,0)</math>, <math>(0,1)</math>, and <math>(1,1)</math> of the unit square. The interpolated surface can then be written as
<math display="block">p(x,y) = \sum\limits_{i=0}^3 \sum_{j=0}^3 a_{ij} x^i y^j.</math>

The interpolation problem consists of determining the 16 coefficients <math>a_{ij}</math>.
Matching <math>p(x,y)</math> with the function values yields four equations:
# <math>f(0,0)      = p(0,0)   = a_{00},</math>
# <math>f(1,0)      = p(1,0)   = a_{00} + a_{10} + a_{20} + a_{30},</math>
# <math>f(0,1)      = p(0,1)   = a_{00} + a_{01} + a_{02} + a_{03},</math>
# <math>f(1,1)      = p(1,1)   = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=0}^3 a_{ij}.</math>

Likewise, eight equations for the derivatives in the <math>x</math> and the <math>y</math> directions:
# <math>f_x(0,0)    = p_x(0,0) = a_{10},</math>
# <math>f_x(1,0)    = p_x(1,0) =  a_{10} + 2a_{20} + 3a_{30},</math>
# <math>f_x(0,1)    = p_x(0,1) = a_{10} + a_{11} + a_{12} + a_{13},</math>
# <math>f_x(1,1)    = p_x(1,1) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 a_{ij} i,</math>
# <math>f_y(0,0)    = p_y(0,0) = a_{01},</math>
# <math>f_y(1,0)    = p_y(1,0) = a_{01} + a_{11} + a_{21} + a_{31},</math>
# <math>f_y(0,1)    = p_y(0,1) = a_{01} + 2a_{02} + 3a_{03},</math>
# <math>f_y(1,1)    = p_y(1,1) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 a_{ij} j.</math>

And four equations for the <math>xy</math> [[mixed partial derivative]]:
# <math>f_{xy}(0,0) = p_{xy}(0,0) = a_{11},</math>
# <math>f_{xy}(1,0) = p_{xy}(1,0) = a_{11} + 2a_{21} + 3a_{31},</math>
# <math>f_{xy}(0,1) = p_{xy}(0,1) = a_{11} + 2a_{12} + 3a_{13},</math>
# <math>f_{xy}(1,1) = p_{xy}(1,1) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 a_{ij} i j.</math>

The expressions above have used the following identities:
<math display="block">p_x(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 a_{ij} i x^{i-1} y^j,</math>
<math display="block">p_y(x,y) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 a_{ij} x^i j y^{j-1},</math>
<math display="block">p_{xy}(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 a_{ij} i x^{i-1} j y^{j-1}.</math>

This procedure yields a surface <math>p(x,y)</math> on the [[unit square]] <math>[0,1] \times [0,1]</math> that is continuous and has continuous derivatives. Bicubic interpolation on an arbitrarily sized [[regular grid]] can then be accomplished by patching together such bicubic surfaces, ensuring that the derivatives match on the boundaries.

Grouping the unknown parameters <math>a_{ij}</math> in a vector
<math display="block">\alpha=\left[\begin{smallmatrix}a_{00}&a_{10}&a_{20}&a_{30}&a_{01}&a_{11}&a_{21}&a_{31}&a_{02}&a_{12}&a_{22}&a_{32}&a_{03}&a_{13}&a_{23}&a_{33}\end{smallmatrix}\right]^T</math>
and letting
<math display="block">x=\left[\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&f_x(0,0)&f_x(1,0)&f_x(0,1)&f_x(1,1)&f_y(0,0)&f_y(1,0)&f_y(0,1)&f_y(1,1)&f_{xy}(0,0)&f_{xy}(1,0)&f_{xy}(0,1)&f_{xy}(1,1)\end{smallmatrix}\right]^T,</math>
the above system of equations can be reformulated into a matrix for the linear equation <math>A\alpha=x</math>.

Inverting the matrix gives the more useful linear equation <math>A^{-1}x=\alpha</math>, where
<math display="block">A^{-1}=\left[\begin{smallmatrix}\begin{array}{rrrrrrrrrrrrrrrr}
 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -3 & 3 & 0 & 0 & -2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 2 & -2 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 3 & 0 & 0 & -2 & -1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & -2 & 0 & 0 & 1 & 1 & 0 & 0 \\
 -3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & -3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & -1 & 0 \\
 9 & -9 & -9 & 9 & 6 & 3 & -6 & -3 & 6 & -6 & 3 & -3 & 4 & 2 & 2 & 1 \\
 -6 & 6 & 6 & -6 & -3 & -3 & 3 & 3 & -4 & 4 & -2 & 2 & -2 & -2 & -1 & -1 \\
 2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\
 -6 & 6 & 6 & -6 & -4 & -2 & 4 & 2 & -3 & 3 & -3 & 3 & -2 & -1 & -2 & -1 \\
 4 & -4 & -4 & 4 & 2 & 2 & -2 & -2 & 2 & -2 & 2 & -2 & 1 & 1 & 1 & 1
\end{array}\end{smallmatrix}\right],</math>
which allows <math>\alpha</math> to be calculated quickly and easily.

There can be another concise matrix form for 16 coefficients:
<math display="block">\begin{bmatrix}
f(0,0)&f(0,1)&f_y (0,0)&f_y (0,1)\\f(1,0)&f(1,1)&f_y (1,0)&f_y (1,1)\\f_x (0,0)&f_x (0,1)&f_{xy} (0,0)&f_{xy} (0,1)\\f_x (1,0)&f_x (1,1)&f_{xy} (1,0)&f_{xy} (1,1)
\end{bmatrix} =
\begin{bmatrix}
1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&3
\end{bmatrix}
\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
\begin{bmatrix}
1&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3
\end{bmatrix},</math>
or
<math display="block">
\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
=
\begin{bmatrix}
1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1
\end{bmatrix}
\begin{bmatrix}
f(0,0)&f(0,1)&f_y (0,0)&f_y (0,1)\\f(1,0)&f(1,1)&f_y (1,0)&f_y (1,1)\\f_x (0,0)&f_x (0,1)&f_{xy} (0,0)&f_{xy} (0,1)\\f_x (1,0)&f_x (1,1)&f_{xy} (1,0)&f_{xy} (1,1)
\end{bmatrix}
\begin{bmatrix}
1&0&-3&2\\0&0&3&-2\\0&1&-2&1\\0&0&-1&1
\end{bmatrix},
</math>
where
<math display="block">p(x,y)=\begin{bmatrix}1
&x&x^2&x^3\end{bmatrix}
\begin{bmatrix}
a_{00}&a_{01}&a_{02}&a_{03}\\a_{10}&a_{11}&a_{12}&a_{13}\\a_{20}&a_{21}&a_{22}&a_{23}\\a_{30}&a_{31}&a_{32}&a_{33}
\end{bmatrix}
\begin{bmatrix}1\\y\\y^2\\y^3\end{bmatrix}.</math>

==Extension to rectilinear grids==

Often, applications call for bicubic interpolation using data on a rectilinear grid, rather than the unit square. In this case, the identities for <math>p_x, p_y,</math> and <math>p_{xy}</math> become
<math display="block">p_x(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=0}^3 \frac{a_{ij} i x^{i-1} y^j}{\Delta x},</math>
<math display="block">p_y(x,y) = \textstyle \sum\limits_{i=0}^3 \sum\limits_{j=1}^3 \frac{a_{ij} x^i j y^{j-1}}{\Delta y},</math>
<math display="block">p_{xy}(x,y) = \textstyle \sum\limits_{i=1}^3 \sum\limits_{j=1}^3 \frac{a_{ij} i x^{i-1} j y^{j-1}}{\Delta x \Delta y},</math>
where <math>\Delta x</math> is the <math>x</math> spacing of the cell containing the point <math>(x,y)</math> and similar for <math>\Delta y</math>.
In this case, the most practical approach to computing the coefficients <math>\alpha</math> is to let
<math display="block">x=\left[\begin{smallmatrix}f(0,0)&f(1,0)&f(0,1)&f(1,1)&\Delta x f_x(0,0)&\Delta xf_x(1,0)&\Delta x f_x(0,1)&\Delta x f_x(1,1)&\Delta y f_y(0,0)&\Delta y f_y(1,0)&\Delta y f_y(0,1)&\Delta y f_y(1,1)&\Delta x \Delta y f_{xy}(0,0)&\Delta x \Delta y f_{xy}(1,0)&\Delta x \Delta y f_{xy}(0,1)&\Delta x \Delta y f_{xy}(1,1)\end{smallmatrix}\right]^T,</math>
then to solve <math>\alpha=A^{-1}x</math> with <math>A</math> as before. Next, the normalized interpolating variables are computed as
<math display="block">\begin{align}
\overline{x} &= \frac{x-x_0}{x_1-x_0}, \\
\overline{y} &= \frac{y-y_0}{y_1-y_0}
\end{align}</math>
where <math>x_0, x_1, y_0,</math> and <math>y_1</math> are the <math>x</math> and <math>y</math> coordinates of the grid points surrounding the point <math>(x,y)</math>. Then, the interpolating surface becomes
<math display="block">p(x,y) = \sum\limits_{i=0}^3 \sum_{j=0}^3 a_{ij} {\overline{x}}^i {\overline{y}}^j.</math>

==Finding derivatives from function values==

If the derivatives are unknown, they are typically approximated from the function values at points neighbouring the corners of the unit square, e.g. using [[finite differences]].

To find either of the single derivatives, <math>f_x</math> or <math>f_y</math>, using that method, find the slope between the two ''surrounding'' points in the appropriate axis. For example, to calculate <math>f_x</math> for one of the points, find <math>f(x,y)</math> for the points to the left and right of the target point and calculate their slope, and similarly for <math>f_y</math>.

To find the cross derivative <math>f_{xy}</math>, take the derivative in both axes, one at a time. For example, one can first use the <math>f_x</math> procedure to find the <math>x</math> derivatives of the points above and below the target point, then use the <math>f_y</math> procedure on those values (rather than, as usual, the values of <math>f</math> for those points) to obtain the value of <math>f_{xy}(x,y)</math> for the target point.  (Or one can do it in the opposite direction, first calculating <math>f_y</math> and then <math>f_x</math> from those.  The two give equivalent results.)

At the edges of the dataset, when one is missing some of the surrounding points, the missing points can be approximated by a number of methods. A simple and common method is to assume that the slope from the existing point to the target point continues without further change, and using this to calculate a hypothetical value for the missing point.

==Bicubic convolution algorithm==

Bicubic spline interpolation requires the solution of the linear system described above for each grid cell. An interpolator with similar properties can be obtained by applying a [[convolution]] with the following kernel in both dimensions:
<math display="block">W(x) = 
\begin{cases}
 (a+2)|x|^3-(a+3)|x|^2+1 & \text{for } |x| \leq 1, \\
 a|x|^3-5a|x|^2+8a|x|-4a & \text{for } 1 < |x| < 2, \\
 0                       & \text{otherwise},
\end{cases}
</math>
where <math>a</math> is usually set to −0.5 or −0.75. Note that <math>W(0)=1</math> and <math>W(n)=0</math> for all nonzero integers <math>n</math>.

This approach was proposed by Keys, who showed that <math>a=-0.5</math> produces third-order convergence with respect to the sampling interval of the original function.<ref name=Keys>{{cite journal
 | author = R. Keys
 | year = 1981
 | title = Cubic convolution interpolation for digital image processing
 | journal = IEEE Transactions on Acoustics, Speech, and Signal Processing
 | doi = 10.1109/TASSP.1981.1163711
 | volume = 29
 | pages = 1153–1160
 | issue = 6
| bibcode = 1981ITASS..29.1153K
 | citeseerx = 10.1.1.320.776
 }}</ref>

If we use the matrix notation for the common case <math>a = -0.5</math>, we can express the equation in a more friendly manner:
<math display="block">p(t) =
\tfrac{1}{2}
\begin{bmatrix}
1 & t & t^2 & t^3
\end{bmatrix}

\begin{bmatrix}
0 & 2 & 0 & 0 \\
-1 & 0 & 1 & 0 \\
2 & -5 & 4 & -1 \\
-1 & 3 & -3 & 1
\end{bmatrix}

\begin{bmatrix}
f_{-1} \\
f_0 \\
f_1 \\
f_2
\end{bmatrix}
</math>
for <math>t</math> between 0 and 1 for one dimension. Note that for 1-dimensional cubic convolution interpolation 4 sample points are required. For each inquiry two samples are located on its left and two samples on the right. These points are indexed from −1 to 2 in this text. The distance from the point indexed with 0 to the inquiry point is denoted by <math>t</math> here.

For two dimensions first applied once in <math>x</math> and again in <math>y</math>:
<math display="block">\begin{align}
b_{-1}&= p(t_x, f_{(-1,-1)}, f_{(0,-1)}, f_{(1,-1)}, f_{(2,-1)}), \\[1ex]
b_{0} &= p(t_x, f_{(-1,0)}, f_{(0,0)}, f_{(1,0)}, f_{(2,0)}), \\[1ex]
b_{1} &= p(t_x, f_{(-1,1)}, f_{(0,1)}, f_{(1,1)}, f_{(2,1)}), \\[1ex]
b_{2} &= p(t_x, f_{(-1,2)}, f_{(0,2)}, f_{(1,2)}, f_{(2,2)}),
\end{align}</math>
<math display="block">p(x,y) = p(t_y, b_{-1}, b_{0}, b_{1}, b_{2}).</math>

==Use in computer graphics==

[[Image:Accutance.svg|250px|thumb|The lower half of this figure is a magnification of the upper half, showing how the apparent sharpness of the left-hand line is created. Bicubic interpolation causes overshoot, which increases [[acutance]].]]<!--Don't scale the image, it will greatly reduce the effect-->

The bicubic algorithm is frequently used for scaling images and video for display (see [[Resampling (bitmap)|bitmap resampling]]). It preserves fine detail better than the common [[bilinear filtering|bilinear]] algorithm.

However, due to the negative lobes on the kernel, it causes [[overshoot (signal)|overshoot]] (haloing). This can cause [[Clipping (signal processing)|clipping]], and is an artifact (see also [[ringing artifacts]]), but it increases [[acutance]] (apparent sharpness), and can be desirable.

==See also==

{{portal|Mathematics}}

* [[Spatial anti-aliasing]]
* [[Bézier surface]]
* [[Bilinear interpolation]]
* [[Cubic Hermite spline]], the one-dimensional analogue of bicubic spline
* [[Lanczos resampling]]
* [[Natural neighbor interpolation]]
* [[Sinc filter]]
* [[Spline interpolation]]
* [[Tricubic interpolation]]
* [[Directional Cubic Convolution Interpolation]]

==References==
{{reflist}}

==External links==
* [https://web.archive.org/web/20051024202307/http://www.geovista.psu.edu/sites/geocomp99/Gc99/082/gc_082.htm Application of interpolation to elevation samples]
* [https://sepwww.stanford.edu/data/media/public/docs/sep107/paper_html/node20.html Interpolation theory]
* [http://www.paulinternet.nl/?page=bicubic Explanation and Java/C++ implementation of (bi)cubic interpolation]
* [http://mathformeremortals.wordpress.com/2013/01/15/bicubic-interpolation-excel-worksheet-function/ Excel Worksheet Function for Bicubic Lagrange Interpolation]

{{DEFAULTSORT:Bicubic Interpolation}}
[[Category:Image processing]]
[[Category:Multivariate interpolation]]