圖形可以用一些參數進行表示,標准霍夫變換的原理就是把圖像空間轉換成參數空間(即霍夫空間),例如霍夫變換的直線檢測就是在距離-角度空間內進行檢測。圓可以表示成:
(x-a)2+(y-b)2=r2 (1)
其中a和b表示圓心坐標,r表示圓半徑,因此霍夫變換的圓檢測就是在這三個參數組成的三維空間內進行檢測。
原則上,霍夫變換可以檢測任何形狀。但復雜的形狀需要的參數就多,霍夫空間的維數就多,因此在程序實現上所需的內存空間以及運行效率上都不利於把標准霍夫變換應用於實際復雜圖形的檢測中。所以一些改進的霍夫變換就相繼提出,它們的基本原理就是盡可能減小霍夫空間的維數。
HoughCircles函數實現了圓形檢測,它使用的算法也是改進的霍夫變換——2-1霍夫變換(21HT)。也就是把霍夫變換分為兩個階段,從而減小了霍夫空間的維數。第一階段用於檢測圓心,第二階段從圓心推導出圓半徑。檢測圓心的原理是圓心是它所在圓周所有法線的交匯處,因此只要找到這個交點,即可確定圓心,該方法所用的霍夫空間與圖像空間的性質相同,因此它僅僅是二維空間。檢測圓半徑的方法是從圓心到圓周上的任意一點的距離(即半徑)是相同,只要確定一個閾值,只要相同距離的數量大於該閾值,我們就認為該距離就是該圓心所對應的圓半徑,該方法只需要計算半徑直方圖,不使用霍夫空間。圓心和圓半徑都得到了,那么通過公式1一個圓形就得到了。從上面的分析可以看出,2-1霍夫變換把標准霍夫變換的三維霍夫空間縮小為二維霍夫空間,因此無論在內存的使用上還是在運行效率上,2-1霍夫變換都遠遠優於標准霍夫變換。但該算法有一個不足之處就是由於圓半徑的檢測完全取決於圓心的檢測,因此如果圓心檢測出現偏差,那么圓半徑的檢測肯定也是錯誤的。2-1霍夫變換的具體步驟為:
第一階段:檢測圓心
1.1、對輸入圖像邊緣檢測;
1.2、計算圖形的梯度,並確定圓周線,其中圓周的梯度就是它的法線;
1.3、在二維霍夫空間內,繪出所有圖形的梯度直線,某坐標點上累加和的值越大,說明在該點上直線相交的次數越多,也就是越有可能是圓心;
1.4、在霍夫空間的4鄰域內進行非最大值抑制;
1.5、設定一個閾值,霍夫空間內累加和大於該閾值的點就對應於圓心。
第二階段:檢測圓半徑
2.1、計算某一個圓心到所有圓周線的距離,這些距離中就有該圓心所對應的圓的半徑的值,這些半徑值當然是相等的,並且這些圓半徑的數量要遠遠大於其他距離值相等的數量;
2.2、設定兩個閾值,定義為最大半徑和最小半徑,保留距離在這兩個半徑之間的值,這意味着我們檢測的圓不能太大,也不能太小;
2.3、對保留下來的距離進行排序;
2.4、找到距離相同的那些值,並計算相同值的數量;
2.5、設定一個閾值,只有相同值的數量大於該閾值,才認為該值是該圓心對應的圓半徑;
2.6、對每一個圓心,完成上面的2.1~2.5步驟,得到所有的圓半徑。
HoughCircles函數的原型為:
void HoughCircles(InputArray image,OutputArray circles, int method, double dp, double minDist, double param1=100, double param2=100, int minRadius=0,int maxRadius=0 )
image為輸入圖像,要求是灰度圖像
circles為輸出圓向量,每個向量包括三個浮點型的元素——圓心橫坐標,圓心縱坐標和圓半徑
method為使用霍夫變換圓檢測的算法,Opencv2.4.9只實現了2-1霍夫變換,它的參數是CV_HOUGH_GRADIENT
dp為第一階段所使用的霍夫空間的分辨率,dp=1時表示霍夫空間與輸入圖像空間的大小一致,dp=2時霍夫空間是輸入圖像空間的一半,以此類推
minDist為圓心之間的最小距離,如果檢測到的兩個圓心之間距離小於該值,則認為它們是同一個圓心
param1為邊緣檢測時使用Canny算子的高閾值
param2為步驟1.5和步驟2.5中所共有的閾值
minRadius和maxRadius為所檢測到的圓半徑的最小值和最大值
HoughCircles函數在sources/modules/imgproc/src/hough.cpp文件內被定義:
void cv::HoughCircles( InputArray _image, OutputArray _circles,
int method, double dp, double min_dist,
double param1, double param2,
int minRadius, int maxRadius )
{
//定義一段內存
Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
Mat image = _image.getMat(); //提取輸入圖像矩陣
CvMat c_image = image; //矩陣轉換
//調用cvHoughCircles函數
CvSeq* seq = cvHoughCircles( &c_image, storage, method,
dp, min_dist, param1, param2, minRadius, maxRadius );
//把序列轉換為矩陣
seqToMat(seq, _circles);
}
cvHoughCircles函數為:
CV_IMPL CvSeq*
cvHoughCircles( CvArr* src_image, void* circle_storage,
int method, double dp, double min_dist,
double param1, double param2,
int min_radius, int max_radius )
{
CvSeq* result = 0;
CvMat stub, *img = (CvMat*)src_image;
CvMat* mat = 0;
CvSeq* circles = 0;
CvSeq circles_header;
CvSeqBlock circles_block;
int circles_max = INT_MAX; //輸出最多圓形的數量,設為無窮多
//canny邊緣檢測中雙閾值中的高閾值
int canny_threshold = cvRound(param1);
//累加器閾值
int acc_threshold = cvRound(param2);
img = cvGetMat( img, &stub );
//確保輸入圖像是灰度圖像
if( !CV_IS_MASK_ARR(img))
CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
//內存空間是否存在
if( !circle_storage )
CV_Error( CV_StsNullPtr, "NULL destination" );
//確保參數的正確性
if( dp <= 0 || min_dist <= 0 || canny_threshold <= 0 || acc_threshold <= 0 )
CV_Error( CV_StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" );
//圓的最小半徑要大於0
min_radius = MAX( min_radius, 0 );
//圓的最大半徑如果小於0,則設最大半徑為圖像寬和長度的最大值,
//如果最大半徑小於最小半徑,則設最大半徑為最小半徑加兩個像素的寬度
if( max_radius <= 0 )
max_radius = MAX( img->rows, img->cols );
else if( max_radius <= min_radius )
max_radius = min_radius + 2;
if( CV_IS_STORAGE( circle_storage ))
{
circles = cvCreateSeq( CV_32FC3, sizeof(CvSeq),
sizeof(float)*3, (CvMemStorage*)circle_storage );
}
else if( CV_IS_MAT( circle_storage ))
{
mat = (CvMat*)circle_storage;
if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) ||
CV_MAT_TYPE(mat->type) != CV_32FC3 )
CV_Error( CV_StsBadArg,
"The destination matrix should be continuous and have a single row or a single column" );
circles = cvMakeSeqHeaderForArray( CV_32FC3, sizeof(CvSeq), sizeof(float)*3,
mat->data.ptr, mat->rows + mat->cols - 1, &circles_header, &circles_block );
circles_max = circles->total;
cvClearSeq( circles );
}
else
CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );
//選擇哪種算法檢測圓,目前只有2-1霍夫變換
switch( method )
{
case CV_HOUGH_GRADIENT:
//調用icvHoughCirclesGradient函數
icvHoughCirclesGradient( img, (float)dp, (float)min_dist,
min_radius, max_radius, canny_threshold,
acc_threshold, circles, circles_max );
break;
default:
CV_Error( CV_StsBadArg, "Unrecognized method id" );
}
if( mat )
{
if( mat->cols > mat->rows )
mat->cols = circles->total;
else
mat->rows = circles->total;
}
else
result = circles;
//輸出圓
return result;
}
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 )
{
//為了提高運算精度,定義一個數值的位移量
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;
//nz表示圓周序列,centers表示圓心序列
CvSeq *nz, *centers;
float idp, dr;
CvSeqReader reader;
//創建一個邊緣圖像矩陣
edges = cvCreateMat( img->rows, img->cols, CV_8UC1 );
//第一階段
//步驟1.1,用canny邊緣檢測算法得到輸入圖像的邊緣圖像
cvCanny( img, edges, MAX(canny_threshold/2,1), canny_threshold, 3 );
//創建輸入圖像的水平梯度圖像和垂直梯度圖像
dx = cvCreateMat( img->rows, img->cols, CV_16SC1 );
dy = cvCreateMat( img->rows, img->cols, CV_16SC1 );
//步驟1.2,用Sobel算子法計算水平梯度和垂直梯度
cvSobel( img, dx, 1, 0, 3 );
cvSobel( img, dy, 0, 1, 3 );
/確保累加器矩陣的分辨率不小於1
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 );
//初始化累加器為0
cvZero(accum);
//創建兩個序列,
storage = cvCreateMemStorage();
nz = cvCreateSeq( CV_32SC2, sizeof(CvSeq), sizeof(CvPoint), storage );
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
//步驟1.3,對邊緣圖像計算累加和
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];
//如果當前的像素不是邊緣點,或者水平梯度值和垂直梯度值都為0,則繼續循環。因為如果滿足上面條件,該點一定不是圓周上的點
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);
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;
//在梯度的方向上位移
// r <= max_radius保證了所檢測的圓的半徑一定是小於最大半徑
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;
//在累加器的相應位置上加1
adata[y2*astep + x2]++;
}
//把位移量設置為反方向
sx = -sx; sy = -sy;
}
//把輸入圖像中的當前點(即圓周上的點)的坐標壓入序列圓周序列nz中
pt.x = x; pt.y = y;
cvSeqPush( nz, &pt );
}
}
//計算圓周點的總數
nz_count = nz->total;
//如果總數為0,說明沒有檢測到圓,則退出該函數
if( !nz_count )
return;
//Find possible circle centers
//步驟1.4和1.5,遍歷整個累加器矩陣,找到可能的圓心
for( y = 1; y < arows - 1; y++ )
{
for( x = 1; x < acols - 1; x++ )
{
int base = y*(acols+2) + x;
//如果當前的值大於閾值,並在4鄰域內它是最大值,則該點被認為是圓心
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] )
//把當前點的地址壓入圓心序列centers中
cvSeqPush(centers, &base);
}
}
//計算圓心的總數
center_count = centers->total;
//如果總數為0,說明沒有檢測到圓,則退出該函數
if( !center_count )
return;
//定義排序向量的大小
sort_buf.resize( MAX(center_count,nz_count) );
//把圓心序列放入排序向量中
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 );
//把排好序的圓心重新放入圓心序列中
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 );
//得到圓心在累加器中的坐標位置
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++ )
{
//從序列中提取出圓心
float* c = (float*)cvGetSeqElem( circles, j );
//計算當前圓心與提取出的圓心之間的距離,如果兩者距離小於所設的閾值,則認為兩個圓心是同一個圓心,退出循環
if( (c[0] - cx)*(c[0] - cx) + (c[1] - cy)*(c[1] - cy) < min_dist )
break;
}
//如果j < circles->total,說明當前的圓心已被認為與之前確定作為輸出的圓心是同一個圓心,則拋棄該圓心,返回上面的for循環
if( j < circles->total )
continue;
// Estimate best radius
//第二階段
//開始讀取圓周序列nz
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;
//步驟2.1,計算圓周上的點與當前圓心的距離,即半徑
_r2 = _dx*_dx + _dy*_dy;
//步驟2.2,如果半徑在所設置的最大半徑和最小半徑之間
if(min_radius2 <= _r2 && _r2 <= max_radius2 )
{
//把半徑存入dist_buf內
ddata[k] = _r2;
sort_buf[k] = k;
k++;
}
}
//k表示一共有多少個圓周上的點
int nz_count1 = k, start_idx = nz_count1 - 1;
//nz_count1等於0也就是k等於0,說明當前的圓心沒有所對應的圓,意味着當前圓心不是真正的圓心,所以拋棄該圓心,返回上面的for循環
if( nz_count1 == 0 )
continue;
dist_buf->cols = nz_count1; //得到圓周上點的個數
cvPow( dist_buf, dist_buf, 0.5 ); //求平方根,得到真正的圓半徑
//步驟2.3,對圓半徑進行排序
icvHoughSortDescent32s( &sort_buf[0], nz_count1, (int*)ddata );
dist_sum = start_dist = ddata[sort_buf[nz_count1-1]];
//步驟2.4
for( j = nz_count1 - 2; j >= 0; j-- )
{
float d = ddata[sort_buf[j]];
if( d > max_radius )
break;
//d表示當前半徑值,start_dist表示上一次通過下面if語句更新后的半徑值,dr表示半徑距離分辨率,如果這兩個半徑距離之差大於距離分辨率,說明這兩個半徑一定不屬於同一個圓,而兩次滿足if語句條件之間的那些半徑值可以認為是相等的,即是屬於同一個圓
if( d - start_dist > dr )
{
//start_idx表示上一次進入if語句時更新的半徑距離排序的序號
// start_idx – j表示當前得到的相同半徑距離的數量
//(j + start_idx)/2表示j和start_idx中間的數
//取中間的數所對應的半徑值作為當前半徑值r_cur,也就是取那些半徑值相同的值
float r_cur = ddata[sort_buf[(j + start_idx)/2]];
//如果當前得到的半徑相同的數量大於最大值max_count,則進入if語句
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;
}
// Check if the circle has enough support
//步驟2.5,最終確定輸出
//如果相同半徑的數量大於所設閾值
if( max_count > acc_threshold )
{
float c[3];
c[0] = cx; //圓心的橫坐標
c[1] = cy; //圓心的縱坐標
c[2] = (float)r_best; //所對應的圓的半徑
cvSeqPush( circles, c ); //壓入序列circles內
//如果得到的圓大於閾值,則退出該函數
if( circles->total > circles_max )
return;
}
}
}
下面是用HoughCircles函數進行霍夫變換圓檢測的實例。由於HoughCircles函數內是調用Canny函數進行邊緣檢測,opencv的Canny函數是不包括平滑濾波這一步的,因此為了增強抗干擾能力,在使用HoughCircles函數之前,我們先對原圖進行濾波處理,我們使用的是高斯模糊方法。
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;
using namespace std;
int main( int argc, char** argv )
{
Mat src, gray;
src=imread("coins.jpg");
if( !src.data )
return -1;
cvtColor( src, gray, CV_BGR2GRAY );
//高斯模糊平滑
GaussianBlur( gray, gray, Size(9, 9), 2, 2 );
vector<Vec3f> circles;
//霍夫變換
HoughCircles( gray, circles, CV_HOUGH_GRADIENT, 1, gray.rows/20, 100, 60, 0, 0 );
//在原圖中畫出圓心和圓
for( size_t i = 0; i < circles.size(); i++ )
{
//提取出圓心坐標
Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
//提取出圓半徑
int radius = cvRound(circles[i][2]);
//圓心
circle( src, center, 3, Scalar(0,255,0), -1, 8, 0 );
//圓
circle( src, center, radius, Scalar(0,0,255), 3, 8, 0 );
}
namedWindow( "Circle", CV_WINDOW_AUTOSIZE );
imshow( "Circle", src );
waitKey(0);
return 0;
}
下圖為圓檢測的結果。
從實際運行的結果來看,我們發現HoughCircles函數不足之處是所需要的參數較多,而且每個參數的改變對結果影響都很大,即漏檢和錯檢的幾率很大。
