Click here to Skip to main content
Click here to Skip to main content

Curve Fitting using Lagrange Interpolation

By , 12 Sep 2008
 

Introduction

In this article, I will explain curve fitting using the Lagrange interpolation polynomial. Curve fitting is used in a wide spectrum in engineering applications such as cars and air crafts surface design. The main problem is, given a set of points in the plan, we want to fit them in a smooth curve that passes through these points. The order of the curve f(x) depends on the number of points given. For example, if we have two points, the function will be of the first order, and the curve will be the line that passes through these two points, while if you have three points, the function will be of the second order f(x) = x2 . Let's first explain the Lagrange polynomial, then we will proceed to the algorithm and the implementation. In this article, I am using C# for coding.

Background

Lagrange Polynomial

An interpolation on two points, (x0, y0) and (x1, y1), results in a linear equation or a straight line. The standard form of a linear equation is given by y = mx + c, where m is the gradient of the line and c is the y-intercept.

m = y1 − y0 / x1 − x0
c = y0 − mx0

which results in:

y = (y1 − y0 / x1 − x0) * x + (x1 y0 − x0 y1 )/ (x1 − x0)

The linear equation is rewritten so that the two interpolated points, (x0, y0) and (x1, y1), are directly represented.

With this in mind, the linear equation is rewritten as:

P1(x) = a0(x − x1) + a1(x − x0)

where a0 and a1 are constants . The points x0 and x1 in the factors of the above equation are called the centers. Applying the equation at (x0, y0), we obtain:

y0 = a0(x0 − x1) + a1(x0 − x0)

or

a0 = y0/ x0−x1 

At (x1, y1), we get:

y1 = a0(x1 − x1) + a1(x1 − x0), or a1 = y1/ x1−x0

Therefore, the linear equation becomes:

P1(x) = y0 (x− x1) /(x0 − x1) + y1 (x− x0) / (x1 − x0)

The quadratic form of the Lagrange polynomial interpolates three points, (x0, y0), (x1, y1), and (x2, y2). The polynomial has the form:

P2(x) = a0(x− x1)(x− x2) + a1(x− x0) (x− x2) + a2(x− x0)(x− x1)

with centers at x0, x1, and x2. At (x0, y0):

y0 = a0(x0 − x1)(x0 − x2) + a1(x0 − x0)(x0 − x2) + a2(x0 − x0)(x0 − x1),

or

a0 = y0/ (x0 − x1)(x0 − x2)
a1 = y1 / (x1 − x0)(x1 − x2)
a2 = y2/ (x2 − x0)(x2 − x1)
P2(x) = y0 (x− x1)(x− x2) / (x0 − x1)(x0 − x2) + y1 (x − x0)(x − x2)/ 
       (x1 − x0)(x1 − x2) + y2 (x− x0)(x− x1)/ (x2 − x0)(x2 − x1)

In general, a Lagrange polynomial of degree n is a polynomial that is produced from an interpolation over a set of points, (xi , yi ) for i = 0, 1, . . ., n, as follows:

Pn(x) = y0L0(x) + y1L1(x) + ··· + ynLn(x)

Using the Code

The Algorithm

Given the interpolating points (xi , yi ) for i =0, 1, . . . ,n;

for i = 0 to n
  //the cumulative multiplication from k = 1 (and k not equal i) to n 
  Evaluate Li (x) = ∏nk=1,k != i (x−xk ) / (xi−xk ) ;
endfor 

Evaluate Pn(x) = y0L0(x) + y1L1(x) + ··· + ynLn(x);

Download the source code from the top of this page to view the code.

Points of Interest

Numerical computing is a very interesting field in software development. This article is related to that field ... And, I will post more articles soon about Computational Geometry ...such as Convex Hull algorithm and Triangulation.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

About the Author

Fady Aladdin
Software Developer
Egypt Egypt
Member
No Biography provided

Sign Up to vote   Poor Excellent
Add a reason or comment to your vote: x
Votes of 3 or less require a comment

Comments and Discussions

 
You must Sign In to use this message board.
Search this forum  
    Spacing  Noise  Layout  Per page   
GeneralA compact versionmembergustavtr20 Feb '12 - 11:24 
This article was very helpful. Nice work! Here is a compact version.
 
/// <summary>
/// Uses Lagrange interpolation to predict y.
/// Note! Make sure that the input point list does not contain points with the same x.
/// </summary>
public int PredictY(List<Point> points, int x)
{
    if (points.Count == 0) return 0;
 
    double[] lagrange = new double[points.Count];
    for (int i = 0; i < points.Count; i++)
    {
        lagrange[i] = 1;
 
        for (int k = 0; k < points.Count; k++)
        {
            if (i != k)
            {
                lagrange[i] *= ((double)(x - points[k].X)) / (points[i].X - points[k].X);
            }
        }
    }
 
    double y = 0;
    for (int i = 0; i < points.Count; i++)
    {
        y += lagrange[i] * points[i].Y;
    }
 
    return Convert.ToInt32(y);
}
Regards,
Gustav

