Opencv--minEnclosingCircle源码--求最小包围圆的算法

在Opencv中，求最小包围圆的函数：

他的方法是什么呢？怎么找出这个最小包围圆的呢？

为了方便理解源码，我新建了一个工程，把相关的源文件复制过来，再输入一张测试图，让代码跑起来看中间过程。

工程文件包括：

其中 my_minEnclosingCircle是main函数入口，其他文件从C:\...\opencv\sources\modules\... 路径下搜索，然后复制到该工程文件夹内。这样在代码中使用“转到定义”的时候就可以查看相关源码了，还可以用不同的测试图去测试。

main函数：

#include <opencv2/opencv.hpp>using namespace cv;using namespace std;int main(){Mat gray = imread("3.jpg",0);Mat  src;threshold(gray, src, 200, 255, THRESH_BINARY);imshow("src", src);vector< vector<Point> > contours;vector<Vec4i> hiearachy;Mat src_cp = src.clone();findContours(src_cp, contours, hiearachy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE);Point2f center;float radius;minEnclosingCircle(contours[0], center, radius);Mat dst(src.size(), CV_8UC3, Scalar::all(0));Point center1;center1.x  = cvRound (center.x);center1.y = cvRound(center.y);int radius1 = cvRound(radius);circle(dst, center1, radius1, Scalar(0,0,255), 2, 8);dst.setTo(Scalar(255,255,255), src);imwrite("dst.jpg", dst);waitKey(0);return 0;}

测试图（左）、自己做一些标记（中）、最小包围圆（右）

   

minEnclosingCircle 函数定义：

void cv::minEnclosingCircle( InputArray _points,Point2f& center, float& radius ){Mat points = _points.getMat();CV_Assert(points.checkVector(2) >= 0 && (points.depth() == CV_32F || points.depth() == CV_32S));CvMat _cpoints = points;cvMinEnclosingCircle( &_cpoints, (CvPoint2D32f*)¢er, &radius );}

c++的接口，然后转到c函数中去。cvMinEnclosingCircle函数定义：

