6
\$\begingroup\$

http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

We have implemented it using the formula from Wikipedia :

public class MonotoneCubicSplineInterpolation
    {
        public static double[] Calc(double[] xs, double[] ys, double[] x_interp)
        {
            var length = xs.Length;

            // Deal with length issues
            if (length != ys.Length)
            {
                IPDevLoggerWrapper.Error("Need an equal count of xs and ys");
                throw new Exception("Need an equal count of xs and ys");
            }
            if (length == 0)
            {
                return null;
            }
            if (length == 1)
            {
                return new double[] {ys[0]};
            }


           // Get consecutive differences and slopes
            var delta = new double[length - 1];
            var m = new double[length];

            for (int i = 0; i < length - 1; i++)
            {
                delta[i] = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]);
                if (i > 0)
                {
                    m[i] = (delta[i - 1] + delta[i]) / 2;
                }
            }
            var toFix = new List<int>();
            for (int i = 1; i < length - 1; i++)
            {
                if ((delta[i] > 0 && delta[i - 1] < 0) || (delta[i] < 0 && delta[i - 1] > 0))
                {
                    toFix.Add(i);
                }
            }
            foreach (var val in toFix)
            {
                m[val] = 0;
            }

            m[0] = delta[0];
            m[length - 1] = delta[length - 2];

            toFix.Clear();
            for (int i = 0; i < length - 1; i++)
            {
                if (delta[i] == 0)
                {
                    toFix.Add(i);
                }
            }
            foreach (var val in toFix)
            {
                m[val] = 0;
                m[val + 1] = 0;
            }

            var alpha = new double[length - 1];
            var beta = new double[length - 1];
            var dist = new double[length - 1];
            var tau = new double[length - 1];
            for (int i = 0; i < length - 1; i++)
            {
                alpha[i] = m[i] / delta[i];
                beta[i] = m[i + 1] / delta[i];
                dist[i] = Math.Pow(alpha[i], 2) + Math.Pow(beta[i], 2);
                tau[i] = 3/Math.Sqrt(dist[i]);
            }

            toFix.Clear();
            for (int i = 0; i < length - 1; i++)
            {
                if (dist[i] > 9)
                {
                    toFix.Add(i);
                }
            }

            foreach (var val in toFix)
            {
                m[val] = tau[val] * alpha[val] * delta[val];
                m[val + 1] = tau[val] * beta[val] * delta[val];
            }

            var y_interp = new double[x_interp.Length];
            int ind = 0;

            foreach (var x in x_interp)
            {
                int i;
                for (i = xs.Length - 2; i >= 0; --i)
                {
                    if (xs[i] <= x)
                    {
                        break;
                    }
                }
                var h = xs[i + 1] - xs[i];
                var t = (x - xs[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;
                y_interp[ind++] = h00*ys[i] + h10*h*m[i] + h01*ys[i + 1] + h11*h*m[i + 1];

                continue;
            }

            return y_interp;
        }
    }

Please comment about style, correctness and complexity.

\$\endgroup\$
  • \$\begingroup\$ Regarding correctness: Where are your unit tests? This is trivial to write unit tests for. \$\endgroup\$ – Daniel Mann Dec 14 '14 at 18:19
  • \$\begingroup\$ @DanielMann it is tested the results are ok. i'm talking in a more general way. \$\endgroup\$ – Gilad Dec 14 '14 at 21:11
  • \$\begingroup\$ @Gilad - can you tell what is this x_interp and what values should this be? Thanks. \$\endgroup\$ – SpaceDog Apr 14 '19 at 18:12
  • \$\begingroup\$ @SpaceDog if I recall correctly it is an array of points which represent how many new points you would like to have, for example when fitting into a curve.. \$\endgroup\$ – Gilad Apr 15 '19 at 6:23
  • \$\begingroup\$ ah, thanks for the info. \$\endgroup\$ – SpaceDog Apr 15 '19 at 15:46
3
\$\begingroup\$

Don't throw Exception

throw new Exception("Need an equal count of xs and ys");

It forces client code to catch any subclass of Exception. In this case I would throw an ArgumentException.


The continue at the end of the last loop is redundant.


Here you're using double.Equals

if (delta[i] == 0)

From MSDN

The Equals method should be used with caution, because two apparently equivalent values can be unequal due to the differing precision of the two values.

That link covers two techniques for dealing with this.


As far as I can tell, toFix can be removed. For example,

var toFix = new List<int>();
for (int i = 1; i < length - 1; i++)
{
    if ((delta[i] > 0 && delta[i - 1] < 0) || (delta[i] < 0 && delta[i - 1] > 0))
    {
        toFix.Add(i);
    }
}
foreach (var val in toFix)
{
    m[val] = 0;
}

Can be rewritten as

for (int i = 1; i < length - 1; i++)
{
    if ((delta[i] > 0 && delta[i - 1] < 0) || (delta[i] < 0 && delta[i - 1] > 0))
    {
        m[i] = 0;
    }
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.