JetCracker

Life-time learner's blog

Constrained Cubic Spline Interpolation in Java

Representation of numeric data is a very common problem. On PC we usually use different data visualisation software from Microsoft Excel or OpenOffice Calc to more advanced utilities for charting, such as QtPlot or Origin. These apps are designed for processing of data and you can easily draw a printable scatter plot or a pie chart. But what if you need to draw a plot in your own application?

Example of charts in CardioMood

There are plenty of libraries that can draw different kinds of charts. To draw charts in the web, you might use Google Charts API, or jQuery Flot for example. In my Android app, I use GraphView and ShinobiCharts.

Spline overshoot example

Line charts look nice when they are smooth. The most popular way to make your chart look smooth is to apply a cubic spline interpolation. This type of interpolation is very common and you have it already implemented in many open source mathematic libraries. However a simple cubic spline has a few flaws.

The result function is smooth (it is actually very smooth – twice continuously differentiable), but it tends to overshoot, especially in places where you don’t expect. This is illustrated in the picture above. When we draw a line chart, we usually have some knowledge about the process we are trying to visualise. For instance, the resulting function must be monotonic, or the function value should never be negative. If this information is critical, we must use more advanced interpolation algorithms.

If the function is monotonic, in Android you can just use MonotoneCubicSpline:

import android.util.Spline;

float x[] = ...
float y[[] = ...
Spline f = Spline.createMonotoneCubicSpline(x, y);

// this will evaluate y0 = f(x0)
float y0 = f.interpolate(x0);

If we want to avoid overshoot for an arbitrary function, we should find some other interpolation algorithm. One of the methods is called Constrained Cubic Spline Interpolation. It was proposed by CJC Kruger in his article: http://www.korf.co.uk/spline.pdf

The algorithms is based on a classic cubic spline algorithm. The key step in it is the calculation of the slope (first derivative) at each point. Intuitively, the slope will be between the slopes of the adjacent straight lines (can be a mean value of the two slopes), but it also should approach zero if the slope of either line approaches zero.

Unfortunately, I couldn’t find an implementation of the algorithm in Java. That’s why I decided to implement Constrained Spline for my project.

I use Apache Commons Math library. To make the interpolation compatible with my project, I’ve overridden UnivariateInterpolator class. Here is a full code of my implementation.

import org.apache.commons.math3.analysis.interpolation.UnivariateInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.MathArrays;

/**
 * Constrained Spline Interpolation algorithm (CJC Kruger).
 * Source: http://www.korf.co.uk/spline.pdf
 *
 * Created by Anton Danshin on 24/12/14.
 */
public class ConstrainedSplineInterpolator implements UnivariateInterpolator {

    @Override
    public PolynomialSplineFunction interpolate(double x[], double y[])
            throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException
    {
        if (x.length != y.length) {
            throw new DimensionMismatchException(x.length, y.length);
        }

        if (x.length < 3) {
            throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
                    x.length, 3, true);
        }

        // Number of intervals.  The number of data points is n + 1.
        final int n = x.length - 1;

        MathArrays.checkOrder(x);

        // Differences between knot points
        final double dx[] = new double[n];
        final double dy[] = new double[n];
        for (int i = 0; i < n; i++) {
            dx[i] = x[i + 1] - x[i];
            dy[i] = y[i + 1] - y[i];
        }

        final double f1[] = new double[n+1]; // F'(x[i])
        for (int i=1; i<n; i++) {
            double slope = dy[i-1]*dy[i];
            if (slope > 0d) {
                // doesn't change sign
                f1[i] = 2d / (dx[i]/dy[i] + dx[i-1]/dy[i-1]);
            } else if (slope <= 0d) {
                // changes sign
                f1[i] = 0d;
            }
        }
        f1[0] = 3d * dy[0] / (2d * (dx[0])) - f1[1]/2d;
        f1[n] = 3d * dy[n-1] / (2d * (dx[n-1])) - f1[n-1]/2d;

        // cubic spline coefficients -- a contains constants, b is linear, c quadratic, d is cubic
        final double a[] = new double[n+1];
        final double b[] = new double[n+1];
        final double c[] = new double[n+1];
        final double d[] = new double[n+1];

        for (int i = 1; i <= n; i++) {
            final double f2a = - 2d * (f1[i] + 2 * f1[i-1]) / dx[i-1] + 6d * dy[i-1] / (dx[i-1]*dx[i-1]);
            final double f2b = 2d * (2d * f1[i] + f1[i-1]) / dx[i-1] - 6d * dy[i-1] / (dx[i-1]*dx[i-1]);
            d[i] = (f2b - f2a) / (6d * dx[i-1]);
            c[i] = (x[i] * f2a - x[i-1] * f2b) / (2d * dx[i-1]);
            b[i] = (dy[i-1] -
                    c[i] * (x[i]*x[i] - x[i-1]*x[i-1]) -
                    d[i] * (x[i]*x[i]*x[i] - x[i-1]*x[i-1]*x[i-1])
            ) / dx[i-1];
            a[i] = y[i-1] - b[i]*x[i-1] - c[i]*x[i-1]*x[i-1] - d[i]*x[i-1]*x[i-1]*x[i-1];
        }

        final PolynomialFunction polynomials[] = new PolynomialFunction[n];
        final double coefficients[] = new double[4];
        for (int i = 1; i <= n; i++) {
            coefficients[0] = a[i];
            coefficients[1] = b[i];
            coefficients[2] = c[i];
            coefficients[3] = d[i];
            final double x0 = x[i-1];
            polynomials[i-1] = new PolynomialFunction(coefficients) {
                @Override
                public double value(double x) {
                    // bypass the standard Apache Commons behavior
                    return super.value(x + x0);
                }
            };
        }

        return new PolynomialSplineFunction(x, polynomials);
    }

} 

Feel free to use this in your own code! 🙂

Advertisements

3 responses to “Constrained Cubic Spline Interpolation in Java

  1. Robert Schmid February 19, 2015 at 19:52

    I was just about to do this. You just saved me a couple of days of work. Have you considered submitting this to commons-math?

  2. V July 28, 2015 at 09:39

    Hello, this is a very interesting implementation. Would it be possible to have a sample code of the implementation? More specifically, of the arguments required and the values returned.
    Taking the case of the natural cubic spline interpolation, how is one to fetch the values from the interpolation?

    double[] x = new double[] {1, 5, 7, 10, 12, 15, 18, 20};
    double[] y = new double[] {229.49, 233.36, 235.64, 242.09, 240.29, 259.16, 254.85, 261.27};

    SplineInterpolator interpolator = new SplineInterpolator();
    PolynomialSplineFunction splines = interpolator.interpolate(x, y);

    Thanks in advance.

    • jetcracker September 9, 2015 at 23:39

      The result of interpolation is an instance of PolynomialSplineFunction which is a UnivariateRealFunction. So you can calculate the value within the range of your source data like this:

      splines.value(x)
      

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: