#313 Spline plot using C# Windows Forms App

Add ScottPlot.WinForms using Manage NuGet packages.

Form1.cs

using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using ScottPlot;

namespace Spline
{
    public partial class Form1 : Form
    {
        public static class Cubic
        {
            /// <summary>
            /// Generate a smooth (interpolated) curve that follows the path of the given X/Y points
            /// </summary>
            public static (double[] xs, double[] ys) InterpolateXY(double[] xs, double[] ys, int count)
            {
                if (xs is null || ys is null || xs.Length != ys.Length)
                    throw new ArgumentException($"{nameof(xs)} and {nameof(ys)} must have same length");

                int inputPointCount = xs.Length;
                double[] inputDistances = new double[inputPointCount];
                for (int i = 1; i < inputPointCount; i++)
                {
                    double dx = xs[i] - xs[i - 1];
                    double dy = ys[i] - ys[i - 1];
                    double distance = Math.Sqrt(dx * dx + dy * dy);
                    inputDistances[i] = inputDistances[i - 1] + distance;
                }

                double meanDistance = inputDistances.Last() / (count - 1);
                double[] evenDistances = Enumerable.Range(0, count).Select(x => x * meanDistance).ToArray();
                double[] xsOut = Interpolate(inputDistances, xs, evenDistances);
                double[] ysOut = Interpolate(inputDistances, ys, evenDistances);
                return (xsOut, ysOut);
            }

            private static double[] Interpolate(double[] xOrig, double[] yOrig, double[] xInterp)
            {
                (double[] a, double[] b) = FitMatrix(xOrig, yOrig);

                double[] yInterp = new double[xInterp.Length];
                for (int i = 0; i < yInterp.Length; i++)
                {
                    int j;
                    for (j = 0; j < xOrig.Length - 2; j++)
                        if (xInterp[i] <= xOrig[j + 1])
                            break;

                    double dx = xOrig[j + 1] - xOrig[j];
                    double t = (xInterp[i] - xOrig[j]) / dx;
                    double y = (1 - t) * yOrig[j] + t * yOrig[j + 1] +
                        t * (1 - t) * (a[j] * (1 - t) + b[j] * t);
                    yInterp[i] = y;
                }
                return yInterp;
            }
            private static (double[] a, double[] b) FitMatrix(double[] x, double[] y)
            {
                int n = x.Length;
                double[] a = new double[n - 1];
                double[] b = new double[n - 1];
                double[] r = new double[n];
                double[] A = new double[n];
                double[] B = new double[n];
                double[] C = new double[n];

                double dx1, dx2, dy1, dy2;

                dx1 = x[1] - x[0];
                C[0] = 1.0f / dx1;
                B[0] = 2.0f * C[0];
                r[0] = 3 * (y[1] - y[0]) / (dx1 * dx1);

                for (int i = 1; i < n - 1; i++)
                {
                    dx1 = x[i] - x[i - 1];
                    dx2 = x[i + 1] - x[i];
                    A[i] = 1.0f / dx1;
                    C[i] = 1.0f / dx2;
                    B[i] = 2.0f * (A[i] + C[i]);
                    dy1 = y[i] - y[i - 1];
                    dy2 = y[i + 1] - y[i];
                    r[i] = 3 * (dy1 / (dx1 * dx1) + dy2 / (dx2 * dx2));
                }

                dx1 = x[n - 1] - x[n - 2];
                dy1 = y[n - 1] - y[n - 2];
                A[n - 1] = 1.0f / dx1;
                B[n - 1] = 2.0f * A[n - 1];
                r[n - 1] = 3 * (dy1 / (dx1 * dx1));

                double[] cPrime = new double[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]);

                double[] dPrime = new double[n];
                dPrime[0] = r[0] / B[0];
                for (int i = 1; i < n; i++)
                    dPrime[i] = (r[i] - dPrime[i - 1] * A[i]) / (B[i] - cPrime[i - 1] * A[i]);

                double[] k = new double[n];
                k[n - 1] = dPrime[n - 1];
                for (int i = n - 2; i >= 0; i--)
                    k[i] = dPrime[i] - cPrime[i] * k[i + 1];

                for (int i = 1; i < n; i++)
                {
                    dx1 = x[i] - x[i - 1];
                    dy1 = y[i] - y[i - 1];
                    a[i - 1] = k[i - 1] * dx1 - dy1;
                    b[i - 1] = -k[i] * dx1 + dy1;
                }

                return (a, b);
            }
        }

        public Form1()
        {
            InitializeComponent();
        }

        private void Form1_Load(object sender, EventArgs e)
        {
            Random rand = new(1268);
            int pointCount = 15;
            double[] xs1 = DataGen.Consecutive(pointCount);
            double[] ys1 = new double[pointCount];
            for (int i = 1; i < pointCount; i++)
            {
                //xs1[i] = xs1[i - 1] + rand.NextDouble() - .5;
                ys1[i] = ys1[i - 1] + rand.NextDouble() - .5;
            }

            // Use cubic interpolation to smooth the original data
            (double[] xs2, double[] ys2) = Cubic.InterpolateXY(xs1, ys1, 200);

            // Plot the original vs. interpolated data
            var plt = new ScottPlot.Plot(900, 500);
            plt.AddScatter(xs1, ys1, label: "original", markerSize: 7);
            plt.AddScatter(xs2, ys2, label: "interpolated", markerSize: 3);
            plt.Style(Style.Seaborn);
            plt.Title("Cubic Spline");
            plt.XLabel("x axis");
            plt.YLabel("y axis");
            plt.Legend();
            plt.SaveFig("interpolation.png");
            pictureBox1.ImageLocation = "interpolation.png";
        }
    }
}