java - 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

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