Saturday, March 20, 2010

Cubic spline interpolation

I spent a while today casting around for a simple, fast and just-good-enough cubic interpolation algorithm to use in the GeoTools library where it is to serve as a 'filter function' that can be used when rendering maps.

The ACM digitial library has a mind-boggling collection of references on interpolation methods, most of which were fascinating but also over the top for this application.

In the end I put together the following code based on a Wikipedia article about cubic Hermite splines. Here it is...



/**
* Cubic hermite spline interpolation. This is adapted from the description of the
* algorithm at: http://en.wikipedia.org/wiki/Cubic_Hermite_spline.
* Tangent caculations are done with simple finite differencing in the interests
* of speed.
* <p>
* The input arrays xi and yi contain the coordinates of four interpolation
* points defining three segments with the middle segment containing the point
* for which we seek an interpolated value.
*
* @param x x ordinate of the point for which we seek an interpolated value
* and which lies between xi[1] and xi[2]
* @param xi interpolation point x ordinates
* @param yi interpolation point y ordinates
*
* @return interpolated value for x
*/
public double cubic(double x, double[] xi, double[] yi) {
double span01 = xi[1] - xi[0];
double span12 = xi[2] - xi[1];
double span23 = xi[3] - xi[2];

double t = (x - xi[1]) / span12;
double t2 = t*t;
double t3 = t2 * t;

double m1 = 0.5 * ((yi[2] - yi[1]) / span12 + (yi[1] - yi[0]) / span01);
double m2 = 0.5 * ((yi[3] - yi[2]) / span23 + (yi[2] - yi[1]) / span12);

double y = (2*t3 - 3*t2 + 1) * yi[1] +
(t3 - 2*t2 + t) * span12 * m1 +
(-2*t3 + 3*t2) * yi[2] +
(t3 - t2) * span12 * m2;

return y;
}


I make no claims about the elegance or efficiency of this code other than suspecting that it rates badly in both areas, but it does a perfectly acceptable job when I throw some test data at it...