拉格朗日插值編程實現

拉格朗日插值原理：

拉格朗日插值的具體介紹網址：https://zh.wikipedia.org/wiki/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E6%8F%92%E5%80%BC%E6%B3%95

翻譯成人話就是，該曲線是由多個n次多項式的和構成的，n是參與插值的點的個數。每個n次多項式的計算方法如上圖所示。轉化成程序的話，就是要保存每個多項式中（x－xi）中的每一項xi。然后就是系數yi/(xj-x0)(xj-x1)...(xj-xk)

程序實現：

 

#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;
}

 

上述代碼使用到了opencv畫圖，運行結果圖如下：