## Bicubic Surface Interpolation

**Summary**

This blog describes a formulation of the bicubic interpolation method. The bicubic method aims to produce a best fit surface over a 2D area when the function values are known only at a set of discrete points. Application to image processing is also demonstrated.

**Problem Definition**

Assume Æ is a bounded 2D area, for example, a rectangle. We place a regular grid over Æ, as shown below,

There is a function z = f(x,y), for which, at the set of grid points, the function values and partial derivatives
∂z/∂x, ∂z/∂y, ∂^{2}z/∂x∂y,
are known. What's the best fit function value at any given point (x, y) ∈ Æ?

The bicubic interpolation aims to solve this problem by using the full cubic function to represent the surface, that is,

f(x, y) = ∑∑ A_{ij}x^{i}y^{j}(i=0,3 j=0,3) ------ (1)

**Formulation**

Rearrange equation (1) in matrix form, as follows:

z = f(x,y) = [1 x xwhere [A] is a matrix^{2}x^{3}][A]{1 y y^{2}y^{3}} ------ (2)

aFor a unit square of (0,0) to (1,1), the values of f(x,y) ∂f(x,y)/∂x, ∂f(x,y)/∂y, ∂_{00}a_{01}a_{02}a_{03}[A] = [ a_{10}a_{11}a_{12}a_{13}] ------ (3) a_{20}a_{21}a_{22}a_{23}a_{30}a_{31}a_{32}a_{33}

^{2}f(x,y)/∂x∂y are known.

Therefore,

f(x,y) = [1 x x^{2}x^{3}][A]{1 y y^{2}y^{3}} f(0,0) = [1 0 0 0][A]{1 0 0 0} f(0,1) = [1 0 0 0][A]{1 1 1 1} f(1,0) = [1 1 1 1][A]{1 0 0 0} f(1,1) = [1 1 1 1][A]{1 1 1 1}

f_{x}(x,y) = [0 1 2x 3x^{2}][A]{1 y y^{2}y^{3}} f_{x}(0,0) = [0 1 0 0][A]{1 0 0 0} f_{x}(0,1) = [0 1 0 0][A]{1 1 1 1} f_{x}(1,0) = [0 1 2 3][A]{1 0 0 0} f_{x}(1,1) = [0 1 2 3][A]{1 1 1 1}

f_{y}(x,y) = [1 x x^{2}x^{3}][A]{0 1 2y 3y^{2}} f_{y}(0,0) = [1 0 0 0][A]{0 1 0 0} f_{y}(0,1) = [1 0 0 0][A]{0 1 2 3} f_{y}(1,0) = [1 1 1 1][A]{0 1 0 0} f_{y}(1,1) = [1 1 1 1][A]{0 1 2 3}

fThe above 16 equations can be rearranged in the following way :_{xy}(x,y) = [0 1 2x 3x^{2}][A]{0 1 2y 3y^{2}} f_{xy}(0,0) = [0 1 0 0][A]{0 1 0 0} f_{xy}(0,1) = [0 1 0 0][A]{0 1 2 3} f_{xy}(1,0) = [0 1 2 3][A]{0 1 0 0} f_{xy}(1,1) = [0 1 2 3][A]{0 1 2 3}

f(0,0) = [1 0 0 0][A]{1 0 0 0} f(1,0) = [1 1 1 1][A]{1 0 0 0} fRearranged to_{x}(0,0) = [0 1 0 0][A]{1 0 0 0} f_{x}(1,0) = [0 1 2 3][A]{1 0 0 0}

f(0,0) 1 0 0 0 1 { f(1,0) } = [ 1 1 1 1 ] [A] { 0 } ----- (4) f_{x}(0,0) 0 1 0 0 0 f_{x}(1,0) 0 1 2 3 0

f(0,1) = [1 0 0 0][A]{1 1 1 1} f(1,1) = [1 1 1 1][A]{1 1 1 1} fRearranged to_{x}(0,1) = [0 1 0 0][A]{1 1 1 1} f_{x}(1,1) = [0 1 2 3][A]{1 1 1 1}

f(0,1) 1 0 0 0 1 { f(1,1) } = [ 1 1 1 1 ] [A] { 1 } ----- (5) f_{x}(0,1) 0 1 0 0 1 f_{x}(1,1) 0 1 2 3 1

