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

Comments

Circle packing with R

To visualize the results of a simulation model of woodland trees within R, I needed an algorithm that could arrange a large number of circles within a rectangle such that no two circles overlapped by more than a specified amount. A colleague had approached this problem earlier by sorting the circles in order of descending size, then randomly dropping each one into the rectangle repeatedly until it landed in a position with acceptable overlap. I suspected a faster and more robust algorithm could be constructed using some kind of "jiggling the circles" approach. Luckily for me, I discovered that Sean McCullough had written a really nice example of circles packing into a cluster using the Processing language. Sean's program is based on an iterative pair-repulsion algorithm in which overlapping circles move away from each other. Based on this, and modifying the algorithm a little, I came up with an R function to produce constrained random layouts of a given set of circles. ...

Circle packing in R (again)

Back in 2010 I posted some R code for circle packing . Now, just five years later, I've ported the code to Rcpp and created a little package which you can find at GitHub . The main function is circleLayout which takes a set of overlapping circles and tries to find a non-overlapping arrangement for them. Here's an example: And here's the code: # Create some random circles, positioned within the central portion # of a bounding square, with smaller circles being more common than # larger ones. ncircles

Homogeneity analysis of hierarchical classifications

I've spent more years than I care to remember analysing vegetation survey data (typically species abundances in plots) using a variety of software including my own algorithms coded in FORTRAN and C++. A recent query on the r-help list, about how to determine the number of groups to define in a hierarchical classification produced with the hclust function, prompted me to unearth one of these algorithms, homogeneity analysis 1 , which can help to visualize how different levels of grouping partition the variability in a distance matrix. This algorithm is extremely simple. The classification is progressively divided into groups, with all groups being defined at the same dendrogram height. At each level of grouping, the average of within-group pairwise distances is calculated. Homogeneity is then defined as: H = 1 - Dav within-group - Dav total where Dav total is the average pairwise distance in the dataset as a whole. For data were there exist well-defined clusters of values, a...

Last Resort Software

Search This Blog