圆检测（续）- RANSAC

继之前提到的两种方法之后，这里再列出基于RANSAC的圆检测，RANSAC(Random Sample Consensus)随机抽样一致性，略不同于霍夫圆变换那种基于投票的策略，这是一种对观测数据进行最大化模型检验的方法。下面来简单介绍一下它的原理：
1、原理
最小二乘法通常用在线性拟合参数中，但一旦最小二乘法输入的观测数据中包含有大量分散的干扰点时，它拟合出来的效果可能并不好，如可能会出现这样的情况：

可以看到，拟合出来的直线与期望有效点之间的重合率不高，也就代表着它的代价函数Cost(m,b)=∑ni=1|yi−(mxi+b)|虽然已经是最小的了，期望(不是概率论里面的期望)函数的值却不是最大的。
Ransac的思路是随机通过几个点用最小二乘法给出一个假设的直线，然后计算在直线内的inliers和在直线范围外的outliers。对所有可能的直线中找出inliers数目最多的那个，也就能找到最好的直线。
算法步骤：
（1） 随机地抽取出所需要数目的点去拟合模型：

绿色的点代表采样的点
（2）用样本求出模型参数：

（3）在设定好直线的阈值范围中区分出inliers和outliers，并求内点占观测数据的比：

重复步骤1-3直到找出置信度最高的模型
RANSAC几点要注意的：
① 只有outliers%<50%时，得到的模型才是有保证的。
②inliers的阈值δ跟我们期望的模型抗噪能力相关，阈值越大，抗噪能力越弱，通常我们选取3个像素偏差的高斯模型作为噪声模型。
③ 重复1-3步骤的次数跟模型的outliers占比和我们所需要多高的置信度有关，可以用下面公式表示：

S=log(1−P)log(1−pk)

其中，S是所需最小试验的次数，P是置信度，p是inliers占的百分比数，k是随机采样的数目。

2、实际例子
这里用网友Micka的代码来举例：

