先上Github地址

線性擬合的應(yīng)用領(lǐng)域比較廣泛, 如運(yùn)動(dòng)軌跡計(jì)算、數(shù)據(jù)分析、圖像處理等領(lǐng)域, 故在此寫一篇學(xué)習(xí)相關(guān)文章.

直線擬合

在說(shuō)直線擬合前我們先來(lái)復(fù)習(xí)一下直線方程. 直線方程的表達(dá)式有以下幾種(引用自百度):

1. 一般式:Ax+By+C=0(A、B不同時(shí)為0)【適用于所有直線】
   A1/A2=B1/B2≠C1/C2←→兩直線平行
   A1/A2=B1/B2=C1/C2←→兩直線重合
   橫截距a=-C/A
   縱截距b=-C/B
2. 點(diǎn)斜式:y-y0=k(x-x0) 【適用于不垂直于x軸的直線】
   表示斜率為k，且過(guò)（x0,y0）的直線
3. 截距式:x/a+y/b=1【適用于不過(guò)原點(diǎn)或不垂直于x軸、y軸的直線】
   表示與x軸、y軸相交，且x軸截距為a，y軸截距為b的直線
4. 斜截式:y=kx+b【適用于不垂直于x軸的直線】
   表示斜率為k且y軸截距為b的直線
5. 兩點(diǎn)式:【適用于不垂直于x軸、y軸的直線】
   表示過(guò)（x1,y1）和(x2,y2)的直線
   **(y-y1)/(y2-y1)=(x-x1)/(x2-x1) **(****x1≠x2，y1≠y2****)****
6. 交點(diǎn)式:f1(x,y) *m+f2(x,y)=0 【適用于任何直線】
   表示過(guò)直線f1(x,y)=0與直線f2(x,y)=0的交點(diǎn)的直線
7. 點(diǎn)平式:f(x,y) -f(x0,y0)=0【適用于任何直線】
   表示過(guò)點(diǎn)（x0,y0）且與直線f（x,y）=0平行的直線
8. 法線式：x·cosα+ysinα-p=0【適用于不平行于坐標(biāo)軸的直線】
   過(guò)原點(diǎn)向直線做一條的垂線段，該垂線段所在直線的傾斜角為α，p是該線段的長(zhǎng)度
9. 點(diǎn)向式：(x-x0)/u=(y-y0)/v (u≠0,v≠0)【適用于任何直線】
   表示過(guò)點(diǎn)(x0,y0)且方向向量為（u,v ）的直線
10. 法向式：a（x-x0）+b（y-y0）=0【適用于任何直線】****
    表示過(guò)點(diǎn)（x0，y0）且與向量（a，b）垂直的直線

在OpenCV中使用的是點(diǎn)向式：(x-x0)/u=(y-y0)/v (u≠0,v≠0)

在OpenCV中我們可以主要使用fitLine函數(shù)進(jìn)行直線擬合, 其算法基于最小二乘法實(shí)現(xiàn)的, 關(guān)于最小二乘法后續(xù)再作詳細(xì)介紹. fitLine函數(shù)定義如下:

void fitLine( InputArray points, OutputArray line, int distType, double param, double reps, double aeps );

該方法可以對(duì)二維或三維的點(diǎn)集進(jìn)行擬合, 在二維模式下返回一個(gè)Vec4f的元素, 分別是點(diǎn)向式中的(u,v,x,y)

參數(shù)介紹:
points
輸入一個(gè)二維或三維的向量, 數(shù)據(jù)結(jié)構(gòu)可以是std::vector<> 或 Mat(推薦使用).

line
二維模式下返回Vec4f的元素(u,v,x,y), 三維模式下返回Vec6f的元素(vx, vy, vz, x0, y0, z0)

distType
擬合算法, 這里的擬合算法基于M-estimators實(shí)現(xiàn), 有以下幾種:

enum DistanceTypes {
    DIST_USER    = -1,  //!< User defined distance
    DIST_L1      = 1,   //!< distance = |x1-x2| + |y1-y2|
    DIST_L2      = 2,   //!< the simple euclidean distance
    DIST_C       = 3,   //!< distance = max(|x1-x2|,|y1-y2|)
    DIST_L12     = 4,   //!< L1-L2 metric: distance = 2(sqrt(1+x*x/2) - 1))
    DIST_FAIR    = 5,   //!< distance = c^2(|x|/c-log(1+|x|/c)), c = 1.3998
    DIST_WELSCH  = 6,   //!< distance = c^2/2(1-exp(-(x/c)^2)), c = 2.9846
    DIST_HUBER   = 7    //!< distance = |x|<c ? x^2/2 : c(|x|-c/2), c=1.345
};

其算法主要使用的是加權(quán)最小二乘法, 分別采用如下算法:

DistanceTypes

其中CV_DIST_L2為最簡(jiǎn)單快速的最小二乘法, 推薦使用.

param
擬合算法中的的C參, 設(shè)置為0時(shí)系統(tǒng)自動(dòng)選用最優(yōu)值.

reps
算法精度radius, 官方推薦0.01

aeps
算法精度angle, 官方推薦0.01

曲線擬合

首先，曲線的函數(shù)類型有:雙曲線型、對(duì)數(shù)型、指數(shù)型、多項(xiàng)式型等。

對(duì)于最常見的曲線, 其特征接近于多項(xiàng)式型, 所以這里的曲線擬合問(wèn)題就變成了多項(xiàng)式回歸求解.
一般情況下,對(duì)實(shí)際的點(diǎn)數(shù)據(jù)進(jìn)行求解是無(wú)解的, 所以我們需要引入最小二乘法來(lái)求解.(對(duì)于最小二乘法的原理和求解過(guò)程,后續(xù)再詳細(xì)介紹)

在OpenCV中沒(méi)有找到擬合曲線的算子, 這里使用C++來(lái)實(shí)現(xiàn). 此算法迭代次數(shù)大于10時(shí)求解均為0. 若有更優(yōu)的算法可以在評(píng)論中交流.

bool fittingCurve(vector<Point2f> &vec,int times,float *p) // 輸入點(diǎn)數(shù)據(jù), 迭代次數(shù), 輸出多項(xiàng)式參數(shù)
{
    float *py = new float[vec.size()];
    float *px = new float[times*vec.size()];
    int i = 0;
    
    for(vector<Point2f>::iterator itr = vec.begin();itr!=vec.end();++itr) {
        
        py[i] = (*itr).y;
        int j=0;
        while (j<times)
        {
            px[times*i+j] = pow(((*itr).x),float(j));
            j++;
        }
        i++;
    }
    
    Mat X = Mat((int)vec.size(),times,CV_32FC1,px);
    Mat X_T;
    transpose(X,X_T);
    
    Mat Y = Mat((int)vec.size(),1,CV_32FC1,py);
    Mat para = ((X_T*X).inv())*X_T*Y;
    
    for (int s = 0; s<times;s++){
        p[s] = para.at<float>(s);
    }
    
    delete[] px;
    delete[] py;
    
    return true;
}

下面我們來(lái)驗(yàn)證一組數(shù)據(jù):

        int x[]={50,100,150,200,250,300,350,400,450,500,550,600,650,700,750};
        int y[]={428,454,480,506,532,458,384,210,636,662,688,778,504,430,456};

當(dāng)?shù)螖?shù)為3的時(shí)候求解結(jié)果為:

381.534
0.575073
-0.000505627

即:y = 381.534 + 0.575073x + -0.000505627x2, 圖像如下:

當(dāng)?shù)螖?shù)為7的時(shí)候求解結(jié)果為:

-210.757
17.1697
-0.132375
0.000430976
-6.52814e-07
4.45947e-10
-1.06088e-13

圖像如下:

以下是一些自己測(cè)試的數(shù)據(jù):