fRearranged to_{y}(0,0) = [1 0 0 0][A]{0 1 0 0} f_{y}(1,0) = [1 1 1 1][A]{0 1 0 0} f_{xy}(0,0) = [0 1 0 0][A]{0 1 0 0} f_{xy}(1,0) = [0 1 2 3][A]{0 1 0 0}

f_{y}(0,0) 1 0 0 0 0 { f_{y}(1,0) } = [ 1 1 1 1 ] [A] { 1 } ----- (6) f_{xy}(0,0) 0 1 0 0 0 f_{xy}(1,0) 0 1 2 3 0

fRearranged to_{y}(0,1) = [1 0 0 0][A]{0 1 2 3} f_{y}(1,1) = [1 1 1 1][A]{0 1 2 3} f_{xy}(0,1) = [0 1 0 0][A]{0 1 2 3} f_{xy}(1,1) = [0 1 2 3][A]{0 1 2 3}

fObserving that the 4x4 matrix on the left of [A] in equations 4,5,6,7 are identical. Therefore, the columuns in these equations can be arranged into a single 4x4 matrix as follows:_{y}(0,1) 1 0 0 0 0 { f_{y}(1,1) } = [ 1 1 1 1 ] [A] { 1 } ----- (7) f_{xy}(0,1) 0 1 0 0 2 f_{xy}(1,1) 0 1 2 3 3

f(0,0) f(0,1) fSolve the above equations for [A],_{y}(0,0) f_{y}(0,1) 1 0 0 0 1 1 0 0 [ f(1,0) f(1,1) f_{y}(1,0) f_{y}(1,1) ] = [ 1 1 1 1 ] [A] [ 0 1 1 1 ] f_{x}(0,0) f_{x}(0,1) f_{xy}(0,0) f_{xy}(0,1) 0 1 0 0 0 1 0 2 f_{x}(1,0) f_{x}(1,1) f_{xy}(1,0) f_{xy}(1,1) 0 1 2 3 0 1 0 3

1 0 0 0 1 0 -3 2 [A] = [ 0 0 1 0 ] [F] [ 0 0 3 -2 ] ----- (8) -3 3 -2 -1 0 1 -2 1 2 -2 1 1 0 0 -1 1where

f(0,0) f(0,1) f_{y}(0,0) f_{y}(0,1) [F] = [ f(1,0) f(1,1) f_{y}(1,0) f_{y}(1,1) ] ---- (9) 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)

**Application to Imaging**

In the following diagram, the function values at all grid points are known. But the derivatives are not. We need to first find the approximate derivative values at all nodes.

Consider node (i,j), the derivatives can be approximated using the ajacent node values as follows:

∂f(i,j)/∂x = (f(i+1,j) - f(i-1,j)) / 2Δx; ∂f(i,j)/∂y = (f(i,j+1) - f(i,j-1)) / 2Δy;For an edge node, we only use node value - the inner side node value. For example,

Left edge: ∂f(0,j)/∂x = (f(1,j) - f(0,j)) / Δx; Right edge: ∂f(n,j)/∂y = (f(n,j) - f(n-1,j)) / Δy;Δx, Δy above are the intervals in x and y direction. Since we are using normalized coordinates for equatin (9), these can be set to Δx = Δy = 1

Using the same rule, we can work out ∂f^{2}/∂x∂y. For example,

∂f^{2}(i,j)/∂x∂y = (∂f(i+1,j)/∂x - ∂f(i-1,j)/∂x) / 2Δy;

Assume the rectangle bounded by 4 corner points (i, i+1, j, j+1) is mapped to the image space (Xs,Ys). dxs is the number of pixels in x-dir. when mapped to the screen space. dys is the number of pixels in y-dir.

We need to find the values for all pixels within the rectangle for which matrix [A] is a constant and calculated once using equation (8). Once a is calculated, function value is calculate using equation (2), that is

f(x,y) = [1 x xSince (x,y) are normalized coordinates between 0 and 1, for point (px,py) in image space,^{2}x^{3}][A]{1 y y^{2}y^{3}}

x = px / dxs y = (dys-py) / dysOnce a function value is calculated, we can assign it a color value and set to the pixel. As an example, one simple way is to asign red (255,0,0) color to max value and blue (0,0,255) color to the min value. Any thing between uses linear interpolation. Then our colur scheme is

red = 255 * (f(px,py) - fMin) / (fMax - fMin); green = 0 blue = 255 * (fMax - f(px,py) ) / (fMax - fMin);

END