# Lagrange Interpolation over quadrilateral points

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?

• A link or explanation of Lagrange interpolation would be nice. Jun 13 '18 at 16:08
• Done! I'm also looking for some plots but I can't find any...
– Dan
Jun 13 '18 at 16:41
• 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. Jun 13 '18 at 17:41
• 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)?
– user73941
Jun 14 '18 at 8:59
• Done! Let me know if there is something else I could improve.
– Dan
Jun 18 '18 at 8:11

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.

• Even better would be the new readonly ref struct :-o or the usage of the new in parameter Reference semantics with value types Jun 13 '18 at 17:35