JetCracker

Life-time learner's blog

Tag Archives: Constrained Spline

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! 🙂