算法的整体思路：
1.根据设定的阈值canny_threshold，使用canny边缘检测得到可能为圆边缘的点edges。
2.分别计算x,y方向的sobel梯度，用来判断edges点的边缘梯度方向，一个点需计算正、反两个梯度方向(由源码中的k1控制)。
3.遍历所有的edges点，根据设定的圆半径范围[min_radius,max_radius]，在累加投票图中统计圆心出现的次数(参数dp控制累加投票图尺寸)。一个edges点需在所有可能的梯度方向、半径上投票，即需2*(max_radius-min_radius)次投票。最终得到累加投票图adata，同时保存edges点到nz。
4.从第3步的adata中筛选出可能的圆心，判断方法为：圆心投票次数大于acc_threshold，且大于上下左右四个点的投票次数。得到圆心索引centers。
5.按投票次数，从高到低对圆心索引centers排序。此时，adata[centers[0]]即为最高的投票次数，centers[0]对应为圆心位置的索引。
6.遍历centers圆心坐标，计算每个圆心的半径，同时对圆心进行第二次筛选。半径计算方法：找出圆上点最多的半径。圆心第二次筛选方法：圆上点的个数需大于acc_threshold。
源码及注释：
static void
icvHoughCirclesGradient( CvMat* img, float dp, float min_dist,
int min_radius, int max_radius,
int canny_threshold, int acc_threshold,
CvSeq* circles, int circles_max )
{
const int SHIFT = 10, ONE = 1 << SHIFT;//控制圆心坐标计算精度
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;
edges = cvCreateMat( img->rows, img->cols, CV_8UC1 );
cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 );//canny得到可能的圆边缘点
dx = cvCreateMat( img->rows, img->cols, CV_16SC1 );
dy = cvCreateMat( img->rows, img->cols, CV_16SC1 );
cvSobel( img, dx, 1, 0, 3 );//sobel结果用于计算点所在边缘梯度方向，结合半径来推算圆心的位置
cvSobel( img, dy, 0, 1, 3 );
if( dp < 1.f )//控制累加投票图与原图尺寸比例
dp = 1.f;
idp = 1.f/dp;
accum = 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 );//用于存储canny得到的圆边缘点
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]);
// Accumulate circle evidence for each edge pixel
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;
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 = cvRound((vx*idp)*ONE/mag);//梯度方向，乘idp使x1映射投票图坐标
sy = cvRound((vy*idp)*ONE/mag);
x0 = cvRound((x*idp)*ONE);//坐标映射后的圆边缘点
y0 = cvRound((y*idp)*ONE);
// Step from min_radius to max_radius in both directions of the gradient
for(int k1 = 0; k1 < 2; k1++ )
{
x1 = x0 + min_radius * sx;//圆心坐标，投票图坐标系
y1 = y0 + min_radius * sy;
for( 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;
//Find possible circle 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) );
cvCvtSeqToArray( centers, &sort_buf[0] );
icvHoughSortDescent32s( &sort_buf[0], center_count, adata );//按投票次数进行排序
cvClearSeq( 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;
// For each found possible center
// Estimate radius and check support
for( i = 0; i < centers->total; i++ )
{
int ofs = *(int*)cvGetSeqElem( centers, i );//遍历圆心索引，并根据索引值计算圆心坐标x,y
y = ofs/(acols+2);
x = ofs - (y)*(acols+2);
//Calculate circle's center in pixels
float cx = (float)((x + 0.5f)*dp), cy = (float)(( y + 0.5f )*dp);
float start_dist, dist_sum;
float r_best = 0;
int max_count = 0;
// Check distance with previously detected circles
for( j = 0; j < circles->total; j++ )//与已经确定输出的圆比较，距离小于min_dist则排除掉
{
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;
// Estimate best radius
cvStartReadSeq( nz, &reader );
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 );
icvHoughSortDescent32s( &sort_buf[0], nz_count1, (int*)ddata );//按距离值，由大到小排序
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]];//r_cur=j到start_idx的平均距离值
if( (start_idx - j)*r_best >= max_count*r_cur ||//start_idx=当前圆边缘点个数，max_cout=历史最大圆边缘点个数，除以各自的半径，减弱半径越大圆边缘点越多的影响
(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;
}
// Check if the circle has enough support
if( max_count > acc_threshold )
{
float c[3];
c[0] = cx;
c[1] = cy;
c[2] = (float)r_best;
cvSeqPush( circles, c );
if( circles->total > circles_max )
return;
}
}
}
扫一扫
378
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
