Proper implementation of cubic spline interpolation

Question

I was using one of the proposed algorithms out there but the results are very bad.

I implemented the wiki algorithm

in Java (code below). x(0) is points.get(0), y(0) is values[points.get(0)], α is alfa and μ is mi. The rest is the same as in the wiki pseudocode.

    public void createSpline(double[] values, ArrayList<Integer> points){
    a = new double[points.size()+1];

    for (int i=0; i <points.size();i++)
    {
        a[i] = values[points.get(i)];

    }

    b = new double[points.size()];
    d = new double[points.size()];
    h = new double[points.size()];

    for (int i=0; i<points.size()-1; i++){
        h[i] = points.get(i+1) - points.get(i);
    }

    alfa = new double[points.size()];

    for (int i=1; i <points.size()-1; i++){
        alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
    }

    c = new double[points.size()+1];
    l = new double[points.size()+1];
    mi = new double[points.size()+1];
    z = new double[points.size()+1];

    l[0] = 1; mi[0] = z[0] = 0;

    for (int i =1; i<points.size()-1;i++){
        l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
        mi[i] = h[i]/l[i];
        z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
    }

    l[points.size()] = 1;
    z[points.size()] = c[points.size()] = 0;

    for (int j=points.size()-1; j >0; j--)
    {
        c[j] = z[j] - mi[j]*c[j-1];
        b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
        d[j] = (c[j+1]-c[j])/((double)3*h[j]);
    }


    for (int i=0; i<points.size()-1;i++){
        for (int j = points.get(i); j<points.get(i+1);j++){
            //                fk[j] = values[points.get(i)];
            functionResult[j] = a[i] + b[i] * (j - points.get(i)) 
                                + c[i] * Math.pow((j - points.get(i)),2)
                                + d[i] * Math.pow((j - points.get(i)),3);
        }
    }

}

The result that I get is the following:

but it should be similar to this:

---

I'm also trying to implement the algorithm in another way according to: http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

At first they show how to do linear spline and it's pretty easy. I create functions that calculate A and B coefficients. Then they extend linear spline by adding second derivative. C and D coefficients are easy to calculate too.

But the problems starts when I attempt to calculate the second derivative. I do not understand how they calculate them.

So I implemented only linear interpolation. The result is:

Does anyone know how to fix the first algoritm or explain me how to calculate the second derivative in the second algorithm?

Wallaby

Try it here: codereview.stackexchange.com – Fulcrum

@Fulcrum they say me that they reviewing code quality not code results. – Wallaby

@Fulcrum Because the results here are not the desired results, this question is not a good fit for Code Review. – Oliva

will be nice to see the control points in your charts to see what is really wrong and what not – Marismarisa

One should perhaps not use outdated wikipedia articles that are flagged as confusing. The recent article on spline interpolation is, surprisingly, at spline interpolation. – Forefather

Answer 1

Cubic b-spline has been recently descried in a series of papers by Unser, Thévenaz et al., see among others

M. Unser, A. Aldroubi, M. Eden, "Fast B-Spline Transforms for Continuous Image Representation and Interpolation", IEEE Trans. Pattern Anal. Machine Intell., vol. 13, n. 3, pp. 277-285, Mar. 1991.

M. Unser, "Splines, a Perfect Fit for Signal and Image Processing", IEEE Signal Proc. Mag., pp. 22- 38, Nov. 1999.

and

P. Thévenaz, T. Blu, M. Unser, "Interpolation Revisited," IEEE Trans. on Medical Imaging, vol. 19, no. 7, pp. 739-758, July 2000.

Here are some guidelines.

What are splines?

Splines are piecewise polynomials that are smoothly connected together. For a spline of degree n, each segment is a polynomial of degree n. The pieces are connected so that the spline is continuous up to its derivative of degree n-1 at the knots, namely, the joining points of the polynomial pieces.

How can splines be constructed?

The zero-th order spline is the following

