Click here to Skip to main content
15,867,328 members
Articles / Programming Languages / C#

C# Cubic Spline Interpolation

Rate me:
Please Sign up or sign in to vote.
4.92/5 (72 votes)
1 Apr 2016CPOL7 min read 277.8K   16.6K   114   66
This article presents a from-scratch implementation of cubic spline interpolation in C#.

Image 1

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:

C#
// 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:

C#
// 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:

C#
// 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:

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

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

Image 2

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:

C#
// 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:

Image 3

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:

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

Image 4

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

License

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


Written By
United States United States
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

 
GeneralMy vote of 5 Pin
Nicolae Fieraru22-Jul-21 16:23
Nicolae Fieraru22-Jul-21 16:23 
GeneralMy vote of 3 Pin
Ángel Gaspar16-Jan-21 7:11
Ángel Gaspar16-Jan-21 7:11 
GeneralMy vote of 3 Pin
Ángel Gaspar16-Jan-21 7:10
Ángel Gaspar16-Jan-21 7:10 
QuestionParabolic runouts Pin
Member 1359788928-Oct-20 10:06
Member 1359788928-Oct-20 10:06 
GeneralMy vote of 5 Pin
Member 1470970214-Jan-20 10:22
Member 1470970214-Jan-20 10:22 
QuestionClosed curve Pin
Member 1348614010-Sep-19 22:41
Member 1348614010-Sep-19 22:41 
Questionwhat for 3D points Pin
MayurVadher4-May-19 0:35
MayurVadher4-May-19 0:35 
Question3d Spline? Pin
JamesParsons5-Feb-19 19:17
JamesParsons5-Feb-19 19:17 
AnswerRe: 3d Spline? Pin
xiahangli10-Feb-20 8:24
xiahangli10-Feb-20 8:24 
GeneralRe: 3d Spline? Pin
Member 1172068130-Jan-21 6:33
Member 1172068130-Jan-21 6:33 
QuestionEvaluation at Fixed Intervals Pin
rem45acp8-Jan-19 8:59
rem45acp8-Jan-19 8:59 
GeneralMy vote of 5 Pin
BillWoodruff11-Nov-17 14:07
professionalBillWoodruff11-Nov-17 14:07 
QuestionYou saved my life thanks for the source Pin
amagitech23-Mar-17 22:39
amagitech23-Mar-17 22:39 
QuestionGood article. Pin
Douglas Smallish10-Feb-17 10:04
Douglas Smallish10-Feb-17 10:04 
QuestionWow, translated to C++, saved a tons of work Pin
Klik12-Dec-16 11:25
Klik12-Dec-16 11:25 
QuestionBugs evaluating beyond provided endpoints? Pin
Dave N 230-Oct-16 7:45
Dave N 230-Oct-16 7:45 
QuestionSplines with periodical boundary conditions Pin
LionAM1-Aug-16 23:13
LionAM1-Aug-16 23:13 
SuggestionPossible improvements Pin
LionAM13-Jul-16 0:09
LionAM13-Jul-16 0:09 
QuestionSmooth noisy data..? Pin
yltsa3-Apr-16 23:35
yltsa3-Apr-16 23:35 
Hi,

Is it possible to use this spline so that it is not going through the original points? So that we can for example smooth noisy data?

It that case we need an input parameter, which defines the strength of the smoothing.

BR,
Ilkka
AnswerRe: Smooth noisy data..? Pin
Ryan Seghers4-Apr-16 17:16
Ryan Seghers4-Apr-16 17:16 
GeneralRe: Smooth noisy data..? Pin
adamhancock19-Jan-17 23:10
adamhancock19-Jan-17 23:10 
QuestionSetting start and end slopes - Geometric fitting? Pin
lasilve21-Mar-16 15:23
lasilve21-Mar-16 15:23 
AnswerRe: Setting start and end slopes - Geometric fitting? Pin
Ryan Seghers2-Apr-16 6:35
Ryan Seghers2-Apr-16 6:35 
QuestionLive spline Pin
Jesús Álvarez26-Jan-16 3:30
Jesús Álvarez26-Jan-16 3:30 
AnswerRe: Live spline Pin
Ryan Seghers2-Apr-16 6:57
Ryan Seghers2-Apr-16 6:57 

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

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