diff --git a/src/Examples/Interpolation/AkimaSpline.cs b/src/Examples/Interpolation/AkimaSpline.cs index 8789d7ee..b98ef15e 100644 --- a/src/Examples/Interpolation/AkimaSpline.cs +++ b/src/Examples/Interpolation/AkimaSpline.cs @@ -72,7 +72,7 @@ namespace Examples.InterpolationExamples Console.WriteLine(); // 2. Create akima spline interpolation - var method = new AkimaSplineInterpolation(points, values); + var method = CubicSpline.InterpolateAkima(points, values); Console.WriteLine(@"2. Create akima spline interpolation based on arbitrary points"); Console.WriteLine(); diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index 72e63eba..ab256bb1 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -54,20 +54,63 @@ namespace MathNet.Numerics } /// - /// Create a linear spline interpolation based on arbitrary points (sorted ascending). + /// Create a linear spline interpolation based on arbitrary points. /// - /// The sample points t, sorted ascending. Supports both lists and arrays. - /// The sample point values x(t). Supports both lists and arrays. + /// The sample points t. Optimized for arrays. + /// The sample point values x(t). Optimized for arrays. /// /// An interpolation scheme optimized for the given sample points and values, /// which can then be used to compute interpolations and extrapolations /// on arbitrary points. /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// public static IInterpolation LinearBetweenPoints(IEnumerable points, IEnumerable values) { return LinearSpline.Interpolate(points, values); } + /// + /// Create an Akima cubic spline interpolation based on arbitrary points. Akima splines are robust to outliers. + /// + /// The sample points t. Optimized for arrays. + /// The sample point values x(t). Optimized for arrays. + /// + /// An interpolation scheme optimized for the given sample points and values, + /// which can then be used to compute interpolations and extrapolations + /// on arbitrary points. + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + public static IInterpolation CubicBetweenPointsRobust(IEnumerable points, IEnumerable values) + { + return CubicSpline.InterpolateAkima(points, values); + } + + /// + /// Create a hermite cubic spline interpolation based on arbitrary points and their slopes/first derivative. + /// + /// The sample points t. Optimized for arrays. + /// The sample point values x(t). Optimized for arrays. + /// The slope at the sample points. Optimized for arrays. + /// + /// An interpolation scheme optimized for the given sample points and values, + /// which can then be used to compute interpolations and extrapolations + /// on arbitrary points. + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + public static IInterpolation CubicBetweenPointsWithDerivatives(IEnumerable points, IEnumerable values, IEnumerable firstDerivatives) + { + return CubicSpline.InterpolateHermite(points, values, firstDerivatives); + } + /// /// Create a floater hormann rational pole-free interpolation based on arbitrary points. /// @@ -102,4 +145,4 @@ namespace MathNet.Numerics return method; } } -} \ No newline at end of file +} diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 7f4f1908..9dba29a9 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -29,6 +29,9 @@ // using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.Interpolation { @@ -52,6 +55,11 @@ namespace MathNet.Numerics.Interpolation /// third order spline coefficients (N) public CubicSpline(double[] x, double[] c0, double[] c1, double[] c2, double[] c3) { + if (x.Length != c0.Length + 1 || x.Length != c1.Length + 1 || x.Length != c2.Length + 1 || x.Length != c3.Length + 1) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + _x = x; _c0 = c0; _c1 = c1; @@ -60,6 +68,123 @@ namespace MathNet.Numerics.Interpolation _indefiniteIntegral = new Lazy(ComputeIndefiniteIntegral); } + /// + /// Create a hermite cubic spline interpolation from a set of (x,y) value pairs and their slope (first derivative). + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + public static CubicSpline InterpolateHermite(IEnumerable x, IEnumerable y, IEnumerable firstDerivatives) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + var dd = (firstDerivatives as double[]) ?? firstDerivatives.ToArray(); + + if (xx.Length != yy.Length || xx.Length != dd.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + Sorting.Sort(xx, yy, dd); + + var c0 = new double[xx.Length - 1]; + var c1 = new double[xx.Length - 1]; + var c2 = new double[xx.Length - 1]; + var c3 = new double[xx.Length - 1]; + for (int i = 0; i < c1.Length; i++) + { + double w = xx[i + 1] - xx[i]; + double w2 = w*w; + c0[i] = yy[i]; + c1[i] = dd[i]; + c2[i] = (3*(yy[i + 1] - yy[i])/w - 2*dd[i] - dd[i + 1])/w; + c3[i] = (2*(yy[i] - yy[i + 1])/w + dd[i] + dd[i + 1])/w2; + } + + return new CubicSpline(xx, c0, c1, c2, c3); + } + + /// + /// Create an Akima cubic spline interpolation from a set of (x,y) value pairs. Akima splines are robust to outliers. + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + public static CubicSpline InterpolateAkima(IEnumerable x, IEnumerable y) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + + if (xx.Length != yy.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + Sorting.Sort(xx, yy); + + /* Prepare divided differences (diff) and weights (w) */ + + var diff = new double[xx.Length - 1]; + var weights = new double[xx.Length - 1]; + + for (int i = 0; i < diff.Length; i++) + { + diff[i] = (yy[i + 1] - yy[i])/(xx[i + 1] - xx[i]); + } + + for (int i = 1; i < weights.Length; i++) + { + weights[i] = Math.Abs(diff[i] - diff[i - 1]); + } + + /* Prepare Hermite interpolation scheme */ + + var dd = new double[xx.Length]; + + for (int i = 2; i < dd.Length - 2; i++) + { + dd[i] = weights[i - 1].AlmostEqual(0.0) && weights[i + 1].AlmostEqual(0.0) + ? (((xx[i + 1] - xx[i])*diff[i - 1]) + ((xx[i] - xx[i - 1])*diff[i]))/(xx[i + 1] - xx[i - 1]) + : ((weights[i + 1]*diff[i - 1]) + (weights[i - 1]*diff[i]))/(weights[i + 1] + weights[i - 1]); + } + + dd[0] = DifferentiateThreePoint(xx, yy, 0, 0, 1, 2); + dd[1] = DifferentiateThreePoint(xx, yy, 1, 0, 1, 2); + dd[xx.Length - 2] = DifferentiateThreePoint(xx, yy, xx.Length - 2, xx.Length - 3, xx.Length - 2, xx.Length - 1); + dd[xx.Length - 1] = DifferentiateThreePoint(xx, yy, xx.Length - 1, xx.Length - 3, xx.Length - 2, xx.Length - 1); + + /* Build Akima spline using Hermite interpolation scheme */ + + return InterpolateHermite(xx, yy, dd); + } + + /// + /// Three-Point Differentiation Helper. + /// + /// Sample Points t. + /// Sample Values x(t). + /// Index of the point of the differentiation. + /// Index of the first sample. + /// Index of the second sample. + /// Index of the third sample. + /// The derivative approximation. + static double DifferentiateThreePoint(double[] xx, double[] yy, int indexT, int index0, int index1, int index2) + { + double x0 = yy[index0]; + double x1 = yy[index1]; + double x2 = yy[index2]; + + double t = xx[indexT] - xx[index0]; + double t1 = xx[index1] - xx[index0]; + double t2 = xx[index2] - xx[index0]; + + double a = (x2 - x0 - (t2/t1*(x1 - x0)))/(t2*t2 - t1*t2); + double b = (x1 - x0 - a*t1*t1)/t1; + return (2*a*t) + b; + } + /// /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). /// diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs index 307252ba..c115ae26 100644 --- a/src/Numerics/Interpolation/LinearSpline.cs +++ b/src/Numerics/Interpolation/LinearSpline.cs @@ -51,6 +51,11 @@ namespace MathNet.Numerics.Interpolation /// Slopes (N) at the sample points (first order coefficients): N public LinearSpline(double[] x, double[] c0, double[] c1) { + if (x.Length != c0.Length + 1 && x.Length != c0.Length || x.Length != c1.Length + 1) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + _x = x; _c0 = c0; _c1 = c1; diff --git a/src/Numerics/Interpolation/QuadraticSpline.cs b/src/Numerics/Interpolation/QuadraticSpline.cs index 948fe518..0fce11af 100644 --- a/src/Numerics/Interpolation/QuadraticSpline.cs +++ b/src/Numerics/Interpolation/QuadraticSpline.cs @@ -29,6 +29,7 @@ // using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.Interpolation { @@ -50,6 +51,11 @@ namespace MathNet.Numerics.Interpolation /// second order spline coefficients (N) public QuadraticSpline(double[] x, double[] c0, double[] c1, double[] c2) { + if (x.Length != c0.Length + 1 || x.Length != c1.Length + 1 || x.Length != c2.Length + 1) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + _x = x; _c0 = c0; _c1 = c1; diff --git a/src/UnitTests/InterpolationTests/AkimaSplineTest.cs b/src/UnitTests/InterpolationTests/AkimaSplineTest.cs index acc97113..c092232f 100644 --- a/src/UnitTests/InterpolationTests/AkimaSplineTest.cs +++ b/src/UnitTests/InterpolationTests/AkimaSplineTest.cs @@ -33,21 +33,11 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.InterpolationTests { - /// - /// AkimaSpline test case. - /// [TestFixture, Category("Interpolation")] public class AkimaSplineTest { - /// - /// Sample points. - /// readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; - - /// - /// Sample values. - /// - readonly double[] _x = { 1.0, 2.0, -1.0, 0.0, 1.0 }; + readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; /// /// Verifies that the interpolation matches the given value at all the provided sample points. @@ -55,11 +45,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [Test] public void FitsAtSamplePoints() { - IInterpolation interpolation = new AkimaSplineInterpolation(_t, _x); - - for (int i = 0; i < _x.Length; i++) + IInterpolation it = CubicSpline.InterpolateAkima(_t, _y); + for (int i = 0; i < _y.Length; i++) { - Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); + Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i); } } @@ -78,12 +67,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests [TestCase(1.2, 0.2, 1e-15)] [TestCase(10.0, 9, 1e-15)] [TestCase(-10.0, -151, 1e-15)] - public void FitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) + public void FitsAtArbitraryPoints(double t, double x, double maxAbsoluteError) { - IInterpolation interpolation = new AkimaSplineInterpolation(_t, _x); + IInterpolation it = CubicSpline.InterpolateAkima(_t, _y); // TODO: Verify the expected values (that they are really the expected ones) - Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); + Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t); } /// @@ -97,10 +86,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { double[] x, y, xtest, ytest; LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); - IInterpolation interpolation = new AkimaSplineInterpolation(x, y); + IInterpolation it = CubicSpline.InterpolateAkima(x, y); for (int i = 0; i < xtest.Length; i++) { - Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i); + Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i); } } }