Welcome to Blogs @ Andrew Qu
Blog Index
All blogs
Search results

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, ∂2z/∂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) = ∑∑ Aij xi yj   (i=0,3 j=0,3)   ------   (1)
Formulation

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

z = f(x,y) = [1 x x2 x3][A]{1 y y2 y3}    ------  (2)
where [A] is a matrix
         a00 a01 a02 a03 
[A] = [ a10 a11 a12 a13 ]                  ------  (3)
        a20 a21 a22 a23
        a30 a31 a32 a33
For a unit square of (0,0) to (1,1), the values of f(x,y) ∂f(x,y)/∂x, ∂f(x,y)/∂y, ∂2f(x,y)/∂x∂y are known.
Therefore,
f(x,y) = [1 x x2 x3][A]{1 y y2 y3}
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}
fx(x,y) = [0 1 2x 3x2][A]{1 y y2 y3}
fx(0,0) = [0 1 0 0][A]{1 0 0 0}
fx(0,1) = [0 1 0 0][A]{1 1 1 1}
fx(1,0) = [0 1 2 3][A]{1 0 0 0}
fx(1,1) = [0 1 2 3][A]{1 1 1 1}
fy(x,y) = [1 x x2 x3][A]{0 1 2y 3y2}
fy(0,0) = [1 0 0 0][A]{0 1 0 0}
fy(0,1) = [1 0 0 0][A]{0 1 2 3}
fy(1,0) = [1 1 1 1][A]{0 1 0 0}
fy(1,1) = [1 1 1 1][A]{0 1 2 3}
fxy(x,y) = [0 1 2x 3x2][A]{0 1 2y 3y2}
fxy(0,0) = [0 1 0 0][A]{0 1 0 0}
fxy(0,1) = [0 1 0 0][A]{0 1 2 3}
fxy(1,0) = [0 1 2 3][A]{0 1 0 0}
fxy(1,1) = [0 1 2 3][A]{0 1 2 3}
The above 16 equations can be rearranged in the following way :
f(0,0)  = [1 0 0 0][A]{1 0 0 0}
f(1,0)  = [1 1 1 1][A]{1 0 0 0}
fx(0,0) = [0 1 0 0][A]{1 0 0 0}
fx(1,0) = [0 1 2 3][A]{1 0 0 0}
Rearranged to
  f(0,0)       1 0 0 0         1 
{ f(1,0) } = [ 1 1 1 1 ] [A] { 0 }    ----- (4)
  fx(0,0)      0 1 0 0         0
  fx(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}
fx(0,1) = [0 1 0 0][A]{1 1 1 1}
fx(1,1) = [0 1 2 3][A]{1 1 1 1}
Rearranged to
  f(0,1)       1 0 0 0         1 
{ f(1,1) } = [ 1 1 1 1 ] [A] { 1 }   ----- (5)
  fx(0,1)      0 1 0 0         1
  fx(1,1)      0 1 2 3         1
fy(0,0)  = [1 0 0 0][A]{0 1 0 0}
fy(1,0)  = [1 1 1 1][A]{0 1 0 0}
fxy(0,0) = [0 1 0 0][A]{0 1 0 0}
fxy(1,0) = [0 1 2 3][A]{0 1 0 0}
Rearranged to
  fy(0,0)       1 0 0 0         0 
{ fy(1,0) } = [ 1 1 1 1 ] [A] { 1 }   ----- (6)
  fxy(0,0)      0 1 0 0         0
  fxy(1,0)      0 1 2 3         0
fy(0,1)  = [1 0 0 0][A]{0 1 2 3}
fy(1,1)  = [1 1 1 1][A]{0 1 2 3}
fxy(0,1) = [0 1 0 0][A]{0 1 2 3}
fxy(1,1) = [0 1 2 3][A]{0 1 2 3}
Rearranged to
  fy(0,1)       1 0 0 0         0 
{ fy(1,1) } = [ 1 1 1 1 ] [A] { 1 }   ----- (7)
  fxy(0,1)      0 1 0 0         2
  fxy(1,1)      0 1 2 3         3
Observing 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:
  f(0,0)  f(0,1)  fy(0,0)  fy(0,1)       1 0 0 0         1 1 0 0
[ f(1,0)  f(1,1)  fy(1,0)  fy(1,1) ] = [ 1 1 1 1 ] [A] [ 0 1 1 1 ]
  fx(0,0) fx(0,1) fxy(0,0) fxy(0,1)       0 1 0 0         0 1 0 2
  fx(1,0) fx(1,1) fxy(1,0) fxy(1,1)       0 1 2 3         0 1 0 3
Solve the above equations for [A],
         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  1
where
        f(0,0)  f(0,1)   fy(0,0)   fy(0,1)
[F] = [ f(1,0)  f(1,1)   fy(1,0)   fy(1,1) ]   ---- (9)
        fx(0,0) fx(0,1)  fxy(0,0)  fxy(0,1)
        fx(1,0) fx(1,1)  fxy(1,0)  fxy(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 ∂f2/∂x∂y. For example,

∂f2(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 x2 x3][A]{1 y y2 y3}
Since (x,y) are normalized coordinates between 0 and 1, for point (px,py) in image space,
x = px / dxs
y = (dys-py) / dys
Once 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


Ads from Googlee
Dr Li Anchor Profi
www.anchorprofi.de
Engineering anchorage plate design system
©Andrew Qu, 2014. All rights reserved. Code snippets may be used "AS IS" without any kind of warranty. DIY tips may be followed at your own risk.