Bug in cubic spline interpolator?

[swharden]

@Marco-Sanguinetti recently commented the following on a blog post of mine swharden/SWHarden.com#27

Now arrays A, B, and r are defined for index "n-1" but array C is not. However, C[n-1] is required in the indicated for loop.

The referenced code is also in ScottPlot, here

ScottPlot/src/ScottPlot4/ScottPlot/Statistics/Interpolation/Cubic.cs

Lines 99 to 106 in 093994e

	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]);

This issue will track better understanding this issue and coming up with a fix. I'll ping @drolevar in case they're interested in providing feedback, but they recently contributed some great related work in #3623 #3629

[swharden]

Adding more context, I think the issue may be that C[n-1] never gets set in the loop above, so it's always 0

ScottPlot/src/ScottPlot4/ScottPlot/Statistics/Interpolation/Cubic.cs

Lines 85 to 95 in 093994e

	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));
	}

Perhaps the fix is to add:

C[n - 1] = 1.0f / (x[n] - x[n-1]);

[swharden]

Actually I'm not going to worry about this. That code is only in ScottPlot 4. There is a more elegant method in ScottPlot 5 that leans on SkiaSharp

ScottPlot/src/ScottPlot5/ScottPlot5/PathStrategies/CubicSpline.cs

Lines 1 to 85 in 34393d9

	namespace ScottPlot.PathStrategies;
	using System.Linq;
	/// <summary>
	/// Connect points with Catmull-Rom cubic splines, see https://doi.org/10.1007/s42979-021-00770-x
	/// NaN values will be skipped, producing a gap in the path.
	/// </summary>
	public class CubicSpline : IPathStrategy
	{
	public double Tension { get; set; } = 1.0f;
	/// <summary>
	/// Organize into segments of connected points padded with empty points at the front and back
	/// </summary>
	/// <param name="pixels">Pixels</param>
	/// <returns>IEnumerable of arrays representing the segments</returns>
	private IEnumerable<Pixel[]> GetSegments(IEnumerable<Pixel> pixels)
	{
	// Reserve an empty element at the front
	List<Pixel> segment = [Pixel.NaN];
	foreach (var pixel in pixels)
	{
	if (float.IsNaN(pixel.X) || float.IsNaN(pixel.Y))
	{
	if (segment.Count > 1)
	{
	// And one at the back
	segment.Add(Pixel.NaN);
	yield return segment.ToArray();
	segment.Clear();
	segment.Add(Pixel.NaN);
	}
	}
	else
	{
	segment.Add(pixel);
	}
	}
	if (segment.Count > 1)
	{
	segment.Add(Pixel.NaN);
	yield return segment.ToArray();
	}
	}
	public SKPath GetPath(IEnumerable<Pixel> pixels)
	{
	SKPath path = new();
	// Organize into segments of connected points
	foreach (var segment in GetSegments(pixels).Where(seg => seg.Length >= 4))
	{
	Pixel first = segment[1];
	Pixel second = segment[2];
	Pixel nextLast = segment[^3];
	Pixel last = segment[^2];
	bool closed = first == last;
	// Fill in the empty padded points
	segment[0] = closed ? nextLast : first;
	segment[^1] = closed ? second : last;
	path.MoveTo(first.ToSKPoint());
	for (int i = 2; i < segment.Length - 1; ++i)
	{
	Pixel p0 = segment[i - 2];
	Pixel p1 = segment[i - 1];
	Pixel p2 = segment[i];
	Pixel p3 = segment[i + 1];
	Pixel c1 = p1 + (p2 - p0) / (6.0f * (float)Tension);
	Pixel c2 = p2 - (p3 - p1) / (6.0f * (float)Tension);
	path.CubicTo(c1.ToSKPoint(), c2.ToSKPoint(), p2.ToSKPoint());
	}
	}
	return path;
	}
	}