Interpolation with Cubic Splines

Adobe Photoshop, Gimp and many other image editors have a tool called Curves. There are several "knots" and while you move them, it computes a curve in 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 a 2D space.

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

I recommend you to read about Splines on Wikipedia..

The algorithm for cubic interpolation is also described on Wikipedia, so I just copied it.

Implementation

In acubic splines interpolation, the input is the set of knots + first derivatives for each knot. I decided to represent it with three arrays: an array of X values (xs), an array of Y values (ys) and an array of derivative values (ks). Here is the 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 a 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, the first and the second derivative of a cubic spline should be continuous. Below is the function, which generates the array of "ks", which have that property. It also sets the second derivative of the first and the last point to zero (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, which you can find here: CSPL - small Javascript library for cubic spline interpolation.