C# Cubic Spline Interpolation

• Download source

Introduction

This is an implementation of cubic spline interpolation based on the Wikipedia articles Spline Interpolation and Tridiagonal Matrix Algorithm. My goal in creating this was to provide a simple, clear implementation that matches the formulas in the Wikipedia articles closely, rather than an optimized implementation.

Background

Spline algorithms are a way to fit data points with a set of connecting curves (each one is called a spline) such that values between data points can be computed (interpolated). There are various types/orders of equations that can be used to specify the splines including linear, quadratic, cubic, etc.

A set of N points requires N-1 splines to connect them. Each spline is described by an equation (e.g., a polynomial). The coefficients in those polynomials are initially unknown, and the spline algorithm computes them. Once the coefficients are known, to compute the Y for any particular X you find the spline that is relevant to that X and use those coefficients to evaluate the polynomial at that X.

To compute the splines' coefficients, the algorithm uses the standard linear algebra procedure of creating a system of equations and then solving that system. The equations embody the constraints on the splines which are that they intersect the points to be fitted, and that at each junction between splines the prior splines' first and second derivatives equal the successor splines' first and second derivatives. That doesn't provide enough constraints yet to fully determine the system so two more are added that specify that the second derivatives of the first and last points are zero. Those final constraints are called the "natural" spline style. Alternatively, user Quience (below) has provided code such that the caller can specify the slopes at the first and last points explicitly. With either set of additional constraints, the system becomes solvable.

Fortunately, the system of equations is constrained to be tridiagonal, so that solving the system is O(N) using the Thomas algorithm, rather than the O(N^3) of Gaussian elimination for example. Here is the tridiagonal matrix corresponding to the above image which is computed to solve the Wikipedia article's example problem:

2.0000, 1.0000, 0.0000
1.0000, 2.6667, 0.3333
0.0000, 0.3333, 0.6667

This example is a little too small to illustrate a tridiagonal matrix well, but tridiagonal means that only the main diagonal (row == col) elements, the diagonal below the main (row == col + 1), and the diagonal above the main (row == col - 1) can have non-zero values.

The linked Wikipedia articles do a good job of explaining splines so I won't provide any more redundant explanation here.

Using the Code

I followed the Wikipedia articles closely for the implementation, including variable names where possible. Some comments refer to equation numbers in the articles (which unfortunately will break at some point when the articles change).

The primary classes are CubicSpline and TriDiagonalMatrix. The CubicSpline class uses TriDiagonalMatrix to solve the system of equations. Generally, if you are just using this code to fit curves, you would just use CubicSpline.

In the source, there is also Program.cs which has a Main() that runs several functions that exercise the classes.

Here is an example use of CubicSpline from the TestSplineOnWikipediaExample() method in Program.cs. This part of the code sets up the test data to be fitted:

// Create the test points.
float[] x = new float[] { -1.0f, 0.0f, 3.0f };
float[] y = new float[] { 0.5f, 0.0f, 3.0f };

Console.WriteLine("x: {0}", ArrayUtil.ToString(x));
Console.WriteLine("y: {0}", ArrayUtil.ToString(y));

The next step is to create a float[] (here called xs) which contains the values of X for which you want to determine the Y values. In other words, these are the X values for which you want to plot points on the graph:

// Create the upsampled X values to interpolate
int n = 20;
float[] xs = new float[n];
float stepSize = (x[x.Length - 1] - x[0]) / (n - 1);

for (int i = 0; i < n; i++)
{
    xs[i] = x[0] + i * stepSize;
}

And next is to use CubicSpline to fit the data and then evaluate the fitted spline at the xs values generated above:

// Solve
CubicSpline spline = new CubicSpline();
float[] ys = spline.FitAndEval(x, y, xs);

Console.WriteLine("xs: {0}", ArrayUtil.ToString(xs));
Console.WriteLine("ys: {0}", ArrayUtil.ToString(ys));

// Plot
string path = @"..\..\spline-wikipedia.png";
PlotSplineSolution("Cubic Spline Interpolation - Wikipedia Example", x, y, xs, ys, path);

This displays the following console output:

x: -1, 0, 3
y: 0.5, 0, 3
xs: -1, -0.7894737, -0.5789474, -0.368421, -0.1578947, 0.05263158, ...
ys: 0.5, 0.3570127, 0.2245225, 0.1130267, 0.0330223, -0.005029888, ...

The PlotSplineSolution() method uses the .NET Framework's System.Windows.Forms.DataVisualization.Charting classes to plot a graph of the input points and the fitted and interpolated spline curve. The image created by this is the one included in this article (above).

Implementation

The core of the spline fitting function sets up the tridiagonal matrix and then uses it to solve the system of equations. The tridiagonal matrix is not represented as a matrix but rather three 1-d arrays, A, B, and C. Array A is the sub-diagonal, B is the diagonal, and C is the super-diagonal, to match the Wikipedia article names. Once the tridiagonal matrix is set up, the spline fitting function calls TriDiagonalMatrix.Solve().

Here is the tridiagonal matrix solver function from TriDiagonalMatrix.cs:

// Solve the system of equations this*x=d given the specified d.
public float[] Solve(float[] d)
{
    int n = this.N;

    if (d.Length != n)
    {
        throw new ArgumentException("The input d is not the same size as this matrix.");
    }

    // cPrime
    float[] cPrime = new float[n];
    cPrime[0] = C[0] / B[0];

    for (int i = 1; i < n; i++)
    {
        cPrime[i] = C[i] / (B[i] - cPrime[i-1] * A[i]);
    }

    // dPrime
    float[] dPrime = new float[n];
    dPrime[0] = d[0] / B[0];

    for (int i = 1; i < n; i++)
    {
        dPrime[i] = (d[i] - dPrime[i-1]*A[i]) / (B[i] - cPrime[i - 1] * A[i]);
    }

    // Back substitution
    float[] x = new float[n];
    x[n - 1] = dPrime[n - 1];

    for (int i = n-2; i >= 0; i--)
    {
        x[i] = dPrime[i] - cPrime[i] * x[i + 1];
    }

    return x;
}

The following is the core of the CubicSpline.Eval() method. This method takes a float[] x which contains the x values you want to compute y for using the fitted splines. The loop is stepping through each value of x computing the corresponding value of y.

for (int i = 0; i < n; i++)
{
    // Find which spline can be used to compute this x
    int j = GetNextXIndex(x[i]);

    // Evaluate using j'th spline
    float t = (x[i] - xOrig[j]) / (xOrig[j + 1] - xOrig[j]);
    // equation 9 in the Wiki
    y[i] = (1 - t) * yOrig[j] + t * yOrig[j + 1] + t * (1 - t) * (a[j] * (1 - t) + b[j] * t); 
}

Another Example

Here is another chart. This one was created by the TestSpline() method in Program.cs:

Computing Slope

In addition to computing the Y for each X, one can also compute the slope of the spline. It is specified as q' in the Wikipedia article and the formula is equation 5. The method to evaluate slope is CubicSpline.EvalSlope() and it returns a float array (one slope value for each X you provide). You must have called either Fit() (or FitEndEval()) before calling this. Here is an example use:

// Try slope (spline is already computed at this point, see above code example)
float[] slope = spline.EvalSlope(xs); // Same xs as first example above
string slopePath = @"..\..\spline-wikipedia-slope.png";
PlotSplineSolution("Cubic Spline Interpolation - Wikipedia Example - Slope",
                   x, y, xs, ys, slopePath, slope);

Here is the resulting chart:

Input Constraints and Parametric Fitting

The normal cubic spline algorithm works on 2-d points where y is a function of x, i.e. y=f(x), and y has a single value for each x. However, user LutzL in the comments below has pointed out a clever way to use splines to fit sequences of points that do not fit this definition:

[You] invent a third time coordinate that increases monotonically, T=0,1,2,3 or start with T=0 and increment by the distance from the current point to the next point. Use then the (T,X) and (T,Y) pairs to compute two cubic splines x(t) and y(t) and draw the curve (x(t),y(t)) as result.

I implemented this as a new static method FitParametric() in the CubicSpline class. Here is an example use:

float[] x = { 0.5f, 2.0f, 3.0f, 4.5f, 3.0f, 2.0f };
float[] y = { 4.0f, 2.0f, 6.0f, 4.0f, 3.0f, 5.0f };
float[] xs, ys;
CubicSpline.FitParametric(x, y, 100, out xs, out ys);

Here is the resulting chart:

Thanks to user YvesDaoust for the correct term for this method.

Points of Interest

Surprisingly, this simple implementation is more than twice as fast as a C# port of the very terse Numerical Recipes in C++ (Press et al, Cambridge University Press, 2002) version in the scenario I benchmarked. I only benchmarked one scenario and did not investigate the results much so shouldn't draw too many conclusions, but it is interesting anyway. I benchmarked the implementations using a Release build by Visual Studio 2012, running on a Win 7 Core i5 2500k. I used 10,000 random points and interpolated 100,000 points, with 100 repetitions. See the Program.cs method TestPerf(). I did not include the Numerical Recipes in C++ version in the code due to licensing and copyright issues. [Update: This performance difference is mostly likely due to the fact that the NR in C implementation does a search for the appropriate spline for each evaluation, whereas I require that the points to be evaluated are sorted and therefore I can do a simultaneous traverse.]

The Chart class in the System.Windows.Forms.DataVisualization.Charting namespace worked well for my simple charts. The classes in there make it very easy to render a chart and save to a file.

It is unfortunate that C#/.NET generics (still) cannot handle creating a generic type constraint that allows you to handle both float and double with the same code. I implemented this with floats because I typically start with floats until it's proven that I need doubles.

History

• March 10, 2013 - First version
• September 20, 2013 - Added description of first and last point constraints including explicit slope arguments, and misc other minor changes
• July 23, 2014 - Updated code to add slope computation and added the Computing Slope section to the article
• July 25, 2014 - Added Input Constraints and Geometric Fitting section, and updated code
• April 1, 2016 - Renamed 'geometric' to 'parametric' and added start and end slope specification to parametric fitting