# C# Cubic Spline Interpolation

, 25 Jul 2014 CPOL
 Rate this:
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.

## Share

United States
No Biography provided

 First Prev Next
 Terminology YvesDaoust 2-Apr-15 0:34
 Get value of Y based on X Member 11204817 6-Nov-14 3:33
 Re: Get value of Y based on X Ryan Seghers 7-Nov-14 7:31
 Re: Get value of Y based on X Member 11204817 13-Nov-14 3:11
 Little suggestion Alexandr Stefek 23-Jul-14 21:37
 Re: Little suggestion Ryan Seghers 25-Jul-14 14:10
 Re: Little suggestion Alexandr Stefek 28-Jul-14 22:30
 Evaluate the spline at the specified Y coordinates Menelay 11-Jul-14 5:59
 Re: Evaluate the spline at the specified Y coordinates Ryan Seghers 25-Jul-14 14:16
 Re: Evaluate the spline at the specified Y coordinates Menelay 28-Jul-14 21:32
 My vote of 5 MarkBoreham 30-Jan-14 8:40
 problem when x1 < x0 Member 10494555 29-Jan-14 4:06
 Re: problem when x1 < x0 Ryan Seghers 29-Jan-14 7:41
 Re: problem when x1 < x0 Member 10494555 29-Jan-14 9:05
 Re: problem when x1 < x0 Ryan Seghers 29-Jan-14 13:44
 Re: problem when x1 < x0 [modified] Member 10494555 29-Jan-14 22:25
 Re: problem when x1 < x0 Ryan Seghers 30-Jan-14 11:47
 Re: problem when x1 < x0 LutzL 24-Feb-14 23:37
 Re: problem when x1 < x0 Ryan Seghers 25-Jul-14 14:11
 My vote of 5 Grump 30-Sep-13 8:28
 Re: My vote of 5 Ryan Seghers 30-Sep-13 10:16
 Re: My vote of 5 Grump 30-Sep-13 22:06
 Re: My vote of 5 Ryan Seghers 6-Oct-13 10:14
 My vote of 5 [modified] Douglas Smallish 21-Sep-13 4:24
 My vote of 5 flobadob1975 1-Jul-13 11:27
 My vote of 5 dhatfield 21-Jun-13 4:44
 Re: My vote of 5 Ryan Seghers 21-Jun-13 6: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.
 My vote of 4 sam.hill 11-Mar-13 17:40
 My vote of 5 Kenneth Haugland 11-Mar-13 8:25
 Last Visit: 31-Dec-99 19:00     Last Update: 17-Apr-15 14:00 Refresh 1