#include <opencv2/opencv.hpp>#include <ctime>float verifyCircle(cv::Mat dt, cv::Point2f center, float radius, std::vector<cv::Point2f> & inlierSet){    unsigned int counter = 0;    unsigned int inlier = 0;    float minInlierDist = 2.0f;    float maxInlierDistMax = 100.0f;    float maxInlierDist = radius/25.0f;    if(maxInlierDist<minInlierDist) maxInlierDist = minInlierDist;    if(maxInlierDist>maxInlierDistMax) maxInlierDist = maxInlierDistMax;    // choose samples along the circle and count inlier percentage    for(float t =0; t<2*3.14159265359f; t+= 0.05f)    {        counter++;        float cX = radius*cos(t) + center.x;        float cY = radius*sin(t) + center.y;        if(cX < dt.cols)            if(cX >= 0)                if(cY < dt.rows)                    if(cY >= 0)                        if(dt.at<float>(cY,cX) < maxInlierDist)                        {                            inlier++;                            inlierSet.push_back(cv::Point2f(cX,cY));                        }    }    return (float)inlier/float(counter);}inline void getCircle(cv::Point2f& p1,cv::Point2f& p2,cv::Point2f& p3, cv::Point2f& center, float& radius){    float x1 = p1.x;    float x2 = p2.x;    float x3 = p3.x;    float y1 = p1.y;    float y2 = p2.y;    float y3 = p3.y;    // PLEASE CHECK FOR TYPOS IN THE FORMULA :)    center.x = (x1*x1+y1*y1)*(y2-y3) + (x2*x2+y2*y2)*(y3-y1) + (x3*x3+y3*y3)*(y1-y2);    center.x /= ( 2*(x1*(y2-y3) - y1*(x2-x3) + x2*y3 - x3*y2) );    center.y = (x1*x1 + y1*y1)*(x3-x2) + (x2*x2+y2*y2)*(x1-x3) + (x3*x3 + y3*y3)*(x2-x1);    center.y /= ( 2*(x1*(y2-y3) - y1*(x2-x3) + x2*y3 - x3*y2) );    radius = sqrt((center.x-x1)*(center.x-x1) + (center.y-y1)*(center.y-y1));}std::vector<cv::Point2f> getPointPositions(cv::Mat binaryImage){    std::vector<cv::Point2f> pointPositions;    for(unsigned int y=0; y<binaryImage.rows; ++y)    {        //unsigned char* rowPtr = binaryImage.ptr<unsigned char>(y);        for(unsigned int x=0; x<binaryImage.cols; ++x)        {            //if(rowPtr[x] > 0) pointPositions.push_back(cv::Point2i(x,y));            if(binaryImage.at<unsigned char>(y,x) > 0) pointPositions.push_back(cv::Point2f(x,y));        }    }    return pointPositions;}int main(){    clock_t starttime, endtime;    starttime = clock();    cv::Mat color = cv::imread("1.jpg");    cv::Mat gray;    // convert to grayscale    // you could load as grayscale if you want, but I used it for (colored) output too    cv::cvtColor(color, gray, CV_BGR2GRAY);    cv::Mat mask;    float canny1 = 100;    float canny2 = 20;    cv::Mat canny;    cv::Canny(gray, canny, canny1,canny2);    //cv::imshow("canny",canny);    mask = canny;    std::vector<cv::Point2f> edgePositions;    edgePositions = getPointPositions(mask);    // create distance transform to efficiently evaluate distance to nearest edge    cv::Mat dt;    cv::distanceTransform(255-mask, dt,CV_DIST_L1, 3);    //TODO: maybe seed random variable for real random numbers.    unsigned int nIterations = 0;    cv::Point2f bestCircleCenter;    float bestCircleRadius;    float bestCirclePercentage = 0;    float minRadius = 10;   // TODO: ADJUST THIS PARAMETER TO YOUR NEEDS, otherwise smaller circles wont be detected or "small noise circles" will have a high percentage of completion    //float minCirclePercentage = 0.2f;    float minCirclePercentage = 0.05f;  // at least 5% of a circle must be present? maybe more...    int maxNrOfIterations = edgePositions.size();   // TODO: adjust this parameter or include some real ransac criteria with inlier/outlier percentages to decide when to stop    printf("%d\n", maxNrOfIterations);    for(unsigned int its=0; its< maxNrOfIterations; ++its)    {        //RANSAC: randomly choose 3 point and create a circle:        //TODO: choose randomly but more intelligent,         //so that it is more likely to choose three points of a circle.         //For example if there are many small circles, it is unlikely to randomly choose 3 points of the same circle.        unsigned int idx1 = rand()%edgePositions.size();        unsigned int idx2 = rand()%edgePositions.size();        unsigned int idx3 = rand()%edgePositions.size();        // we need 3 different samples:        if(idx1 == idx2) continue;        if(idx1 == idx3) continue;        if(idx3 == idx2) continue;        // create circle from 3 points:        cv::Point2f center; float radius;        getCircle(edgePositions[idx1],edgePositions[idx2],edgePositions[idx3],center,radius);        // inlier set unused at the moment but could be used to approximate a (more robust) circle from alle inlier        std::vector<cv::Point2f> inlierSet;        //verify or falsify the circle by inlier counting:        float cPerc = verifyCircle(dt,center,radius, inlierSet);        // update best circle information if necessary        if(cPerc >= bestCirclePercentage)            if(radius >= minRadius)            {                bestCirclePercentage = cPerc;                bestCircleRadius = radius;                bestCircleCenter = center;            }    }    std::cout << "bestCirclePerc: " << bestCirclePercentage << std::endl;    std::cout << "bestCircleRadius: " << bestCircleRadius << std::endl;    // draw if good circle was found    if(bestCirclePercentage >= minCirclePercentage)        if(bestCircleRadius >= minRadius);    cv::circle(color, bestCircleCenter,bestCircleRadius, cv::Scalar(255,255,0),1);    std::cout << "the used time is: "<<clock()-starttime <<std::endl;    cv::imshow("output",color);    cv::imshow("mask",mask);    //cv::imwrite("../outputData/1_circle_normalized.png", normalized);    cv::waitKey(0);    return 0;}

他的思路是：
1. 用Canny提取边缘点， 用distanceTransform得到距离边缘点的距离图；
2. 随机抽取三个不同的点解方程，三个方程三个未知数，有解；
3. 将2得到的圆周上的点与1中对应位置的点进行比较，看是否属于inliers，随后输出百分比；
4. 找出最大百分比对应的圆就是RANSAC得到的圆。

3、比较霍夫变换跟RANSAC：
鲁棒性来说，霍夫变换要稳定一点；
速度来说，霍夫变换要快，而且其所需时间变化不大，100ms左右能够完成；
RANSAC跟HoughTranform的参数调节都很麻烦，相对来说，霍夫变换更加简单一点；
RANSAC拟合的程度可能会更高，但受到outliers%<50%这个条件限制。
所以综合来说，HoughTransform的应用更广，效率更高，某些情况下，它不能很好地找到合理的圆，这时可以将RANSAC加进去进行优化，可能精度会高很多。

参考资料： RANSAC Kavita Bala