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?

Answer 1

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);
        }