CV_IMPL intcvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius ){    const int max_iters = 100;    const float eps = FLT_EPSILON*2;    CvPoint2D32f center = { 0, 0 };    float radius = 0;    int result = 0;    if( _center )        _center->x = _center->y = 0.f;    if( _radius )        *_radius = 0;    CvSeqReader reader;    int k, count;    CvPoint2D32f pts[8];    CvContour contour_header;    CvSeqBlock block;    CvSeq* sequence = 0;    int is_float;    if( !_center || !_radius )        CV_Error( CV_StsNullPtr, "Null center or radius pointers" );    if( CV_IS_SEQ(array) )    {        sequence = (CvSeq*)array;        if( !CV_IS_SEQ_POINT_SET( sequence ))            CV_Error( CV_StsBadArg, "The passed sequence is not a valid contour" );    }    else    {        sequence = cvPointSeqFromMat(            CV_SEQ_KIND_GENERIC, array, &contour_header, &block );    }    if( sequence->total <= 0 )        CV_Error( CV_StsBadSize, "" );    cvStartReadSeq( sequence, &reader, 0 );    count = sequence->total;    is_float = CV_SEQ_ELTYPE(sequence) == CV_32FC2;    if( !is_float )    {        CvPoint *pt_left, *pt_right, *pt_top, *pt_bottom;        CvPoint pt;        pt_left = pt_right = pt_top = pt_bottom = (CvPoint *)(reader.ptr);        CV_READ_SEQ_ELEM( pt, reader );        for(int i = 1; i < count; i++ )        {            CvPoint* pt_ptr = (CvPoint*)reader.ptr;            CV_READ_SEQ_ELEM( pt, reader );            if( pt.x < pt_left->x )                pt_left = pt_ptr;            if( pt.x > pt_right->x )                pt_right = pt_ptr;            if( pt.y < pt_top->y )                pt_top = pt_ptr;            if( pt.y > pt_bottom->y )                pt_bottom = pt_ptr;        }        pts[0] = cvPointTo32f( *pt_left );        pts[1] = cvPointTo32f( *pt_right );        pts[2] = cvPointTo32f( *pt_top );        pts[3] = cvPointTo32f( *pt_bottom );    }    else    {        CvPoint2D32f *pt_left, *pt_right, *pt_top, *pt_bottom;        CvPoint2D32f pt;        pt_left = pt_right = pt_top = pt_bottom = (CvPoint2D32f *) (reader.ptr);        CV_READ_SEQ_ELEM( pt, reader );        for(int i = 1; i < count; i++ )        {            CvPoint2D32f* pt_ptr = (CvPoint2D32f*)reader.ptr;            CV_READ_SEQ_ELEM( pt, reader );            if( pt.x < pt_left->x )                pt_left = pt_ptr;            if( pt.x > pt_right->x )                pt_right = pt_ptr;            if( pt.y < pt_top->y )                pt_top = pt_ptr;            if( pt.y > pt_bottom->y )                pt_bottom = pt_ptr;        }        pts[0] = *pt_left;        pts[1] = *pt_right;        pts[2] = *pt_top;        pts[3] = *pt_bottom;    }    for( k = 0; k < max_iters; k++ )    {        double min_delta = 0, delta;        CvPoint2D32f ptfl, farAway = { 0, 0};        /*only for first iteration because the alg is repared at the loop's foot*/        if(k==0)            icvFindEnslosingCicle4pts_32f( pts, ¢er, &radius );        cvStartReadSeq( sequence, &reader, 0 );        for(int i = 0; i < count; i++ )        {            if( !is_float )            {                ptfl.x = (float)((CvPoint*)reader.ptr)->x;                ptfl.y = (float)((CvPoint*)reader.ptr)->y;            }            else            {                ptfl = *(CvPoint2D32f*)reader.ptr;            }            CV_NEXT_SEQ_ELEM( sequence->elem_size, reader );            delta = icvIsPtInCircle( ptfl, center, radius );             if( delta < min_delta )            {                min_delta = delta;                farAway = ptfl;            }        }        result = min_delta >= 0;        if( result )            break;        CvPoint2D32f ptsCopy[4];        /* find good replacement partner for the point which is at most far away,        starting with the one that lays in the actual circle (i=3) */        for(int i = 3; i >=0; i-- )        {            for(int j = 0; j < 4; j++ )            {                ptsCopy[j]=(i != j)? pts[j]: farAway;            }            icvFindEnslosingCicle4pts_32f(ptsCopy, ¢er, &radius );            if( icvIsPtInCircle( pts[i], center, radius )>=0){ // replaced one again in the new circle?                pts[i] = farAway;                break;            }        }    }    if( !result )    {        cvStartReadSeq( sequence, &reader, 0 );        radius = 0.f;        for(int i = 0; i < count; i++ )        {            CvPoint2D32f ptfl;            float t, dx, dy;            if( !is_float )            {                ptfl.x = (float)((CvPoint*)reader.ptr)->x;                ptfl.y = (float)((CvPoint*)reader.ptr)->y;            }            else            {                ptfl = *(CvPoint2D32f*)reader.ptr;            }            CV_NEXT_SEQ_ELEM( sequence->elem_size, reader );            dx = center.x - ptfl.x;            dy = center.y - ptfl.y;            t = dx*dx + dy*dy;            radius = MAX(radius,t);        }        radius = (float)(sqrt(radius)*(1 + eps));        result = 1;    }    *_center = center;    *_radius = radius;    return result;}

输入一个点集，输出圆心和半径。总结算法的几个步骤：

• 1.遍历所有点，找出最左边、最右边、最上边、最下边的四个点，分别用A、B、C、D来表示。
• 2.求出包围这四个点的最小圆C1的圆心和半径；
• 3.第一次迭代（k=0）。遍历所有点，检查是否存在点出界，即不在圆C1界内也不在边界上。如果没有出界点，则求出的圆就是最终所求。如果有出界点，转到第四步。
• 4.假设出界点中距离圆C1的圆心最远的点为E，则依次尝试以下四种组合：(1)A/B/C/E  (2)A/B/D/E  (3)A/C/D/E  (4)B/C/D/E。求出该组合中四个点的最小包围圆。假设组合(1)A、B、C、E的最小包围圆是C2，然后检测剩下的点D在不在C2内，比如说不在，那就算组合(2)，算出A、B、D、E四个点的最小包围圆C3，检测剩下的点C，发现在圆C3内，那么记录下点E和圆C3的圆心、半径。
• 5.第二次迭代（k=1）。遍历所有点，检查是否存在出界点，即不在圆C3界内也不在边界上。如果没有出界点，则求出的圆C3就是最终所求。如果有出界点，则转到下一步。
• 6.假设出界点中距离圆C3的圆心最远的点为F，则依次尝试以下四种组合：(1)A/B/D/F  (2)A/B/E/F  (3)A/D/E/F  (4)B/D/E/F  (这里已经没有C什么事了，它已经在上一轮被淘汰了)。依次处理这四个组合，求出该组合中四个点的最小包围圆。假设组合(1)求出圆C4，检测点E在圆C4内，则么记录下点F和圆C4的圆心、半径。
• 7.第三次迭代（k=2）。接下来重复第5步、第6步，直到发现遍历完所有点结果都在最新求出的圆界内，则该圆就是最终所求的圆，退出迭代。

在cvMinEnclosingCircle 函数中，pts[0]、pts[1] 、pts[2] 、pts[3] 分别用来表示最左边、最右边、最上边、最下边的四个点。

最大迭代次数设置为100.  for( k = 0; k < max_iters; k++ )部分一个大括号中完成一次迭代。

icvFindEnslosingCicle4pts_32f( pts, &center, &radius ); 该函数是用来求四个点的最小包围圆。从数学原理上，三个不在一条直线上的点确定一个圆。那四个点的最小包围圆怎么求呢？是否是唯一的呢？参考网上的回答：

根据代码中的办法，首先找到距离最远的两个点，比如B和D，保存在idxs[0]、idxs[1]中，然后求出这两个点的中心P，距离的一半为半径，画出圆C11。然后判断剩下的两个点A和C是否过界，如果不过界则圆C11为最终所求的圆，如果有过界的点，则取四个点中的三个画圆，判断第四个点是否在圆内。四种组合中，如果画出了不止一个圆，则比较半径，取半径最小的那个圆作为最终所求的圆。

贴出icvFindEnslosingCicle4pts_32f( pts, &center, &radius ) 函数代码：

static inticvFindEnslosingCicle4pts_32f( CvPoint2D32f * pts, CvPoint2D32f * _center, float *_radius ){    //在理解代码的过程，假设 pts[0]和pts[2]之间的距离最远，即左边和上面的两个点距离最远，右边点的距离比较远。int shuffles[4][4] = { {0, 1, 2, 3}, {0, 1, 3, 2}, {2, 3, 0, 1}, {2, 3, 1, 0} };    int idxs[4] = { 0, 1, 2, 3 };  //idxs[4]里面装了四个点的序号，初值给一任意值。    int i, j, k = 1, mi = 0;    float max_dist = 0;    CvPoint2D32f center; //初始值为pts[0], 后面计算 最远的两个点的中点，保存在这。    CvPoint2D32f min_center;  //初始值为pts[0],     float radius, min_radius = FLT_MAX;    CvPoint2D32f res_pts[4];    center = min_center = pts[0];    radius = 1.f;//两个for循环为了找到距离最远的两个点，更新idxs[0]、idxs[1]的值为这两个点的序号。    for( i = 0; i < 4; i++ )        for( j = i + 1; j < 4; j++ )        {            float dist = icvDistanceL2_32f( pts[i], pts[j] );            if( max_dist < dist )            {                max_dist = dist;                idxs[0] = i;                idxs[1] = j;            }        }    if( max_dist == 0 )        goto function_exit;// idxs[0]、idxs[1]已经装了最远的两个点的序号；排除这两个序号，把下一个序号装在idxs[2]里面；排除以上三个序号，把下一个序号装在idxs[3]里面。    k = 2;    for( i = 0; i < 4; i++ )    {        for( j = 0; j < k; j++ )            if( i == idxs[j] )                break;        if( j == k )            idxs[k++] = i;    }    center = cvPoint2D32f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f,                           (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f );    radius = (float)(icvDistanceL2_32f( pts[idxs[0]], center )*1.03);  //最远距离的一半，再乘以1.03.    if( radius < 1.f )        radius = 1.f;    if( icvIsPtInCircle( pts[idxs[2]], center, radius ) >= 0 && //  if(   pts[idxs[2]]与center这两个点之间的距离小于radius )        icvIsPtInCircle( pts[idxs[3]], center, radius ) >= 0 )  // if(   pts[idxs[3]]与center这两个点之间的距离小于radius )    {        k = 2; //rand()%2+2;    }    else    {        mi = -1;        for( i = 0; i < 4; i++ )        {            if( icvFindCircle( pts[shuffles[i][0]], pts[shuffles[i][1]],                               pts[shuffles[i][2]], ¢er, &radius ) >= 0 )            {                radius *= 1.03f;                if( radius < 2.f )                    radius = 2.f;                if( icvIsPtInCircle( pts[shuffles[i][3]], center, radius ) >= 0 &&                    min_radius > radius )                {                    min_radius = radius;                    min_center = center;                    mi = i;                }            }        }        assert( mi >= 0 );        if( mi < 0 )            mi = 0;        k = 3;        center = min_center;        radius = min_radius;        for( i = 0; i < 4; i++ )            idxs[i] = shuffles[mi][i];    }  function_exit:    *_center = center;    *_radius = radius;    /* reorder output points */    for( i = 0; i < 4; i++ )        res_pts[i] = pts[idxs[i]];    for( i = 0; i < 4; i++ )    {        pts[i] = res_pts[i];        assert( icvIsPtInCircle( pts[i], center, radius ) >= 0 );    }    return k;}

怎么判断点是否在圆内或者圆上？

答：计算该点与圆心之间的距离，记为d。 圆的半径记为r。 则计算 r^2 - d^2 。这个值如果>0, 则在圆内，如果=0，则在圆上，如果<0，则在圆外。

CV_INLINE double icvIsPtInCircle( CvPoint2D32f pt, CvPoint2D32f center, float radius ){    double dx = pt.x - center.x;    double dy = pt.y - center.y;    return (double)radius*radius - dx*dx - dy*dy;}

三点确定一个圆，怎么用代码表示呢？

答：我没理清楚公式，先放在这里。

static CvStatusicvFindCircle( CvPoint2D32f pt0, CvPoint2D32f pt1,               CvPoint2D32f pt2, CvPoint2D32f * center, float *radius ){    double x1 = (pt0.x + pt1.x) * 0.5;    double dy1 = pt0.x - pt1.x;    double x2 = (pt1.x + pt2.x) * 0.5;    double dy2 = pt1.x - pt2.x;    double y1 = (pt0.y + pt1.y) * 0.5;    double dx1 = pt1.y - pt0.y;    double y2 = (pt1.y + pt2.y) * 0.5;    double dx2 = pt2.y - pt1.y;    double t = 0;    CvStatus result = CV_OK;    if( icvIntersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 )    {        center->x = (float) (x2 + dx2 * t);        center->y = (float) (y2 + dy2 * t);        *radius = (float) icvDistanceL2_32f( *center, pt0 );    }    else    {        center->x = center->y = 0.f;        radius = 0;        result = CV_NOTDEFINED_ERR;    }    return result;}

其中，icvIntersectLines()函数：

inticvIntersectLines( double x1, double dx1, double y1, double dy1,                   double x2, double dx2, double y2, double dy2, double *t2 ){    double d = dx1 * dy2 - dx2 * dy1;    int result = -1;    if( d != 0 )    {        *t2 = ((x2 - x1) * dy1 - (y2 - y1) * dx1) / d;        result = 0;    }    return result;}