opencv的曲線擬合polyfit

推薦一個不錯的網頁，可以直接用solve函數求解方程組：

http://m.blog.csdn.net/u014652390/article/details/52789591

 

4.1 曲線擬合的最小二乘法

求以下擬合函數

擬合條件：擬合曲線與各數據點在y方向的誤差平方和最小.

擬合函數為一元函數時--函數圖形為平面曲線--曲線擬合

解決曲線擬合，最先是確定擬合函數的形式。即適當選取

 

選冪函數{1,x,x², ···,xⁿ}, 則多項式擬合函數φ(x)可表示為：

φ(x)=a₀+a₁*x+a₂*x²+a₃*x³+......+an*xⁿ =[a₀ a₁ a₂ ...... a_n][1 x₁ x₁² ... ... x₁ⁿ]^T       (n+1<m)

a0、a1、a2......an是冪系數，也是擬合所求的未知量。

實際中擬合函數有指數函數、三角函數等，根據數據 的分布特點來選取合適的擬合函數。

將第 i 個樣本點的x坐標帶入φ(x)，得到：

這個就是二次方程，我們期望S最小。此時，方程中的x、y已知，想求的是a₀ a₁ a₂ ...... a_n。

S最小的必要條件是：

整理得到如下正規方程組：

解此方程組得系數a₀ a₁ a₂ ...... a_n，, 得出擬合函數φ(x)

最小二乘法：以殘差平方和最小問題的解來確定擬合函數

二、超定方程組得最小二乘解

 

將

寫成向量內積形式：

a₀ a₁ a₂ ...... a_n為待定系數，滿足：

此m個等式如下建立方程組：

方程數(m)多於未知數個數(n+1)，此類方程組稱為超定方程組。下列正規方程組中k個方程中a_j的系數

 

 

經推導，得到最小二次方，冪函數擬合公式如下：

 Φ^T* Φ*a= Φ^T*y

 其中Φ是樣本點坐標x的超定矩陣，將所有x帶入該向量[1  x  x^2 ... ...  x^n]中，就得到超定矩陣Φ。Φ^T表示Φ的轉置

 

#include <iostream>
#include<opencv2/opencv.hpp>
using namespace std;
using namespace cv;
//下面宏定義CV_MAT_ELEM2為方便快速訪問圖像像素
#define CV_MAT_ELEM2(src,dtype,y,x) \
        (dtype*)(src.data+src.step[0]*y+src.step[1]*x)

Mat polyfit(std::vector<cv::Point2f> &chain,int n)
{
    Mat y(chain.size(),1,CV_32F,Scalar::all(0));
/* ********【預聲明phy超定矩陣】************************/
/* 多項式擬合的函數為多項冪函數
 * f(x)=a0+a1*x+a2*x^2+a3*x^3+......+an*x^n
 *a0、a1、a2......an是冪系數，也是擬合所求的未知量。設有m個抽樣點，則：
 * 超定矩陣phy=1 x1 x1^2 ... ...  x1^n
 *           1 x2 x2^2 ... ...  x2^n
 *           1 x3 x3^2 ... ...  x3^n
 *              ... ... ... ...
 *              ... ... ... ...
 *           1 xm xm^2 ... ...  xm^n
 *
 * *************************************************/
    cv::Mat phy(chain.size(),n,CV_32F,Scalar::all(0));
    for(int i=0;i<phy.rows;i++)
    {
        float* pr=phy.ptr<float>(i);
        for(int j=0;j<phy.cols;j++)
        {
           pr[j]=pow(chain[i].x,j);
        }
        y.at<float>(i)=chain[i].y;
    }
    Mat phy_t=phy.t();
    Mat phyMULphy_t=phy.t()*phy;
    Mat phyMphyInv=phyMULphy_t.inv();
    Mat a=phyMphyInv*phy_t;
    a=a*y;
    return a;
}

int main()
{
    vector<Point2f> sp;
    //設有二次曲線點 g(x)=5+2.6x+2x^3,則:
    float a[]={5,2.6,2};
    Mat image(500,500,CV_32FC1,Scalar(0));
    RNG rng;//預聲明一個隨機變量，用於作為離散點的干擾項
    for(int i=1;i<20;i+=2)
    {
        Point2f p;
        p.x=i;
        for(int k=0;k<sizeof(a);k++)
        {
            p.y +=a[k]*pow(i,k);//
        }

        p.y +=rng.uniform(-1,1);//為理想點位置添加隨機干擾
     /*將上面的p點以圓點的形式繪制到image上，為了觀察方便，
      * 將y坐標做了顛倒，坐標原點在image的左下角*/
        Point2f pi;
        pi.x=p.x;
        pi.y=image.rows-p.y;
        circle(image,pi,3,Scalar(255),-1);
      /*-------------end--------------------*/
        sp.push_back(p);
        cout<<p<<endl;
    }
    image.convertTo(image,CV_8UC1);
    imshow("distributed Points",image);
    Mat am=polyfit(sp,3);
    cout<<am<<endl;
    waitKey();
    return 0;
}