diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index af551b96..7c81c638 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -201,7 +201,7 @@ namespace MathNet.Numerics } /// - /// Create an piecewise cubic Akima spline interpolation based on arbitrary points. + /// Create a piecewise cubic Akima spline interpolation based on arbitrary points. /// Akima splines are robust to outliers. /// /// The sample points t. @@ -221,6 +221,27 @@ namespace MathNet.Numerics return Interpolation.CubicSpline.InterpolateAkima(points, values); } + /// + /// Create a piecewise cubic monotone spline interpolation based on arbitrary points. + /// This is a shape-preserving spline with continuous first derivative. + /// + /// The sample points t. + /// The sample point values x(t). + /// + /// An interpolation scheme optimized for the given sample points and values, + /// which can then be used to compute interpolations and extrapolations + /// on arbitrary points. + /// + /// + /// if your data is already sorted in arrays, consider to use + /// MathNet.Numerics.Interpolation.CubicSpline.InterpolatePchipSorted + /// instead, which is more efficient. + /// + public static IInterpolation CubicSplineMonotone(IEnumerable points, IEnumerable values) + { + return Interpolation.CubicSpline.InterpolatePchip(points, values); + } + /// /// Create a piecewise cubic Hermite spline interpolation based on arbitrary points /// and their slopes/first derivative. diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 5dc52146..cba2c54c 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -244,7 +244,9 @@ namespace MathNet.Numerics.Interpolation var mIs0 = m[i].AlmostEqual(0.0); if (mIs0 || mPrevIs0 || Math.Sign(m[i]) != Math.Sign(m[i - 1])) + { dd[i] = 0; + } else { // Weighted harmonic mean of each slope. @@ -272,10 +274,14 @@ namespace MathNet.Numerics.Interpolation var d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1); if (Math.Sign(d) != Math.Sign(m0)) + { return 0.0; + } if (Math.Sign(m0) != Math.Sign(m1) && (Math.Abs(d) > 3 * Math.Abs(m0))) + { return 3 * m0; + } return d; }