Cubic Hermite spline interpolation is a very important tool yet on trying to implement it for my latest project I discovered that there's very little concrete information out there. It's taken me four full days trawling through the maths, buggy code examples, technical papers and various youtube videos to figure it out and write a simple Java class.
It's not actually that complicated providing that you can visualize polynomials. I used Desmos to check the maths...
- hermite.png (83.35 KiB) Viewed 5854 times
And here's the scary single line of code that does the magic...
return ( y[ i ] * ( 1 + 2 * t ) + m[ i ] * t ) * ( 1 - t ) * ( 1 - t ) + ( y[ i + 1 ] * ( 3 - 2 * t ) + m[ i + 1 ] * ( t - 1 ) ) * t * t;
Here's the code for the class. Feel free to use it without restriction.
Code: Select all
import java.util.Arrays;
// cubic spline interpolation
public class CubicSpline
{
private double[] m; // array of tangents
private double[] y; // array of y values
// precompute tangents for subsequent calls to interpolate()
// yInputs is an array of at least two Y values that represent the control points (aka nodes)
// the x values are implicit integers equidistant from 0 to the length of yInputs - 1
// if monotone is true the tangents are adjusted to prevent non-monotonic artifacts
// i.e. overshoot that doesn't seem appropriate given the input data
// but in practice it can be safely set to false to speed things up
// note the yInput data does not need to be monotonic
public void precompute( double[] yInputs, boolean monotone )
{
// n is the number of control points
int n = yInputs.length;
// make our own copy to decouple from the outside world...
y = Arrays.copyOf( yInputs, n );
double[] d = new double[ n - 1 ]; // temp array of slopes
m = new double[ n ]; // array of tangents
// calc slopes...
for( int i = 0; i < n - 1; i++ )
d[ i ] = y[ i + 1 ] - y[ i ]; // the change in y
// calc tangets...
m[ 0 ] = d[ 0 ];
for( int i = 1; i < n - 1; i++ )
m[ i ] = ( d[ i - 1 ] + d[ i ] ) * 0.5; // the average of the slopes
m[ n - 1 ] = d[ n - 2 ];
// optionally modify tangents to preserve monotonicity...
if( monotone )
{
for( int i = 0; i < n - 1; i++ )
{
if( d[ i ] == 0 )
{
// zero slope
m[ i ] = 0;
m[ i + 1 ] = 0;
}
else
{
// non-zero slope
double a = m[ i ] / d[ i ];
double b = m[ i + 1 ] / d[ i ];
double h = Math.hypot( a, b );
if( h > 9 )
{
// adjust the tangents...
double t = 3 / h;
m[ i ] = t * a * d[ i ];
m[ i + 1 ] = t * b * d[ i ];
}
}
}
}
}
// precompute() must be called before using this method
// interpolate f( x ) for precomputed values using cubic Hermite spline interpolation
// 0 <= x <= number of y values - 1
// out of bounds x handled safely by clamping
// results pass exactly through every control point
public double interpolate( double x )
{
// handle limits...
if( x <= 0 )
return y[ 0 ];
int maxIndex = y.length - 1;
if( x >= maxIndex )
return y[ maxIndex ];
// array index is integer part of x...
int i = (int) Math.floor( x );
// difference is fractional part of x...
double t = x - i;
// compute the cubic Hermite spline interpolation...
// h00( t ) * y0 +
// h10( t ) * m0 +
// h01( t ) * y1 +
// h11( t ) * m1
// where the basis functions are the following polynomials...
// h00( t ) == 2t^3 - 3t^2 + 1
// h10( t ) == t^3 - 2t^2 + t
// h01( t ) == -2t^3 + 3t^2
// h11( t ) == t^3 - t^2
// these can be rearranged as follows...
// h00( t ) == ( 1 + 2 * t ) * ( 1 - t ) ^ 2 )
// h10( t ) == t * ( 1 - t) ^ 2
// h01( t ) == t ^ 2 * ( 3 - 2 * t )
// h11( t ) == t ^ 2 * ( t - 1 )
// and coded efficiently like so...
return ( y[ i ] * ( 1 + 2 * t ) + m[ i ] * t ) * ( 1 - t ) * ( 1 - t )
+ ( y[ i + 1 ] * ( 3 - 2 * t ) + m[ i + 1 ] * ( t - 1 ) ) * t * t;
}
}
I've tried to document the code so that everything is obvious but to test the CubicSpline class you could do something like this...
Code: Select all
// in Initialize()...
testCubic = new CubicSpline();
testCubic.precompute( testData, true );
// in ProcessSample()...
testX += 0.0005;
if( testX >= testData.length - 1 )
testX = 0;
outputSocket.SetValue( testCubic.interpolate( testX ) );
// data...
private double[] testData = { 1, 1.5, 2, 3, 4, 3.5, 3, 4, 4.5, 4, 3.75, 3.5, 3, 2.75, 3, 2 };
private CubicSpline testCubic;
private double testX;
Edited to fix a typo in one of the comments explaining h11()
Thanks for sharing!