MonotoneCubicInterpolation

MonotoneCubicInterpolation.cs

    public class MonotoneCubicSpline
    {
        private readonly double[] _x;
        private readonly double[] _y;
        private readonly int _n;
        private readonly double[] _m;

        public MonotoneCubicSpline(IList<Point> points)
        {
            _x = points.OrderBy(p => p.X).Select(p => p.X).ToArray();
            _y = points.OrderBy(p => p.X).Select(p => p.Y).ToArray();
            _n = points.Count;
            _m = new double[_n];

            var delta = new double[_n - 1];

            for (var k = 0; k < _n - 1; k++)
            {
                // 1. Find the slope of the section.
                delta[k] = (_y[k + 1] - _y[k]) / (_x[k + 1] - _x[k]);

                // 2. Obtain an initial estimate of the control point slope.
                if (k == 0)
                {
                    _m[k] = delta[k];
                }
                else
                {
                    if ((delta[k] > 0 && delta[k - 1] < 0) || (delta[k] < 0 && delta[k - 1] > 0))
                    {
                        _m[k] = 0;
                    }
                    else
                    {
                        _m[k] = (delta[k - 1] + delta[k]) / 2;
                    }
                }
            }

            _m[0] = delta[0];
            _m[_n - 1] = delta[_n - 2];

            // Improve control point slope.
            for (var k = 0; k < _n -1; k++)
            {
                // 3. Make the horizontal part horizontal.
                if (Math.Abs(delta[k]) < 1e-6)
                {
                    _m[k] = 0;
                    _m[k + 1] = 0;

                    continue;
                }

                // 4. Make the extreme points horizontal. (a, b> = 0)
                var a = _m[k] / delta[k];
                var b = _m[k + 1] / delta[k];

                if (a < 0 || b < 0)
                {
                    _m[k] = 0;
                }

                // 5. Reduce overshoot. (a, b <= 3)
                var d = Math.Pow(a, 2) + Math.Pow(b, 2);
                if (d > 9)
                {
                    var tau = 3 / Math.Sqrt(d);

                    _m[k] = tau * a * delta[k];
                    _m[k + 1] = tau * b * delta[k];
                }
            }
        }

        public double Interpolate(double x)
        {
            if (x >= _x[_n - 1])
            {
                return (float) _y[_n - 1];
            }

            if (x <= _x[0])
            {
                return (float) _y[0];
            }

            var i = 0;
            for (i = _x.Length - 1; i >= 0; i--)
            {
                if (_x[i] <= x)
                {
                    break;
                }
            }

            var h = _x[i + 1] - _x[i];
            var t = (x - _x[i]) / h;

            var t2 = Math.Pow(t, 2);
            var t3 = Math.Pow(t, 3);

            var h00 = 2 * t3 - 3 * t2 + 1;
            var h10 = t3 - 2 * t2 + t;
            var h01 = -2 * t3 + 3 * t2;
            var h11 = t3 - t2;

            var v = h00 * _y[i] + h10 * h * _m[i] + h01 * _y[i + 1] + h11 * h * _m[i + 1];

            return v;
        }
    }