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

一、基本原理

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

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

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

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

 二、代码分析

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

 

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 )
{
//参数：
//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的边缘点的个数，初始化为0
	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 );
	
	//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表示sin
			sx = 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;
				
				//在法线方向和反方向上，距离边缘点的距离为输入的最大半径和最小半径范围内找点
				//每找到一个点，该点的累加器计数就加1
				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;
	
	//累加器中，计数大于阀值的点，被认为可能的圆心点。因为计算各点计数过程中，距离有误差，所以
	//在与阀值比较时，该点计数先要与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_buf
	cvCvtSeqToArray( 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;
		}
	}
}