#include <iostream>
#include <opencv2/opencv.hpp>using namespace cv;template<class T>
struct Point
{T x;T y;Point(){x=0;y=0;}Point(T a,T b){x=a;y=b;}Point operator+(Point p){return Point(x+p.x,y+p.y);}
};
/*拉格朗日曲线,包含多个拉格朗日多项式和曲线的x坐标最大最小值*/
struct LagrangeCurve
{vector<double>factor;vector<double>minuend;int minX;int maxX;
};void calLagrangePolynomial(LagrangeCurve&curve,Point<int>*curvePoints,int pointNum)
{int i,j;curve.maxX=curvePoints[0].x;curve.minX=curvePoints[0].x;for(i=0;i<pointNum;++i){curve.minuend.push_back(curvePoints[i].x);curve.factor.push_back(curvePoints[i].y);//获取曲线两端的x值if(curvePoints[i].x<curve.minX)curve.minX=curvePoints[i].x;if(curvePoints[i].x>curve.maxX)curve.maxX=curvePoints[i].x;//计算(xi-x0)*...(xi-x5),并且把x0...x5保存for(j=0;j<pointNum;++j){if(j!=i){curve.factor[i]/=(curvePoints[i].x-curvePoints[j].x+1e-8);}}}
}double getInterpolation(int x,LagrangeCurve& lagrangeCurve)
{int i,j;double y=0;for(i=0;i<lagrangeCurve.minuend.size();++i){double temp=lagrangeCurve.factor[i];for(j=0;j<lagrangeCurve.minuend.size();++j){if(j!=i){temp*=(x-lagrangeCurve.minuend[j]);}}y+=temp;}return y;
}int main(int argc, const char * argv[])
{Mat img(256,256,CV_8UC3);img=Scalar::all(255);Point<int>*points=new Point<int>[3];points[0].x=10;points[0].y=30;points[1].x=100;points[1].y=70;;points[2].x=250;points[2].y=210;circle(img, cvPoint(points[0].x,points[0].y),5,Scalar::all(0),-1);circle(img, cvPoint(points[1].x,points[1].y),5,Scalar::all(0),-1);circle(img, cvPoint(points[2].x,points[2].y),5,Scalar::all(0),-1);LagrangeCurve lagrangeCurve;calLagrangePolynomial(lagrangeCurve,points,3);for(int i=lagrangeCurve.minX;i<lagrangeCurve.maxX;++i){double y=getInterpolation(i,lagrangeCurve);img.at<Vec3b>(static_cast<int>(y+0.5),i)=Vec3b(0,0,0);}imshow("img", img);waitKey(-1);return 0;
}