Click here to Skip to main content
Click here to Skip to main content
Go to top

C# Cubic Spline Interpolation

, 25 Jul 2014
Rate this:
Please Sign up or sign in to vote.
This article presents a from-scratch implementation of cubic spline interpolation in C#.

 

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 Geometric 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 FitGeometric() 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.FitGeometric(x, y, 100, out xs, out ys);

Here is the resulting chart:

I called it "Geometric" for lack of a better term, I'd be happy to hear what the correct mathematical terms are for both categories of points (y=f(x) with single y for each x, versus arbitrary x,y pairs).

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 - Add description of first and last point constraints including explicit slope arguments, and misc other minor changes. 

July 23, 2014 - Update code to add slope computation and add the Computing Slope section to the article.

July 25, 2014 - Add Input Constraints and Geometric Fitting section, and update code.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Ryan Seghers

United States United States
No Biography provided

Comments and Discussions

 
GeneralRe: My vote of 5 PinmemberGrump30-Sep-13 21:06 
GeneralRe: My vote of 5 PinmemberRyan Seghers6-Oct-13 9:14 
GeneralMy vote of 5 [modified] PinmemberDouglas Smallish21-Sep-13 3:24 
QuestionSmall change to enable constraining the slope at either end of the dataset PinmemberQuience13-Sep-13 5:54 
AnswerRe: Small change to enable constraining the slope at either end of the dataset PinmemberRyan Seghers20-Sep-13 18:06 
GeneralMy vote of 5 Pinmemberflobadob19751-Jul-13 10:27 
GeneralMy vote of 5 Pinmemberdhatfield21-Jun-13 3:44 
GeneralRe: My vote of 5 PinmemberRyan Seghers21-Jun-13 5:31 
Other than using the same general concept of a spline, they are quite different. This implementation is 1-d and Phillipe's is 2-d, but more importantly Phillipe uses B-splines whereas this article is about cubic splines. B-splines and cubic splines are very different.
GeneralMy vote of 4 Pinmembersam.hill11-Mar-13 16:40 
GeneralMy vote of 5 PinmemberKenneth Haugland11-Mar-13 7:25 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web03 | 2.8.140916.1 | Last Updated 25 Jul 2014
Article Copyright 2013 by Ryan Seghers
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid