ransac 算法原理介紹與實現

ransac是RANdom SAmple Consensus的簡稱，它是根據一組包含異常數據的樣本數據集，經過迭代的方法計算出數據的數學模型參數，獲得有效樣本數據的非肯定性的算法。它於1981年由 Fischler和Bolles最早提出。

對於RANSAC算法有一個基本的假設：樣本中包含正確數據(inliers，符合模型的數據)和異常數據(Outliers，不符合模型的數據)，即數據集中含有噪聲。這些異常數據多是因爲錯誤的測量、錯誤的假設、錯誤的計算等產生的。同時RANSAC也假設， 給定一組正確的數據，存在能夠計算出符合這些數據的模型參數的方法。

下面是WIKI上一個，也是一個比較經典的例子：平面直線的匹配。

左邊是數據集，右邊是經過RANSAC算法擬合的直線。

利用C實現

void Ransac(IplImage *pload,IplImage *pnew,int *x_cord,int *y_cord,int length, double &c0,double &c1,double &c2,char *file_name_prefix2,int nums,int &min)
{
    int *count = new int[RANSAC_TIMES]; //計數
    memset(count,0x00,sizeof(int)*RANSAC_TIMES);
    CvMat  mat_x, mat_b, mat_c;
    double  p_mat_x[9], p_mat_b[3],p_mat_c[3];//要轉化成數組
    min = length; //求最小值用
    int flag_linear =0;         
    int position = 0; //記錄偏差最小的位置

    uchar *loadimageData = (uchar*)pload->imageData;
    int step = pload->widthStep;
    double a =0,b=0,c =0;
    for(int i =0; i<RANSAC_TIMES; i++)
    {
        int randnums[3];
        int n = 0;

        for (int t=0;t<3;t++)
        {
            n = rand() % length;
            randnums[t] = n;
            p_mat_x[t*3] = 1.0;
            p_mat_x[t*3+1] = x_cord[n];
            p_mat_x[t*3+2] = x_cord[n] * x_cord[n];

            p_mat_b[t] = y_cord[n];
            p_mat_c[t] = 0;
        }
        cvInitMatHeader(&mat_x,3,3,CV_64FC1,p_mat_x);
        cvInitMatHeader(&mat_b,3,1,CV_64FC1,p_mat_b);
        cvInitMatHeader(&mat_c,3,1,CV_64FC1,p_mat_c);

        flag_linear = cvSolve(&mat_x, &mat_b, &mat_c,CV_LU);
        if ( flag_linear == 0)
        {
            continue;
        }
        double *temp = mat_c.data.db;  //結果保存下來。
        c = temp[0];   //常數項   //保留第i次結果
        b = temp[1];   //一次項
        a = temp[2];   //平方項
        if (c2 < 0)
        {
            c2 = -c2;
        }
        double y_value = 0;
        double x_square = 0;
        int y_floor = 0;

        double dis_error = 0;

        for (int j=0;j<length;j++)
        {
            dis_error = fabs(y_cord[j] - ( a*(x_cord[j]*x_cord[j]) + b*x_cord[j] + c )); 
            if ( dis_error > RANSAC_THRESHOLD) 
            {
                count[i]++;
            }
        }

        if (min > count[i])
        {
            min = count[i];
            position = i;
            c0 = c;
            c1 = b;
            c2 = a;
        }   

    }   //50次循環結束。

    double x_square = 0;
    double y_value = 0;
    int y_floor = 0;
    int pixelposition;
    for (int x=0; x<pnew->width; x++)   //自變量取值範圍
    {
        x_square = pow(double(x),2);
        y_value = c2*x_square + c1*x + c0;
        y_floor = cvFloor((y_value));

        if ( y_floor >= 0 && y_floor < pnew->height )
        {
            ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x]=255;
            ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x+1]=0;
            ((uchar*)(pnew->imageData + pnew->widthStep*y_floor))[3*x+2]=0;
        }
        for (int y=0;y<pnew->height;y++)
        {
            pixelposition = y*step+x;
            if ( 0 == loadimageData[pixelposition])
            {
                ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x]= 0;
                ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x+1]=0;
                ((uchar*)(pnew->imageData + pnew->widthStep*y))[3*x+2]=255;
            }
        }
    }

    delete []count;
    cvSaveImage(file_name_prefix2,pnew);

}