opencv代码分析--hough变换识别圆

一、基本原理

opencv中圆识别的基本原理如下：

1、canny算子求图像的单像素二值化边缘

2、假设我们需要找半径为R的所有圆，则对于边缘图中的每一个边缘点，该边缘点的切线的法线方向上(正负两个方向)，寻找到该边缘点距离为R的点，将该点的计数加1(初始化所有点的计数都是0)

3、找到计数值大于门限值的点，即圆心所在的点

 二、代码分析

代码在/modules\imgproc\src\hough.cpp文件icvHoughCirclesGradient函数中

 

static voidicvHoughCirclesGradient( CvMat* img, float dp, float min_dist,                         int min_radius, int max_radius,                         int canny_threshold, int acc_threshold,                         CvSeq* circles, int circles_max ){//参数：//img: 输入图像//dp: 识别精度,1.0表示按照原图精度//min_dist: 圆心点位置识别精度//min_radius: 所需要找的圆的最小半径//max_radius：所需要找的圆的最大半径//canny_threshold：canny算子的高阀值//acc_threshold：累加器阀值，计数大于改阀值的点即被认为是可能的圆心//circles: 保存找到的符合条件的所有圆//circles_max: 最多需要的找到的圆的个数const int SHIFT = 10, ONE = 1 << SHIFT, R_THRESH = 30;cv::Ptr<CvMat> dx, dy;cv::Ptr<CvMat> edges, accum, dist_buf;std::vector<int> sort_buf;cv::Ptr<CvMemStorage> storage;int x, y, i, j, k, center_count, nz_count;float min_radius2 = (float)min_radius*min_radius;float max_radius2 = (float)max_radius*max_radius;int rows, cols, arows, acols;int astep, *adata;float* ddata;CvSeq *nz, *centers;float idp, dr;CvSeqReader reader;//canny算子求单像素二值化边缘，保存在edges变量中edges = cvCreateMat( img->rows, img->cols, CV_8UC1 );cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 );//sobel算子求水平和垂直方向的边缘，用于计算边缘点的法线方向dx = cvCreateMat( img->rows, img->cols, CV_16SC1 );dy = cvCreateMat( img->rows, img->cols, CV_16SC1 );cvSobel( img, dx, 1, 0, 3 );cvSobel( img, dy, 0, 1, 3 );//dp表示识别精度if( dp < 1.f )dp = 1.f;idp = 1.f/dp;//accum用作累加器，包含图像中每一个点的计数。图像中每一个点都有一个计数，点的计数表示每一个canny边缘点法线方向上，//到该点距离为R的边缘点的个数，初始化为0accum = cvCreateMat( cvCeil(img->rows*idp)+2, cvCeil(img->cols*idp)+2, CV_32SC1 );cvZero(accum);storage = cvCreateMemStorage();nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage );//centers用于保存可能的圆心点centers = cvCreateSeq( CV_32SC1, sizeof(CvSeq), sizeof(int), storage );rows = img->rows;cols = img->cols;arows = accum->rows - 2;acols = accum->cols - 2;adata = accum->data.i;astep = accum->step/sizeof(adata[0]);//以下这个循环用于获取所有可能的圆边缘点，存储在nz中，同时设置//累加器中的值for( y = 0; y < rows; y++ ){const uchar* edges_row = edges->data.ptr + y*edges->step;const short* dx_row = (const short*)(dx->data.ptr + y*dx->step);const short* dy_row = (const short*)(dy->data.ptr + y*dy->step);for( x = 0; x < cols; x++ ){float vx, vy;int sx, sy, x0, y0, x1, y1, r, k;CvPoint pt;vx = dx_row[x];vy = dy_row[x];if( !edges_row[x] || (vx == 0 && vy == 0) )continue;float mag = sqrt(vx*vx+vy*vy);assert( mag >= 1 );//sx表示cos, sy表示sinsx = cvRound((vx*idp)*ONE/mag);sy = cvRound((vy*idp)*ONE/mag);x0 = cvRound((x*idp)*ONE);y0 = cvRound((y*idp)*ONE);//循环两次表示需要计算两个方向，法线方向和法线的反方向for( k = 0; k < 2; k++ ){//半径方向的水平增量和垂直增量x1 = x0 + min_radius * sx;y1 = y0 + min_radius * sy;//在法线方向和反方向上，距离边缘点的距离为输入的最大半径和最小半径范围内找点//每找到一个点，该点的累加器计数就加1for( r = min_radius; r <= max_radius; x1 += sx, y1 += sy, r++ ){int x2 = x1 >> SHIFT, y2 = y1 >> SHIFT;if( (unsigned)x2 >= (unsigned)acols || (unsigned)y2 >= (unsigned)arows )break;adata[y2*astep + x2]++;}//反方向sx = -sx; sy = -sy;}//保存可能的圆边缘点pt.x = x; pt.y = y;cvSeqPush( nz, &pt );}}nz_count = nz->total;if( !nz_count )return;//累加器中，计数大于阀值的点，被认为可能的圆心点。因为计算各点计数过程中，距离有误差，所以//在与阀值比较时，该点计数先要与4邻域内的各个点的计数比较，最大者才能和阀值比较。可能的圆心//点保存在centers中for( y = 1; y < arows - 1; y++ ){for( x = 1; x < acols - 1; x++ ){int base = y*(acols+2) + x;if( adata[base] > acc_threshold &&    adata[base] > adata[base-1] && adata[base] > adata[base+1] &&    adata[base] > adata[base-acols-2] && adata[base] > adata[base+acols+2] )cvSeqPush(centers, &base);}}center_count = centers->total;if( !center_count )return;sort_buf.resize( MAX(center_count,nz_count) );//链表结构的certers转化成连续存储结构sort_bufcvCvtSeqToArray( centers, &sort_buf[0] );//经过icvHoughSortDescent32s函数后，以sort_buf中元素作为adata数组下标, //adata中的元素降序排列, 即adata[sort_buf[0]]是adata所有元素中最大的, //adata[sort_buf[center_count-1]]是所有元素中最小的icvHoughSortDescent32s( &sort_buf[0], center_count, adata );cvClearSeq( centers );//经过排序后的元素，重新以链表形式存储到centers中cvSeqPushMulti( centers, &sort_buf[0], center_count );dist_buf = cvCreateMat( 1, nz_count, CV_32FC1 );ddata = dist_buf->data.fl;dr = dp;min_dist = MAX( min_dist, dp );min_dist *= min_dist;//对于每一个可能的圆心点，计算所有边缘点到该圆心点的距离。由于centers中的//元素已经经过前面排序，所以累加器计数最大的可能圆心点最先进行下面的操作for( i = 0; i < centers->total; i++ ){int ofs = *(int*)cvGetSeqElem( centers, i );y = ofs/(acols+2) - 1;x = ofs - (y+1)*(acols+2) - 1;float cx = (float)(x*dp), cy = (float)(y*dp);float start_dist, dist_sum;float r_best = 0, c[3];int max_count = R_THRESH;//如果该可能的圆心点和已经确认的圆心点的距离小于阀值，则表示//这个圆心点和已经确认的圆心点是同一个点for( j = 0; j < circles->total; j++ ){float* c = (float*)cvGetSeqElem( circles, j );if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )break;}if( j < circles->total )continue;cvStartReadSeq( nz, &reader );//求所有边缘点到当前圆心点的距离，符合条件的距离值保存在ddata中for( j = k = 0; j < nz_count; j++ ){CvPoint pt;float _dx, _dy, _r2;CV_READ_SEQ_ELEM( pt, reader );_dx = cx - pt.x; _dy = cy - pt.y;_r2 = _dx*_dx + _dy*_dy;if(min_radius2 <= _r2 && _r2 <= max_radius2 ){ddata[k] = _r2;sort_buf[k] = k;k++;}}int nz_count1 = k, start_idx = nz_count1 - 1;if( nz_count1 == 0 )continue;dist_buf->cols = nz_count1;cvPow( dist_buf, dist_buf, 0.5 );//经过如下处理后，以sort_buf中元素作为ddata数组下标,ddata中的元素降序排列, //即ddata[sort_buf[0]]是ddata所有元素中最大的, ddata[sort_buf[nz_count1-1]]//是所有元素中最小的icvHoughSortDescent32s( &sort_buf[0], nz_count1, (int*)ddata );//对所有的距离值做处理，求出最可能圆半径值，max_count为到圆心的距离为最可能半径值的点的个数dist_sum = start_dist = ddata[sort_buf[nz_count1-1]];for( j = nz_count1 - 2; j >= 0; j-- ){float d = ddata[sort_buf[j]];if( d > max_radius )break;if( d - start_dist > dr ){float r_cur = ddata[sort_buf[(j + start_idx)/2]];if( (start_idx - j)*r_best >= max_count*r_cur ||(r_best < FLT_EPSILON && start_idx - j >= max_count) ) {r_best = r_cur;max_count = start_idx - j;}start_dist = d;start_idx = j;dist_sum = 0;}dist_sum += d;}//max_count大于阀值，表示这几个边缘点构成一个圆if( max_count > R_THRESH ){c[0] = cx;c[1] = cy;c[2] = (float)r_best;cvSeqPush( circles, c );if( circles->total > circles_max )return;}}}