Interpolation with Cubic Splines

Adobe Photoshop, Gimp and many other graphic programs have a tool called Curves. You have several "knots" and while you move them, it computes a curve between. I always wondered, how is this curve computed.

About cubic splines

First, I must say, that we are going to interpolate a function of one real variable, not a curve in 2D space.

Having several ($n$) "knots" - points in 2D space - pairs $(x, f(x))$, we would like to interpolate between them with smooth curve. We can interpolate between them with polynomial of degree $n-1$, but it usually generates too "wavy" curve. The idea with splines is, that we compute unique polynomial for each sub-interval between knots, and these polynomials should look good together - they should have the same first, second, third ... derivatives in knots.

I recommend you to read about Splines on Wikipedia..

Algorithm for cubic interpolation is described also no Wikipedia, so I just copied it.

Implementation

In cubic splines interpolation, the input is knots + first derivatives for each knot. I decided to represent it with 3 arrays: array of X values (xs), array of Y values (ys) and array of derivative values (ks). Here is my function for evaluating a cubic spline for any point X:

CSPL.evalSpline = function(x, xs, ys, ks)
{
    var i = 1;
    while(xs[i]<x) i++;
		
    var t = (x - xs[i-1]) / (xs[i] - xs[i-1]);
		
    var a =  ks[i-1]*(xs[i]-xs[i-1]) - (ys[i]-ys[i-1]);
    var b = -ks[i  ]*(xs[i]-xs[i-1]) + (ys[i]-ys[i-1]);
		
    var q = (1-t)*ys[i-1] + t*ys[i] + t*(1-t)*(a*(1-t)+b*t);
    return q;
}

In code above, we can give the function any first derivatives (ks). But we usually don't have any specific derivatives, we just want the curve to be as smooth as possible.

From the spline definition, first and second derivative of cubic spline should be continuous. Below is the function, which generates the "ks" array, which have that property. It also sets the second derivative of first and last points to 0 (Natural Spline).

CSPL.getNaturalKs = function(xs, ys, ks)    // in x values, in y values, out k values
{
    var n = xs.length-1;
    var A = CSPL._gaussJ.zerosMat(n+1, n+2);
        
    for(var i=1; i<n; i++)    // rows
    {
        A[i][i-1] = 1/(xs[i] - xs[i-1]);
        
        A[i][i  ] = 2 * (1/(xs[i] - xs[i-1]) + 1/(xs[i+1] - xs[i])) ;
        
        A[i][i+1] = 1/(xs[i+1] - xs[i]);
        
        A[i][n+1] = 3*   ( (ys[i]-ys[i-1])/ ((xs[i] - xs[i-1])*(xs[i] - xs[i-1])) 
                         +  (ys[i+1]-ys[i])/ ((xs[i+1] - xs[i])*(xs[i+1] - xs[i])) );
    }
    
    A[0][0  ] = 2/(xs[1] - xs[0]);
    A[0][1  ] = 1/(xs[1] - xs[0]);
    A[0][n+1] = 3 * (ys[1] - ys[0]) / ((xs[1]-xs[0])*(xs[1]-xs[0]));
    
    A[n][n-1] = 1/(xs[n] - xs[n-1]);
    A[n][n  ] = 2/(xs[n] - xs[n-1]);
    A[n][n+1] = 3 * (ys[n] - ys[n-1]) / ((xs[n]-xs[n-1])*(xs[n]-xs[n-1]));
        
    CSPL._gaussJ.solve(A, ks);        
}

Where A is $n \times (n+1)$ matrix and CSPL._gaussJ.solve solves the system of linear equations.

I put all my spline code into a small Javascript library named CSPL, here it is - CSPL - small Javascript library for cubic spline interpolation.

Old comments (closed because of spam)

One Comment

  1. Draw a cubic spline curve through Google Earth Placemarks | netkingcol said:

    [...] cubic spline functions provided by Ivan Kuckir at the Faculty of Mathematics and Physics at Charles University in Prague. [...]

    January 2nd, 2014