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);

}