Spline interpolation in C#

Interpolator.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace Interpolation
{
    /// <summary>
    /// Spline interpolation class.
    /// </summary>
    public class SplineInterpolator
    {
        private readonly double[] _keys;

        private readonly double[] _values;
        
        private readonly double[] _h;
        
        private readonly double[] _a;

        /// <summary>
        /// Class constructor.
        /// </summary>
        /// <param name="nodes">Collection of known points for further interpolation.
        /// Should contain at least two items.</param>
        public SplineInterpolator(IDictionary<double, double> nodes)
        {
            if (nodes == null)
            {
                throw new ArgumentNullException("nodes");
            }

            var n = nodes.Count;

            if (n < 2)
            {
                throw new ArgumentException("At least two point required for interpolation.");
            }

            _keys = nodes.Keys.ToArray();
            _values = nodes.Values.ToArray();
            _a = new double[n];
            _h = new double[n];

            for (int i = 1; i < n; i++)
            {
                _h[i] = _keys[i] - _keys[i - 1];
            }

            if (n > 2)
            {
                var sub = new double[n - 1];
                var diag = new double[n - 1];
                var sup = new double[n - 1];

                for (int i = 1; i <= n - 2; i++)
                {
                    diag[i] = (_h[i] + _h[i + 1]) / 3;
                    sup[i] = _h[i + 1] / 6;
                    sub[i] = _h[i] / 6;
                    _a[i] = (_values[i + 1] - _values[i]) / _h[i + 1] - (_values[i] - _values[i - 1]) / _h[i];
                }

                SolveTridiag(sub, diag, sup, ref _a, n - 2);
            }
        }

        /// <summary>
        /// Gets interpolated value for specified argument.
        /// </summary>
        /// <param name="key">Argument value for interpolation. Must be within 
        /// the interval bounded by lowest ang highest <see cref="_keys"/> values.</param>
        public double GetValue(double key)
        {
            int gap = 0;
            var previous = double.MinValue;

            // At the end of this iteration, "gap" will contain the index of the interval
            // between two known values, which contains the unknown z, and "previous" will
            // contain the biggest z value among the known samples, left of the unknown z
            for (int i = 0; i < _keys.Length; i++)
            {
                if (_keys[i] < key && _keys[i] > previous)
                {
                    previous = _keys[i];
                    gap = i + 1;
                }
            }

            var x1 = key - previous;
            var x2 = _h[gap] - x1;

            return ((-_a[gap - 1] / 6 * (x2 + _h[gap]) * x1 + _values[gap - 1]) * x2 +
                (-_a[gap] / 6 * (x1 + _h[gap]) * x2 + _values[gap]) * x1) / _h[gap];
        }


        /// <summary>
        /// Solve linear system with tridiagonal n*n matrix "a"
        /// using Gaussian elimination without pivoting.
        /// </summary>
        private static void SolveTridiag(double[] sub, double[] diag, double[] sup, ref double[] b, int n)
        {
            int i;

            for (i = 2; i <= n; i++)
            {
                sub[i] = sub[i] / diag[i - 1];
                diag[i] = diag[i] - sub[i] * sup[i - 1];
                b[i] = b[i] - sub[i] * b[i - 1];
            }

            b[n] = b[n] / diag[n];
            
            for (i = n - 1; i >= 1; i--)
            {
                b[i] = (b[i] - sup[i] * b[i + 1]) / diag[i];
            }
        }
    }
}

Program.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace Interpolation
{
    class Program
    {
        static void Main(string[] args)
        {
            var r = 1000;
            var known = new Dictionary<double, double>
                { 
                    { 0.0, 0.0 },
                    { 100.0, 0.50 * r },
                    { 300.0, 0.75 * r },
                    { 500.0, 1.00 * r },
                };

            foreach (var pair in known)
            {
                Debug.WriteLine(String.Format("{0:0.000}\t{1:0.000}", pair.Key, pair.Value));
            }

            var scaler = new SplineInterpolator(known);
            var start = known.First().Key;
            var end = known.Last().Key;
            var step = (end - start) / 50;

            for (var x = start; x <= end; x += step)
            {
                var y = scaler.GetValue(x);
                Debug.WriteLine(String.Format("\t\t{0:0.000}\t{1:0.000}", x, y));
            }
        }
    }
}

EDIT: The more I test this with random sequential data the more I think this code doesn't actually work 😅

@swharden I'm looking for an easy/existing of a cubic spline implementation that I can take apart and reuse. So far the best I've found is this https://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation which still isn't great as it's written weirdly but at least all the math is done and works.

