Bicubic interpolation explained

In mathematics, bicubic interpolation is an extension of cubic spline interpolation (a method of applying cubic interpolation to a data set) for interpolating data points on a two-dimensional regular grid. The interpolated surface (meaning the kernel shape, not the image) is smoother than corresponding surfaces obtained by bilinear interpolation or nearest-neighbor interpolation. Bicubic interpolation can be accomplished using either Lagrange polynomials, cubic splines, or cubic convolution algorithm.

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

Computation

Suppose the function values

and the derivatives

f_x

f_y

and

f_xy

are known at the four corners

(0,0)

(1,0)

(0,1)

, and

(1,1)

of the unit square. The interpolated surface can then be written as

p(x,y) = \sum\limits_^3 \sum_^3 a_ x^i y^j.

The interpolation problem consists of determining the 16 coefficients

a_ij

.Matching

p(x,y)

with the function values yields four equations:

f(0,0)=p(0,0)=a₀₀,

f(1,0)=p(1,0)=a₀₀+a₁₀+a₂₀+a₃₀,

f(0,1)=p(0,1)=a₀₀+a₀₁+a₀₂+a₀₃,

f(1,1)=p(1,1)=

	3
style\sum\limits
	i=0

	3
\sum\limits
	j=0

a_ij.

Likewise, eight equations for the derivatives in the

and the

directions:

f_x(0,0)=p_x(0,0)=a₁₀,

f_x(1,0)=p_x(1,0)=a₁₀+2a₂₀+3a₃₀,

f_x(0,1)=p_x(0,1)=a₁₀+a₁₁+a₁₂+a₁₃,

f_x(1,1)=p_x(1,1)=

	3
style\sum\limits
	i=1

	3
\sum\limits
	j=0

a_iji,

f_y(0,0)=p_y(0,0)=a₀₁,

f_y(1,0)=p_y(1,0)=a₀₁+a₁₁+a₂₁+a₃₁,

f_y(0,1)=p_y(0,1)=a₀₁+2a₀₂+3a₀₃,

f_y(1,1)=p_y(1,1)=

	3
style\sum\limits
	i=0

	3
\sum\limits
	j=1

a_ijj.

And four equations for the

mixed partial derivative:

f_xy(0,0)=p_xy(0,0)=a₁₁,

f_xy(1,0)=p_xy(1,0)=a₁₁+2a₂₁+3a₃₁,

f_xy(0,1)=p_xy(0,1)=a₁₁+2a₁₂+3a₁₃,

f_xy(1,1)=p_xy(1,1)=

	3
style\sum\limits
	i=1

	3
\sum\limits
	j=1

a_ijij.

The expressions above have used the following identities: $p_x(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 a_ i x^ y^j,$ $p_y(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 a_ x^i j y^,$ $p_(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 a_ i x^ j y^.$

This procedure yields a surface

p(x,y)

on the unit square

[0,1] x [0,1]

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

a_ij

in a vector

\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

and letting

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,

the above system of equations can be reformulated into a matrix for the linear equation

A\alpha=x

Inverting the matrix gives the more useful linear equation

A^-1x=\alpha

, where

A^=\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],

which allows

\alpha

to be calculated quickly and easily.

There can be another concise matrix form for 16 coefficients: $\beginf(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_ (0,0)&f_ (0,1)\\f_x (1,0)&f_x (1,1)&f_ (1,0)&f_ (1,1)\end =\begin1&0&0&0\\1&1&1&1\\0&1&0&0\\0&1&2&3\end\begina_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\end\begin1&1&0&0\\0&1&1&1\\0&1&0&2\\0&1&0&3\end,$ or $\begina_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\end=\begin1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end\beginf(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_ (0,0)&f_ (0,1)\\f_x (1,0)&f_x (1,1)&f_ (1,0)&f_ (1,1)\end\begin1&0&-3&2\\0&0&3&-2\\0&1&-2&1\\0&0&-1&1\end,$ where $p(x,y)=\begin1&x&x^2&x^3\end\begina_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\\a_&a_&a_&a_\end\begin1\\y\\y^2\\y^3\end.$

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

p_x,p_y,

and

p_xy

become

p_x(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 \frac,

p_y(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 \frac,

p_(x,y) = \textstyle \sum\limits_^3 \sum\limits_^3 \frac,

where

\Deltax

is the

spacing of the cell containing the point

(x,y)

and similar for

\Deltay

.In this case, the most practical approach to computing the coefficients

\alpha

is to let

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,

then to solve

\alpha=A^-1x

with

as before. Next, the normalized interpolating variables are computed as

\begin\overline &= \frac, \\\overline &= \frac\end

where

x_0,x_1,y_0,

and

y₁

are the

and

coordinates of the grid points surrounding the point

(x,y)

. Then, the interpolating surface becomes

p(x,y) = \sum\limits_^3 \sum_^3 a_ ^i ^j.

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,

f_x

f_y

, using that method, find the slope between the two surrounding points in the appropriate axis. For example, to calculate

f_x

for one of the points, find

f(x,y)

for the points to the left and right of the target point and calculate their slope, and similarly for

f_y

To find the cross derivative

f_xy

, take the derivative in both axes, one at a time. For example, one can first use the

f_x

procedure to find the

derivatives of the points above and below the target point, then use the

f_y

procedure on those values (rather than, as usual, the values of

for those points) to obtain the value of

f_xy(x,y)

for the target point. (Or one can do it in the opposite direction, first calculating

f_y

and then

f_x

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: $W(x) = \begin (a+2)|x|^3-(a+3)|x|^2+1 & \text |x| \leq 1, \\ a|x|^3-5a|x|^2+8a|x|-4a & \text 1 < |x| < 2, \\ 0 & \text,\end$ where

is usually set to −0.5 or −0.75. Note that

W(0)=1

and

W(n)=0

for all nonzero integers

This approach was proposed by Keys, who showed that

a=-0.5

produces third-order convergence with respect to the sampling interval of the original function.^[1]

If we use the matrix notation for the common case

a=-0.5

, we can express the equation in a more friendly manner:

p(t) =\tfrac\begin1 & t & t^2 & t^3\end\begin0 & 2 & 0 & 0 \\-1 & 0 & 1 & 0 \\2 & -5 & 4 & -1 \\-1 & 3 & -3 & 1\end
\beginf_ \\f_0 \\f_1 \\f_2\end

for

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

here.

For two dimensions first applied once in

and again in

\beginb_&= p(t_x, f_, f_, f_, f_), \\[1ex]b_ &= p(t_x, f_, f_, f_, f_), \\[1ex]b_ &= p(t_x, f_, f_, f_, f_), \\[1ex]b_ &= p(t_x, f_, f_, f_, f_),\end

p(x,y) = p(t_y, b_, b_, b_, b_).

Use in computer graphics

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

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

External links

Notes and References

R. Keys . 1981 . Cubic convolution interpolation for digital image processing . IEEE Transactions on Acoustics, Speech, and Signal Processing . 10.1109/TASSP.1981.1163711 . 29 . 1153–1160 . 6. 1981ITASS..29.1153K . 10.1.1.320.776 .