All the other splines can be constructed as

where the convolution is taken n-1 times.

Cubic splines

The most popular splines are cubic splines, whose expression is

Spline interpolation problem

Given a function f(x) sampled at the discrete integer points k, the spline interpolation problem is to determine an approximation s(x) to f(x) expressed in the following way

where the ck's are interpolation coefficients and s(k) = f(k).

Spline prefiltering

Unfortunately, starting from n=3 on,

so that the ck's are not the interpolation coefficients. They can be determined by solving the linear system of equations obtained by forcing s(k) = f(k), namely,

Such an equation can be recast in a convolution form and solved in the transformed z-space as

where

Accordingly,

Proceeding this way is always preferable than affording the solution of a linear system of equations by, for example, LU decomposition.

The solution to the above equation can be determined by noticing that

where

The first fraction is representative of a causal filter, while the second one is representative of an anticausal filter. Both of them are illustrated in the figures below.

Causal filter

Anticausal filter

In the last figure,

The output of the filters can be expressed by the following recursive equations

The above equations can be solved by first determining "initial conditions" for c- and c+. On assuming a periodic, mirrored input sequence fk such that

then it can be shown that the initial condition for c+ can be expressed as

while the initial condition for c+ can be expressed as

Flense

Answer 2

Sorry but Your source code is really a unreadable mess to me so I stick to theory. Here are some hints:

1. SPLINE cubics

   SPLINE is not interpolation but approximation to use them you do not need any derivation. If you have ten points: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9 then cubic spline starts/ends with triple point. If you create function to 'draw' SPLINE cubic curve patch then to assure continuity the call sequence will be like this:

       spline(p0,p0,p0,p1);
       spline(p0,p0,p1,p2);
       spline(p0,p1,p2,p3);
       spline(p1,p2,p3,p4);
       spline(p2,p3,p4,p5);
       spline(p3,p4,p5,p6);
       spline(p4,p5,p6,p7);
       spline(p5,p6,p7,p8);
       spline(p6,p7,p8,p9);
       spline(p7,p8,p9,p9);
       spline(p8,p9,p9,p9);
   

   do not forget that SPLINE curve for p0,p1,p2,p3 draw only curve 'between' p1,p2 !!!

2. BEZIER cubics

   4-point BEZIER coefficients can be computed like this:

       a0=                              (    p0);
       a1=                    (3.0*p1)-(3.0*p0);
       a2=          (3.0*p2)-(6.0*p1)+(3.0*p0);
       a3=(    p3)-(3.0*p2)+(3.0*p1)-(    p0);
   

3. Interpolation

   to use interpolation you must use interpolation polynomials. There are many out there but I prefer to use cubics ... similar to BEZIER/SPLINE/NURBS... so

   • p(t) = a0+a1*t+a2*t*t+a3*t*t*t where t = <0,1>

   The only thing left to do is compute a0,a1,a2,a3. You have 2 equations (p(t) and its derivation by t) and 4 points from the data set. You also must ensure the continuity ... So first derivation for join points must be the same for neighboring curves (t=0,t=1). This leads to system of 4 linear equations (use t=0 and t=1). If you derive it it will create an simple equation depended only on input point coordinates:

       double  d1,d2;
       d1=0.5*(p2-p0);
       d2=0.5*(p3-p1);
       a0=p1;
       a1=d1;
       a2=(3.0*(p2-p1))-(2.0*d1)-d2;
       a3=d1+d2+(2.0*(-p2+p1));
   

   the call sequence is the same as for SPLINE

[Notes]

1. the difference between interpolation and approximation is that:

   interpolation curve goes always through the control points but high order polynomials tend to oscillate and approximation just approaches to control points (in some cases can cross them but usually not).

2. variables:

   • a0,... p0,... are vectors (number of dimensions must match the input points)
   • t is scalar

