# 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? 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. 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 = 3d * dy / (2d * (dx)) - f1/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;
for (int i = 1; i <= n; i++) {
coefficients = a[i];
coefficients = b[i];
coefficients = c[i];
coefficients = 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! 🙂