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 int
cvMinEnclosingCircle( 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 int
icvFindEnslosingCicle4pts_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 CvStatus
icvFindCircle( 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()函数：

int
icvIntersectLines( 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;
}