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

Fitting an ellipse to point data

Some time ago I wrote an R function to fit an ellipse to point data, using an algorithm developed by Radim Halíř and Jan Flusser 1 in Matlab, and posted it to the r-help list . The implementation was a bit hacky, returning odd results for some data. A couple of days ago, an email arrived from John Minter asking for a pointer to the original code. I replied with a link and mentioned that I'd be interested to know if John made any improvements to the code. About ten minutes later, John emailed again with a much improved version ! Not only is it more reliable, but also more efficient. So with many thanks to John, here is the improved code: fit.ellipse <- function (x, y = NULL) { # from: # http://r.789695.n4.nabble.com/Fitting-a-half-ellipse-curve-tp2719037p2720560.html # # Least squares fitting of an ellipse to point data # using the algorithm described in: # Radim Halir & Jan Flusser. 1998. # Numerically stable direct least squares fitting of ellipses

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 <- 200 limits <- c(-50, 50) inset <- diff(limits) / 3 rmax <- 20 xyr <- data.frame( x = runif(ncircles, min(limits) + inset, max(limits) - inset), y = runif(ncircles, min(limits) + inset, max(limits) - inset), r = rbeta(ncircles, 1, 10) * rmax) # Next, we use the `circleLayout` function to try to find a non-overlapping # arrangement, allowing the circles to occupy any part of the bounding square. # The returned value is a list with elements for

Graph-based circle packing

The previous two posts showed examples of a simple circle packing algorithm using the packcircles package (available from CRAN and GitHub ). The algorithm involved iterative pair-repulsion to jiggle the circles until (hopefully) a non-overlapping arrangement emerged. In this post we'll look an alternative approach. An algorithm to find an arrangement of circles satisfying a prior specification of circle sizes and tangencies was described by Collins and Stephenson in their 2003 paper in Computation Geometry Theory and Applications. A version of their algorithm was implemented in Python by David Eppstein as part of his PADS library (see CirclePack.py ). I've now ported David's version to R/Rcpp and included it in the packcircles package. In the figure below, the graph on the left represents the desired pattern of circle tangencies: e.g. circle 7 should touch all of, and only, circles 1, 2, 6 and 8. Circles 5, 7, 8 and 9 are internal , while the remaining circles are exter

Last Resort Software

Search This Blog