Editing Bicubic interpolation (section)

==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>