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...

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...

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...

## Comments

## Post a Comment