GeneralMy vote of 5memberMarkDaniel5 Jan '12 - 12:21 
Thanks for sharing
Questionhow to get all pixels coordinate from curve fitting using Lagrange interpolationan on an imagememberMrKyaw24 Oct '11 - 16:18 
Dear Fady,
 

 
The topic here is quite useful for me and I am tryig to get the x, y coordinates in terms of pixels form the user drawn curve on the image.
 
My targeted output was show in this link
http://imageshack.us/photo/my-images/4/currentdetectedboundary.png/ http://imageshack.us/photo/my-images/196/targetedoutput.png/Seems http://imageshack.us/photo/my-images/98/originalboundarymanualh.png/
 
and integrate or using it tough for me as well as took so much time to do that.
 
Please advise me any alternative simple ways since I was stuck on this stage one week plus and rushing on it.
 

 

 
Thanks and best regards
 
Kyaw
GeneralRadius of curvememberM i s t e r L i s t e r25 Aug '10 - 9:48 
Hey great article and code...
 
My question: Can you get the radius of the curve that is defined by three sequential points along the line generated by your code as well?
Generalneed help plzmemberd_e_e_o6 Aug '10 - 4:33 
i have a data but i can't programming it. with C#. my problem is,i need a C# source code that can calculate and draw a curve with "interpolation polynomial" below from list
(1,1000),(2,1300),(3,2000),(4,3450),(5,4000),(6,n)
and what is the value of n ....
please anyone help me
GeneralCorrected an Overflow errormemberSBGTrading22 Jul '10 - 5:05 
Thank you for the code on this one...very helpful in my research.
 
I did make a change to the "renderCurve" method to prevent an overflow error I was getting on the g.DrawLines(new Pen(Color.Fuchsia), CurveArray); statement.
 
This overflow occurs if your datapoints are too close together in the X direction, and too far apart in the Y direction.
 
Here's the new renderCurve method I've used:
 

private void renderCurve()
{
try{
List<PointF> Curve = new List<PointF>();
for (int i = 0; i < Xs.Count; i++)
{
PointF p = new PointF(Xs[i], Ys[i]);
Curve.Add(p);
}
string msg = string.Empty;
int index = Curve.Count;
PointF[] CurveArray = new PointF[index];
CurveArray = Curve.ToArray();
for(int idx = 0; idx<CurveArray.Length; idx++) {
if(float.IsNegativeInfinity(CurveArray[idx].Y)) {
CurveArray[idx]=new PointF(CurveArray[idx].X,0.0f);
msg = string.Concat(msg,Environment.NewLine,idx.ToString()+"-element was NegativeInfinity");
}
if(float.IsPositiveInfinity(CurveArray[idx].Y)) {
CurveArray[idx]=new PointF(CurveArray[idx].X,g.VisibleClipBounds.Height);
msg = string.Concat(msg,Environment.NewLine,idx.ToString()+"-element was PositiveInfinity");
}
if(float.IsNaN(CurveArray[idx].Y)) {
CurveArray[idx]=new PointF(CurveArray[idx].X,0.0f);
msg = string.Concat(msg,Environment.NewLine,idx.ToString()+"-element was NaN");
}
if(CurveArray[idx].Y<0) {
CurveArray[idx]=new PointF(CurveArray[idx].X,0.0f);
msg = string.Concat(msg,Environment.NewLine,idx.ToString()+"-element was negative");
}
if(CurveArray[idx].Y>g.VisibleClipBounds.Height) {
CurveArray[idx]=new PointF(CurveArray[idx].X,g.VisibleClipBounds.Height);
msg = string.Concat(msg,Environment.NewLine,idx.ToString()+"-element was beyond the visible clip bounds");
}

}
//g.DrawString(msg,new Font("Arial",8),new SolidBrush(Color.Red),10,80);
g.DrawLines(new Pen(Color.Fuchsia), CurveArray);
renderCircules(points);
//points.Clear();
curveRendered = true;
}catch(Exception err) {
g.DrawString("renderCurve: "+err.ToString(),new Font("Arial",10),new SolidBrush(Color.Red),10,5);
}
}

GeneralContactmembericae27 Apr '09 - 5:42 
Hello Fady,
Please write your e-mail.
 
