在http://blog.csdn.net/piaoxuezhong/article/details/58587907中對霍夫變換實現直線檢測進行了匯總，這篇對霍夫變換實現圓形檢測進行匯總~
總體來講，檢測圓形和檢測直線的實現原理相似，在笛卡爾坐標下，圓的表示方程為：(x-a)2+(y-b)2=r2；但在極坐標下，假設已知圓心（x0,y0）,那么圓上的點可以表示為：
所以對于任意一個圓, 假設中心像素點p(x0, y0)像素點已知, 圓半徑已知,則旋轉360度，由極坐標方程可以得到每個點上的坐標。同樣,如果只是知道圖像上像素點, 圓半徑,旋轉360°，則會有一個集中的交點，即圓心，也就是說圓點處的坐標值最強，這正是霍夫變換檢測圓的數學原理。
HoughCircles函數實現了圓形檢測，它使用的算法是改進的霍夫變換，該算法把霍夫變換分為兩個階段，從而減小了霍夫空間的維數。第一階段用于檢測圓心，第二階段從圓心推導出圓半徑。
檢測圓心的方法是圓心是它所在圓周所有法線的交匯處，因此只要找到這個交點，即可確定圓心，該方法所用的霍夫空間與圖像空間的性質相同，因此它僅僅是二維空間。檢測圓半徑的方法是從圓心到圓周上的任意一點的距離相同，首先確定一個閾值，只要計算得到相同距離的數量大于該閾值，就認為該距離就是該圓心所對應的圓半徑，并且該方法只需要計算半徑直方圖，不使用霍夫空間。圓心和圓半徑都得到后，就能確定圓形了。
霍夫變換具體步驟：第一階段：檢測圓心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步驟，得到所有的圓半徑。
opencv函數源碼：
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; CvSeq* seq = cvHoughCircles( &c_image, storage, method, dp, min_dist, param1, param2, minRadius, maxRadius ); seqToMat(seq, _circles);}image為輸入圖像，要求是灰度圖像；circles為輸出圓向量，每個向量包括三個浮點型的元素：圓心橫坐標，圓心縱坐標和圓半徑；method為使用霍夫變換圓檢測的算法，Opencv2.4.13的參數是CV_HOUGH_GRADIENT；dp為第一階段所使用的霍夫空間的分辨率，dp=1時表示霍夫空間與輸入圖像空間的大小一致，dp=2時霍夫空間是輸入圖像空間的一半，以此類推；minDist為圓心之間的最小距離，如果檢測到的兩個圓心之間距離小于該值，則認為它們是同一個圓心；param1為邊緣檢測時使用Canny算子的高閾值；param2為步驟1.5和步驟2.5中所共有的閾值；minRadius和maxRadius為所檢測到的圓半徑的最小值和最大值，默認值0；在opencv檢測圓形調用的函數為：cvHoughCircles跟icvHoughCirclesGradient.具體實現可以參見：點擊打開鏈接，這里只附上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; } } } 使用此函數可以很容易地檢測出圓的圓心，但是它可能找不到合適的圓半徑。通過第八個參數minRadius和第九個參數maxRadius指定最小和最大的圓半徑，來輔助圓檢測的效果?；蛘呖梢灾苯雍雎苑祷匕霃?，因為它們都有著默認值0，單單用HoughCircles函數檢測出來的圓心，然后用額外的一些步驟來進一步確定半徑。霍夫梯度法的缺點（淺墨大神的觀點）<1> 在霍夫梯度法中，使用Sobel導數來計算局部梯度，其可以視作等同于一條局部切線，這不是一個數值穩定的做法。在大多數情況下，這樣做會得到正確的結果，但或許會在輸出中產生一些噪聲。<2> 在邊緣圖像中的整個非0像素集被看做每個中心的候選部分。因此，如果把累加器的閾值設置偏低，算法將要消耗比較長的時間。第三，因為每一個中心只選擇一個圓，如果有同心圓，就只能選擇其中的一個。<3> 因為中心是按照其關聯的累加器值的升序排列的，并且如果新的中心過于接近之前已經接受的中心的話，就不會被保留下來。且當有許多同心圓或者是近似的同心圓時，霍夫梯度法的傾向是保留最大的一個圓?？梢哉f這是一種比較極端的做法，因為在這里默認Sobel導數會產生噪聲，若是對于無窮分辨率的平滑圖像而言的話，這才是必須的。實例：
#include <opencv2/opencv.hpp>#include <opencv2/imgproc/imgproc.hpp>using namespace cv;int main(int argc, char** argv){ Mat srcImage = imread("F://IM_VIDEO//yingbi1.jpg"); Mat midImage, dstImage; imshow("【原始圖】", srcImage); cvtColor(srcImage, midImage, CV_BGR2GRAY);//轉化邊緣檢測后的圖為灰度圖 GaussianBlur(midImage, midImage, Size(9, 9), 2, 2); vector<Vec3f> circles; HoughCircles(midImage, circles, CV_HOUGH_GRADIENT, 1, midImage.rows/20, 100, 100, 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(srcImage, center, 3, Scalar(0, 255, 0), -1, 8, 0); //繪制圓輪廓 circle(srcImage, center, radius, Scalar(155, 50, 255), 3, 8, 0); } imshow("【效果圖】", srcImage); waitKey(0); return 0;}運行結果：
參考：
http://blog.csdn.net/zhaocj/article/details/50454847
http://blog.csdn.net/poem_qianmo/article/details/26977557
新聞熱點
疑難解答
