Zen and the art of cublic spline interpolation

ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Zen and the art of cublic spline interpolation

Post by ColinP »

Here's a little gift.

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
hermite.png (83.35 KiB) Viewed 2447 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;
	}
	
}

The forum software has messed up the tabs and it's written in my own lexical style rather than K&R (which I detest) but you shouldn't have any trouble reading it.

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()
Last edited by ColinP on Fri Jun 16, 2023 2:15 pm, edited 1 time in total.
UrbanCyborg
Posts: 588
Joined: Mon Nov 15, 2021 9:23 pm

Re: Zen and the art of cublic spline interpolation

Post by UrbanCyborg »

As Darth Vader said, "Impressive. Most impressive." ;) Thanks for sharing!

Reid
Cyberwerks Heavy Industries -- viewforum.php?f=76
User avatar
seal58
Posts: 351
Joined: Fri Jul 12, 2019 5:28 pm
Location: Rostock, Germany
Contact:

Re: Zen and the art of cublic spline interpolation

Post by seal58 »

Hello Colin,

finally with the graphics now it is clear to me what spline interpolation is meant to be. Thanks for doing this work and sharing it.

Roland
UrbanCyborg
Posts: 588
Joined: Mon Nov 15, 2021 9:23 pm

Re: Zen and the art of cublic spline interpolation

Post by UrbanCyborg »

Just rummaged through my ancient copy of "Computer Graphics: Principles and Practice, 2nd Ed." by Foley, van Dam, et al. ISBN 0-201-12110-7, and it's got lots on bicubics in general, and Bézier and Hermite patches in particular. One of the original indispensible texts.

Reid
Cyberwerks Heavy Industries -- viewforum.php?f=76
ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Re: Zen and the art of cublic spline interpolation

Post by ColinP »

seal58 wrote: Tue Jun 06, 2023 6:10 am finally with the graphics now it is clear to me what spline interpolation is meant to be. Thanks for doing this work and sharing it.
Yes, it makes far more sense when you see the basis functions plotted. Especially as h01() (the green line) is exactly the same -2x^3+3x^2 polynomial that I'm always raving about and h00() is the same thing mirrored about x = 0.5. It's all really elegant.

Visualization of maths is so important. To this end I recommend everyone looks at Desmos if they haven't already.
UrbanCyborg wrote: Tue Jun 06, 2023 1:32 pm Just rummaged through my ancient copy of "Computer Graphics: Principles and Practice, 2nd Ed." by Foley, van Dam, et al. ISBN 0-201-12110-7, and it's got lots on bicubics in general, and Bézier and Hermite patches in particular. One of the original indispensible texts.

Reid
Yeah, that book used to sit next to Muscial Applications of Microprocessors by Hal Chamberlin on my shelves back in the old days. They're both in storage now but perhaps I should dig them out and read them again.
Steve W
Posts: 758
Joined: Thu Jul 16, 2020 5:55 pm

Re: Zen and the art of cublic spline interpolation

Post by Steve W »

ColinP wrote: Tue Jun 06, 2023 2:17 pm Yeah, that book used to sit next to Muscial Applications of Microprocessors by Hal Chamberlin on my shelves back in the old days. They're both in storage now but perhaps I should dig them out and read them again.
Not sure about the other book, but why risk inhaling the dust when there's this: http://sites.music.columbia.edu/cmc/cou ... berlin.pdf
ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Re: Zen and the art of cublic spline interpolation

Post by ColinP »

Reid, you've mentioned Bezier a few times.

I absolutely love the recursive lerping construction - it's one of those concepts that strike me as intrinsically beautiful and the convex hull property is useful but because the results don't pass through all the control points I couldn't see how to use it for interpolation.

So I'm wondering if I've missed something and there's some efficient way to derive the Bezier control points from the data? I'm sure it can be done but a quick glance brought up only very expensive solutions,
ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Re: Zen and the art of cublic spline interpolation

Post by ColinP »

Cheers Steve. I actually came across that PDF after posting that comment. A quick scan reinforced my opinion that it's an excellent book although the assembly language examples look somewhat quaint now!

I do like physical books and have a good library in storage but it's a pain to access.
ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Re: Zen and the art of cublic spline interpolation

Post by ColinP »

Just had a quick look again at the Bezier thing and found this..

https://stackoverflow.com/questions/192 ... ong-that-c

I estimate one could rewrite that to cost 8 multiplications to derive the two intermediate control points but it would be in the precompute phase anyway so not a really worry. At some point in the future I'll have to look again at using Bezier instead of Hermite but for now the code I'm already using is fast enough for my purposes.
ColinP
Posts: 939
Joined: Mon Aug 03, 2020 7:46 pm

Re: Zen and the art of cublic spline interpolation

Post by ColinP »

The other book that was so good it was known just as Foley and Van Dam (rather like K&R was known as K&R or even just the "White Book") - appears to have been rewritten as a third edition and is available as a PDF if you google it but I'm not sure if it's OK from a copyright POV as I've never heard of the Atlantic International University. It seems a bit strange that you can download a first class 1,263 page book for free so I'm not going to post a link.
Post Reply

Return to “Module Designer”