Lagrange Interpolation over quadrilateral points

Question

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?

Answer 1

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).