## 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

Splines are a way to fit data points with a set of connecting curves (splines) 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`

.
`CubicSpline`

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:

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:

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:

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));
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*:

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.");
}
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]);
}
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]);
}
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++)
{
int j = GetNextXIndex(x[i]);
float t = (x[i] - xOrig[j]) / (xOrig[j + 1] - xOrig[j]);
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*:

## 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.