3. to draw cubic from coefficients a0,..a3

   just do something like this:

       MoveTo(a0.x,a0.y);
       for (t=0.0;t<=1.0;t+=0.01)
        {
        tt=t*t;
        ttt=tt*t;
        p=a0+(a1*t)+(a2*tt)+(a3*ttt);
        LineTo(p.x,p.y);
        }
   

Marismarisa

One can do interpolating cubic splines, where function values, derivatives and second derivatives match at the interpolation points. One just has to solve a large, but banded, system of linear equations. – Forefather

why to heck do it so backwards, cubic curve has c1 connectivity limit (value and first derivation match) if higher connectivity is needed then use more points 5 points -> c2 , 6 points -> c3 and so on ... with 4 points to achieve connectivity c2 then it is not SPLINE anymore ... and not even 4 point cubic because you solving near points also ... Of course I can miss something like new equations discovered in the recent past but i strongly doubt it – Marismarisa

I believe you missed some very old formulas. Yes, given values and derivatives, you can build a piecewise cubic function. If you do not have the derivatives given, then there is much freedom in choosing them. This freedom is reduced by demanding minimized second derivatives. And indeed, this leads to a twice differentiable piecewise cubic with linear ends, which is uniquely determined by a tridiagonal linear system. See spline interpolation – Forefather

And usually, the term "spline" is reserved for those piecewise cubic functions that have minimal curvature (or second derivative) for the given constraints. With cubic splines one gets no better than C2, the third derivative is almost always piecewise constant with jumps. If there are no jumps, then the third derivative is constant and the function is a cubic polynomial. – Forefather

Answer 3

See spline interpolation, although they give only a usable 3x3 example. For more sample points, say N+1 enumerated x[0..N] with values y[0..N] one has to solve the following system for the unknown k[0..N]

              2*    c[0]    * k[0] + c[0] * k[1] == 3*   d[0];
c[0] * k[0] + 2*(c[0]+c[1]) * k[1] + c[1] * k[2] == 3*(d[0]+d[1]);
c[1] * k[1] + 2*(c[1]+c[2]) * k[2] + c[2] * k[3] == 3*(d[1]+d[2]);
c[2] * k[2] + 2*(c[2]+c[3]) * k[3] + c[3] * k[4] == 3*(d[2]+d[3]);
...
c[N-2]*k[N-2]+2*(c[N-2]+c[N-1])*k[N-1]+c[N-1]*k[N] == 3*(d[N-2]+d[N-1]);
c[N-2]*k[N-1] + 2*c[N-1]        *  k[N]          == 3*   d[N-1];

where

c[k]=1/(x[k+1]-x[k]); d[k]=(y[k+1]-y[k])*c[k]*c[k];

This can be solved using the Gauß-Seidel iteration (was this not invented exactly for this system?) or your favorite Krylov space solver.

Forefather

Answer 4

// File: cbspline-test.cpp

// C++ Cubic Spline Interpolation using the Boost library.
// Interpolation is an estimation of a value within two known values in a sequence of values. 

// From a selected test function y(x) = (5.0/(x))*sin(5*x)+(x-6), 21 numbers of (x,y) equal-spaced sequential points were generated. 
// Using the known 21 (x,y) points, 1001 numbers of (x,y) equal-spaced cubic spline interpolated points were calculated.
// Since the test function y(x) is known exactly, the exact y-value for any x-value can be calculated.
// For each point, the interpolation error is calculated as the difference between the exact (x,y) value and the interpolated (x,y) value. 
// This program writes outputs as tab delimited results to both screen and named text file

// COMPILATION: $ g++ -lgmp -lm -std=c++11 -o cbspline-test.xx cbspline-test.cpp
// EXECUTION:   $ ./cbspline-test.xx 

// GNUPLOT 1:   gnuplot> load "cbspline-versus-exact-function-plots.gpl"
// GNUPLOT 2:   gnuplot> load "cbspline-interpolated-point-errors.gpl"

#include <iostream>
#include <fstream>
#include <boost/math/interpolators/cubic_b_spline.hpp>
using namespace std;

