*
You have my permission to use freely, as long
as you keep the attribution. - Ken Perlin
*

**
Note: this Bicubic.html file also
works as a legal Bicubic.java file.
If you save the source under that name, you can just run
javac on it.
**

I created this class because I needed to evaluate a fairly expensive function over x,y many times. It turned out that I was able to rewrite the function as a sum of a "low frequency" part that only varied slowly, and a much less expensive "high frequency" part that just contained the details.

I used the bicubic patches as follows: First I presampled the low frequency part on a coarse grid. During the inner loop of the evaluation, I used this bicubic approximation, and then added in the high frequency part. Through this approach I was able to speed up the computation by a large factor.

If you provide a 4x4 grid of values (ie: all the corners in a 3x3 arrangement of square tiles), this class creates an object that will interpolate a Catmull-Rom spline to give you the value within any point of the center tile.

If you want to create a spline surface, you can make a two dimensional array of such objects, arranged so that adjoining objects overlap by one grid point.

Bicubic(double[][] G)Given 16 values on the grid [-1,0,1,2] X [-1,0,1,2], calculate bicubic coefficients.double eval(double x, double y)Given a point in the square [0,1] X [0,1], return a value.

f(x,y) = X * M * G * M^{T}* Y^{T}, where:X = (xThe constructor^{3}x^{2}x 1) ,

Y = (y^{3}y^{2}y 1) ,

M is the Catmull-Rom basis matrix.Bicubic()calculates the matrix C = M * G * M^{T}The method

eval(x,y)calculates the value X * C * Y^{T}

*/ public class Bicubic { static double[][] M = { //Catmull-Rom basis matrix{-0.5, 1.5, -1.5, 0.5}, { 1 , -2.5, 2 ,-0.5}, {-0.5, 0 , 0.5, 0 }, { 0 , 1 , 0 , 0 } }; double[][] C = new double[4][4]; //coefficients matrixBicubic(double[][] G) { double[][] T = new double[4][4]; for (int i = 0 ; i < 4 ; i++) //T = G * Mfor (int j = 0 ; j < 4 ; j++) for (int k = 0 ; k < 4 ; k++) T[i][j] += G[i][k] * M[j][k]; for (int i = 0 ; i < 4 ; i++) //^{T}C = M * Tfor (int j = 0 ; j < 4 ; j++) for (int k = 0 ; k < 4 ; k++) C[i][j] += M[i][k] * T[k][j]; } double[] X3 = C[0], X2 = C[1], X1 = C[2], X0 = C[3]; public double eval(double x, double y) { return x*(x*(x*(y * (y * (y * X3[0] + X3[1]) + X3[2]) + X3[3]) + (y * (y * (y * X2[0] + X2[1]) + X2[2]) + X2[3])) + (y * (y * (y * X1[0] + X1[1]) + X1[2]) + X1[3])) + (y * (y * (y * X0[0] + X0[1]) + X0[2]) + X0[3]); } }