Tnx
GeneralRe: ContactmemberFady Aladdin27 Apr '09 - 7:19 
Hi ... this is my e-mail
fady_aladdin@hotmail.com
GeneralArea under the curvememberDorisvaldo16 Sep '08 - 15:06 
Hi, I'M looking for sample code for calculations of areas under the curve or integrating peaks.
An application on Chromatography.
thanks for any help
Generalpiecewise beziermemberSarath.14 Sep '08 - 21:54 
Can you differentiate it with piecewise bezier. I'm also into curve fitting algo. What's the advantage on using Lagrange Polynomial?
 
BTW, the application is having lot of flickering. Please consider updating the code.
 
-Sarath.
"Great hopes make everything great possible" - Benjamin Franklin

GeneralRe: piecewise beziermembercatbertfromgilead15 Sep '08 - 9:31 
well, bezier doesn't go through all the points that you specify
GeneralI don't mean to quibble but this looks like interpolationmemberMicroImaging11 Sep '08 - 11:15 
As a Numerical Computer Programmer working with Filters, Data Fitting the pictuer and results look like interpolation not fitting.
Fitting can be described by using Least squares fitting, Least Absolute deviation or using Chebyshev fitting
 
Also you do not mention the most important property of orthogonality of the set of Legendre Polynomials.
I encourage you to be more exact in your use of language if you are going to continue your pursuit of Numerical computation, otherwise
you could lose respect if you misuse many different words.
Further Research on the Subject
http://reference.wolfram.com/mathematica/tutorial/CurveFitting.html[^]
 
http://documents.wolfram.com/mathematica/book/section-3.8.2[^]
 
http://citeseer.comp.nus.edu.sg/403088.html[^]
 
In these and many other cases you will see that fitting of data is different from interpolation.
 
Good luck in your career choice, and always remember in computational math that
a + sqrt(b^2-4ac) != a - sqrt(b^2 - 4ac)

GeneralRe: I don't mean to quibble but this looks like interpolationmemberaxelriet12 Sep '08 - 5:30 
Curve Fitting with Splines[^] seems to be a rather valid sentence.
GeneralRe: I don't mean to quibble but this looks like interpolationmemberMicroImaging12 Sep '08 - 9:20 
If you read through this book at

[^]
You will see that data interpolation and data fitting are differentiated, although the title is somewhat ambiguous. The author clearly talks about fitting using norms, especially around page 44 where the chapter title "Curve Fitting: An Introduction" talks specifically as i did abuot the L2, L1 and LInfinity norms.
Again I encourage you to use or learn, or whatever you wish, to differentiate interpolation versus fitting.
GeneralAwesome! Just what I've been looking for!memberdybs6 Sep '08 - 7:42 
Thanks for the article! I'm actually working on a project now where we need to at least fit a quadratic curve, but may in the future need to fit a curve to an arbitrary number of points. We had some difficulty making sense of what we've found so far on the web, but your article REALLY cleared things up. Thanks again!
 
Dybs
GeneralRe: Awesome! Just what I've been looking for!memberFady Aladdin11 Sep '08 - 2:36 
Hi
I am glade that my article clears things up for u ...
Plz tell me more about ur application ..
Regards ..
Fady
GeneralRe: Awesome! Just what I've been looking for!memberdybs11 Sep '08 - 6:32 
I can't give out too much information since it's a work-related project, but we have a graph with 2 fixed points at the end, and we fit a curve to a point the user clicks on somewhere on the graph. Right now we only use a single click point, but we may extend it to work with a series of user-defined points. Sorry, I can't really say much more than that.
 
Thanks again for the article!
 
Dybs
GeneralRe: Awesome! Just what I've been looking for!memberaxelriet12 Sep '08 - 5:26 
You should find plenty of information about interpolation, splines etc in Numerical Analysis[^] books.
GeneralComputational Geometrymemberabombardier2 Sep '08 - 14:37 
Thanks for your article on Lagrange Multipliers. Looking forward for your artical on computationnal geometry !
 
Regards,
André Bombardier
GeneralRe: Computational Geometrymemberaxelriet12 Sep '08 - 5:10 
>Lagrange Multipliers
 
You are referring to another article, right? Wink | ;-)
GeneralRe: Computational Geometrymemberabombardier12 Sep '08 - 9:26 
Yes. It would be nice to see an article on Triangulation. I read in your article that you have interests in Computational Geometry, or you had plans to write an article on it.
 
Have you seen Jonathan Shewchuk's web site ?
 
http://www.cs.berkeley.edu/~jrs/
Generalgetting back valuesmemberKJJ22 Sep '08 - 10:21 
nice work!
 
how do i pass a new value for x and get the y value back from the curve?

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Permalink | Advertise | Privacy | Mobile
Web03 | 2.6.130523.1 | Last Updated 12 Sep 2008
Article Copyright 2008 by Fady Aladdin
Everything else Copyright © CodeProject, 1999-2013
Terms of Use
Layout: fixed | fluid