// ========================================================
int main(int argc, char* argv[]) {
// ========================================================

    // Vector size 21 = 20 space intervals, size 1001 = 1000 space intervals  
    std::vector<double> x_knot(21), x_exact(1001), x_cbspline(1001);
    std::vector<double> y_knot(21), y_exact(1001), y_cbspline(1001); 
    std::vector<double> y_diff(1001);

    double x;   double x0 = 0.0;    double t0 = 0.0;    
    double xstep1 = 0.5;    double xstep2 = 0.01;       

    ofstream outfile;                                   // Output data file tab delimited values
    outfile.open("cbspline-errors-1000-points.dat");    // y_diff = (y_exact - y_cbspline)

    // Handling zero-point infinity (1/x) when x = 0
    x0 = 1e-18; 
    x_knot[0] = x0; 
    y_knot[0] = (5.0/(x0))*sin(5*x0)+(x0-6);            // Selected test function

    for (int i = 1; i <= 20; i++) { // Note: Loop i begins with 1 not 0, 20 points 
        x = (i*xstep1); 
        x_knot[i] = x;
        y_knot[i] = (5.0/(x))*sin(5*x)+(x-6); 
    }

    // Execution of the cubic spline interpolation from Boost library
    // NOTE: Using xstep1 = 0.5 based on knot points stepsize (for 21 known points) 
    boost::math::cubic_b_spline<double> spline(y_knot.begin(), y_knot.end(), t0, xstep1);

    // Again handling zero-point infinity (1/x) when x = 0
    x_cbspline[0] = x_knot[0];  
    y_cbspline[0] = y_knot[0];

    for (int i = 1; i <= 1000; ++i) {       // 1000 points.
        x = (i*xstep2);  
        x_cbspline[i] = x;
        // NOTE: y = spline(x) is valid for any value of x in xrange [0:10] 
        // meaning, any x within range [x_knot.begin() and x_knot.end()]
        y_cbspline[i] = spline(x);   
    }

    // Write headers for screen display and output file
    std::cout << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  
    outfile   << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;  

    // Again handling zero-point infinity (1/x) when x = 0
    x_exact[0] = x_knot[0]; 
    y_exact[0] = y_knot[0];
     y_diff[0] = (y_exact[0] - y_cbspline[0]);
    std::cout << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to screen
    outfile   << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl;  // Write to text file

    for (int i = 1; i <= 1000; ++i) {       // 1000 points
        x = (i*xstep2);  
        x_exact[i] = x;
        y_exact[i] = (5.0/(x))*sin(5*x)+(x-6); 
         y_diff[i] = (y_exact[i] - y_cbspline[i]);
        std::cout << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to screen
        outfile   << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl;  // Write to text file
    }
    outfile.close();

return(0);
}
// ========================================================
/*
# GNUPLOT script 1
# File: cbspline-versus-exact-function-plots.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-12.0:20.0]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 2.0
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y(x)"
set title "Function points y(x) = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:2 with lines linetype 3 linewidth 1.0 linecolor 'blue' title 'exact function curve', 'cbspline-errors-1000-points.dat' using 1:3 with lines linecolor 'red' linetype 3 linewidth 1.0 title 'cubic spline interpolated curve'

*/
// ========================================================
/*
# GNUPLOT script 2
# File: cbspline-interpolated-point-errors.gpl

set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-2.50:2.5]
# set autoscale         # set xtics automatically
set xtics 0.5
set ytics 0.5
# set xtic auto         # set xtics automatically
# set ytic auto         # set ytics automatically
set grid x
set grid y
set xlabel "x" 
set ylabel "y"
set title "Function points y = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid

plot 'cbspline-errors-1000-points.dat' using 1:4 with lines linetype 3 linewidth 1.0 linecolor 'red' title 'cubic spline interpolated errors'

*/
// ========================================================

Spinoza