8
\$\begingroup\$

I've been working on a program which calculates, given a point and 4 surrounding points, the Lagrange polynomial, in order to interpolate a value. Consider that I'm not a mathematician and I better understand code than formulas.

I've came up with the following code, which works, but I really don't think is general (and it's pretty ugly to me).

using System;
using System.Collections.Generic;
using System.Linq;

namespace LagrangeInterpolation
{
    public class Point : ICloneable
    {
        public double X { get; set; }

        public double Y { get; set; }

        public double Value { get; set; }   

        public Point Clone()
        {
            return (Point)this.MemberwiseClone();
        }           
    }

    public static class Lagrange
    {
        public static Point Interpolate(Point[] controlPoints, Point point)
        {           
            var A = -controlPoints[0].X - controlPoints[1].X + controlPoints[2].X + controlPoints[3].X;
            var B = -controlPoints[0].X + controlPoints[1].X + controlPoints[2].X - controlPoints[3].X;
            var C = +controlPoints[0].X - controlPoints[1].X + controlPoints[2].X - controlPoints[3].X;
            var X = 4 * point.X - controlPoints.Sum(x => x.X);

            var D = -controlPoints[0].Y - controlPoints[1].Y + controlPoints[2].Y + controlPoints[3].Y;
            var E = -controlPoints[0].Y + controlPoints[1].Y + controlPoints[2].Y - controlPoints[3].Y;
            var F = +controlPoints[0].Y - controlPoints[1].Y + controlPoints[2].Y - controlPoints[3].Y;
            var Y = 4 * point.Y - controlPoints.Sum(x => x.Y);

            var r = (X / 4 - B * Y / 4 * E) / (1 - D / 4 * E);
            var s = (Y - D * r) / E;

            var prevR = 0d;
            var prevS = 0d;
            const double precision = 0.00000001;
            while (!(prevR - r < precision && prevS - s < precision))
            {
                prevR = r;
                prevS = s;
                r = (X - B * s) / (A + C * s);
                s = (Y - D * r) / (E + F * r);
            }

            // Interpolate value
            var result = point.Clone();
            result.Value = ((1 - r) * (1 - s) * controlPoints[0].Value + (1 - r) * (1 + s) * controlPoints[1].Value + (1 + r) * (1 + s) * controlPoints[2].Value + (1 + r) * (1 - s) * controlPoints[3].Value) / 4;

            return result;
        }
    }
}

The input parameters are:

  • controlPoints: the 4 points, each one with its Value.
  • point: the point for which we want to calculate the interpolated value

The returned Point is a clone of the Input point with the Value property set. Every instance of Point have X and Y normalised within range -1..1 (I subtract the quadrilateral center from each point).

Example:

controlPoints = new [] {
    new Point() { X = -0.033675000000000566, Y = -0.02564999999999884, Value = 1.2787 },
    new Point() { X = -0.035524999999999807, Y = 0.024329999999999075, Value = 1.329 },
    new Point() { X = 0.03370499999999943, Y = 0.02564999999999884, Value = 1.3376 },
    new Point() { X = 0.035494999999999166, Y = -0.024329999999999075, Value = 1.302 }
}

point = new Point() { X = 0.018148174616284152, Y = -0.014201699949808244 }

Expected result is Point.Value = 1.3044829106888913

Can someone suggest a better way (formally and mathematically) to perform this calculation?

\$\endgroup\$
  • \$\begingroup\$ A link or explanation of Lagrange interpolation would be nice. \$\endgroup\$ – paparazzo Jun 13 '18 at 16:08
  • \$\begingroup\$ Done! I'm also looking for some plots but I can't find any... \$\endgroup\$ – Dan Jun 13 '18 at 16:41
  • \$\begingroup\$ I have rolled back your last edit. Please don't change or add to the code in your question after you have received answers. See What should I do when someone answers my question? Thank you. \$\endgroup\$ – Phrancis Jun 13 '18 at 17:41
  • \$\begingroup\$ What do you mean by "Point have x and y normalized within range -1..1"? Could you update the question with a data set and the desired result (It is allowed to update the question with additional information, but not to change or add code blocks to be reviewed)? \$\endgroup\$ – Henrik Hansen Jun 14 '18 at 8:59
  • \$\begingroup\$ Done! Let me know if there is something else I could improve. \$\endgroup\$ – Dan Jun 18 '18 at 8:11
3
\$\begingroup\$

I would use a struct with copy constructor instead of a class with clone. Performance wise it would be better and would be more readable.

Extract function for your condition in the while.

Add explicit parenthesis on complex arithmetic expressions (normal people don’t read operator precedence fluently). ( I think about this part B * Y / 4 )

Seems like you do six times the same thing with your points. Extract method somehow and name it appropriately (cross or whatever it is).

\$\endgroup\$
  • \$\begingroup\$ Even better would be the new readonly ref struct :-o or the usage of the new in parameter Reference semantics with value types \$\endgroup\$ – t3chb0t Jun 13 '18 at 17:35
  • \$\begingroup\$ About "six times the same thing with your points" I think you refer to the A,B,C and D,E,F variables. I think it's a Vector operation but I don't know the definition. I would like to have it in order to extract something meaningful. For the rest, I'll do my best to clean up the code and update question. Thanks \$\endgroup\$ – Dan Jun 13 '18 at 17:35
  • \$\begingroup\$ @dna2 please do not update your question with new code ;-) see What should I do when someone answers my question?: Do not add an improved version of the code after receiving an answer. Including revised versions of the code makes the question confusing, especially if someone later reviews the newer code. \$\endgroup\$ – t3chb0t Jun 13 '18 at 17:39
  • \$\begingroup\$ @t3chb0t understood, Voted the answer and won't update the code. I think the answer must be improved however. I'm looking for a more general approach to the Lagrange interpolation than syntax/style improvements (which are really welcome). \$\endgroup\$ – Dan Jun 13 '18 at 17:44
  • \$\begingroup\$ I’ll do my best to improve the answer from my desktop. \$\endgroup\$ – Regis Portalez Jun 13 '18 at 18:58

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.