@swharden I'm looking for an easy/existing of a cubic spline implementation that I can take apart and reuse. So far the best I've found is this https://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation which still isn't great as it's written weirdly but at least all the math is done and works.

Hey @crener, since I last wrote that I found one and implemented it in C# as a portable little module. Check out the notes and resources here ScottPlot/ScottPlot#459

The standalone module is here:
https://github.com/swharden/ScottPlot/tree/404105e5d7ae8399b2e40e9bd64b246d3b3b80dd/src/ScottPlot/Statistics/Interpolation

and I use it like this:

double[] xs = { 0, 10, 20, 30 };
double[] ys = { 65, 25, 55, 80 };

var nsi = new ScottPlot.Statistics.Interpolation.NaturalSpline(xs, ys, resolution: 20);
var psi = new ScottPlot.Statistics.Interpolation.PeriodicSpline(xs, ys, resolution: 20);
var esi = new ScottPlot.Statistics.Interpolation.EndSlopeSpline(xs, ys, resolution: 20);

var plt = new ScottPlot.Plot();
plt.PlotScatter(xs, ys, Color.Black, markerSize: 10, lineWidth: 0, label: "Original Data");
plt.PlotScatter(nsi.interpolatedXs, nsi.interpolatedYs, Color.Red, markerSize: 3, label: "Natural Spline");
plt.PlotScatter(psi.interpolatedXs, psi.interpolatedYs, Color.Green, markerSize: 3, label: "Periodic Spline");
plt.PlotScatter(esi.interpolatedXs, esi.interpolatedYs, Color.Blue, markerSize: 3, label: "End Slope Spline");
plt.Legend();

To create this output:

Oh wow, http://www.mosismath.com/NaturalSplines/NaturalSplines.html Really is a way better example (plus the added example of being in the goto graphing library for c#). I hope that that people looking at this will go over there instead (or the search engines put that result a little higher) :)

Thanks

@swharden Thanks for the work on that. I am trying to use your package in Unity3D but the game crashes when I try to run

var psi = new ScottPlot.Statistics.Interpolation.PeriodicSpline(xs, ys, resolution: 20);

I downloaded your package through Nuget and put the .dll with its' dependency at my Unity3D Plugins folder. I also use the ScottPlot.Statistics.Interpolation namespace in the script I want to perform the interpolation.

Do you know if I am doing something wrong or is there something I do not understand?

Thanks a lot in advance.

@iliachry Ha, ironically thats what I did, ended up making this https://github.com/crener/Dots-Spline though I did basically refactor the code a lot to integrate it into the existing spline system I had :)

@iliachry I understand you are trying to use ScottPlot's Statistics.Interpolation module in Unity. You said "the game crashes", but that's not really enough information to facilitate a helpful response. Is an error message displayed? Does ScottPlot throw an exception, or is the problem elsewhere in your code?

If you're not sure how to use general NuGet packages in Unity, that topic is probably outside the scope of spline interpolation and this gist. @crener's Dots-Spline package for Unity may be a simpler starting point.

On the other hand if you found a list of X and Y pairs that causes ScottPlot's spine interpolation method to fail, copy/paste those values here or post those values in an issue on ScottPlot's GitHub and I'd be happy to take a look!

@iliachry one more thought is that if your xs or ys arrays contain double.NaN, double.PositiveInfinity, or double.NegativeInfinity, it might explain an exception thrown during interpolation...

A closer look at what values you're sending into Statistics.Interpolation.PeriodicSpline method may be helpful 👍

@crener thanks a lot. I'll give it a try if I won't make it with ScottPlott.

@swharden finally I found out that Unity crashes because I wanted to Interpolate 10^4 samples. I just now use 5 samples and Unity doesn't crash. Regarding your questions, Unity doesn't throw any exception or displays an error message. However, I don't get a smooth signal yet. I will search a little bit more what's happening and I will reply here if it doesn't finally work. Thanks a lot for your quick response.

@iliachry thanks for following up! I hadn't tested my implementation against a large number of samples like that. I suspect ScottPlot's spline interpolation module isn't optimized for large datasets, so your math problem took so long to solve that it gave the appearance of "Unity crashing". Good luck on your project! 👍

@swharden the problem was that Unity didn't write in the console something I wanted, which was in the same script before the creation of the spline. Anyway, now it works very well. Thank you very much for the implementation.

I recently found a more reliable method for achieving cubic spline interpolation in C#. The code I'm linking to here is simpler, faster, stateless, and doesn't throw exceptions.

https://swharden.com/blog/2022-01-22-spline-interpolation