霍夫（圓）變換(hough Transform/hough cirlce Transform)原理和實現

本文轉載自查看原文 2020-06-25 07:39 2250 圖像處理的嚴肅實驗

一、霍夫（圓）變換的廣泛使用和簡要歷史

霍夫變換是一種特征提取方法，被廣泛應用在圖像處理和計算機視覺應用中。霍夫變換是用來辨別找出物件中的特征，例如：線條。他的算法流程大致如下，給定一個物件、要辨別的形狀的種類，算法會在參數空間中執行投票來決定物體的形狀，而這是由累加空間（accumulator space）里的局部最大值來決定。

現在廣泛使用的霍夫變換是由 Richard Duda 和 Peter Hart 在公元1972年發明，並稱之為廣義霍夫變換，廣義霍夫變換和更早前1962年的Paul Hough 的專利有關，這也是其名稱的又來。經典的霍夫變換是偵測圖片中的直線，之后，霍夫變換不僅能識別直線，也能夠識別任何形狀，常見的有圓形、橢圓形。霍夫變換在1962年申請為專利U.S. Patent 3,069,654，其專利名為"辨識復雜圖案的方法及手段"(Method and Means for Recognizing Complex Patterns)。任一條直線可以由斜率和截距來表示，在該專利中，利用斜率和截距來將一條直線參數化，然而這會導致無界的轉換空間(unbounded transform space)，因為斜率有可能是無限大。1981年，因為Dana H. Ballard 的一篇期刊論文 "Generalizing the Hough transform to detect arbitrary shapes"，讓霍夫變換開始流行於電腦視覺界。

“Hough transform”即使在今天進行搜索，其引用次數也遠大於其他經典計算機視覺發現運算符，例如Sobel filter 或 canny edge detection 。

二、調用方法

2.1標准霍夫線變換和統計概率霍夫線變換

OpenCV實現了以下兩種霍夫線變換:

標准霍夫線變換

原理在上面的部分已經說明了. 它能給我們提供一組參數對 $(\theta, r_{\theta})$ 的集合來表示檢測到的直線

在OpenCV 中通過函數 HoughLines 來實現

統計概率霍夫線變換

這是執行起來效率更高的霍夫線變換. 它輸出檢測到的直線的端點 $(x_{0}, y_{0}, x_{1}, y_{1})$

在OpenCV 中它通過函數 HoughLinesP 來實現

代碼

這個程序是用來做什么的?

加載一幅圖片

對圖片進行標准霍夫線變換或是統計概率霍夫線變換.

分別在兩個窗口顯示原圖像和繪出檢測到直線的圖像.

我們將要說明的例程能從這里下載。一個更高級的版本 (能同時演示標准霍夫線變換和統計概率霍夫線變換並帶有活動條來改變變換的閾值) 能從這里下載。

# include "opencv2/highgui/highgui.hpp" # include "opencv2/imgproc/imgproc.hpp" # include < iostream > using namespace cv; using namespace std; void help(){ cout << "\nThis program demonstrates line finding with the Hough transform.\n" "Usage:\n" "./houghlines <image_name>, Default is pic1.jpg\n" << endl;} int main( int argc, char * * argv){ const char * filename = argc > = 2 ? argv[ 1 ] : "pic1.jpg" ; Mat src = imread(filename, 0 ); if (src.empty()) { help(); cout << "can not open " << filename << endl; return - 1 ; } Mat dst, cdst; Canny(src, dst, 50 , 200 , 3 ); cvtColor(dst, cdst, CV_GRAY2BGR); #if 0 vector<Vec2f> lines; HoughLines(dst, lines, 1, CV_PI/180, 100, 0, 0 );   for( size_t i = 0; i < lines.size(); i++ ) { float rho = lines[i][0], theta = lines[i][1]; Point pt1, pt2; double a = cos(theta), b = sin(theta); double x0 = arho, y0 = brho; pt1.x = cvRound(x0 + 1000(-b)); pt1.y = cvRound(y0 + 1000(a)); pt2.x = cvRound(x0 - 1000(-b)); pt2.y = cvRound(y0 - 1000(a)); line( cdst, pt1, pt2, Scalar(0,0,255), 3, CV_AA); } #else vector<Vec4i> lines; HoughLinesP(dst, lines, 1, CV_PI/180, 50, 50, 10 );   for( size_t i = 0; i < lines.size(); i++ ) { Vec4i l = lines[i]; line( cdst, Point(l[0], l[1]), Point(l[2], l[3]), Scalar(0,0,255), 3, CV_AA); } #endif imshow( "source" , src); imshow( "detected lines" , cdst); waitKey(); return 0 ;}

代碼說明

加載圖片

`Mat src = imread(filename, 0 ); if (src.empty()){ help(); cout << "can not open " << filename << endl; return - 1 ;}`

用Canny算子對圖像進行邊緣檢測

`Canny(src, dst, 50 , 200 , 3 );`

現在我們將要執行霍夫線變換. 我們將會說明怎樣使用OpenCV的函數做到這一點:

標准霍夫線變換

首先, 你要執行變換:

`vector < Vec2f > lines; HoughLines(dst, lines, 1 , CV_PI / 180 , 100 , 0 , 0 );`

帶有以下自變量:

dst: 邊緣檢測的輸出圖像. 它應該是個灰度圖 (但事實上是個二值化圖)

lines: 儲存着檢測到的直線的參數對 $(r,\theta)$ 的容器 * rho : 參數極徑 $r$ 以像素值為單位的分辨率. 我們使用 1 像素.

theta: 參數極角 $\theta$ 以弧度為單位的分辨率. 我們使用 1度 (即CV_PI/180)

threshold: 要”檢測” 一條直線所需最少的的曲線交點

srn and stn: 參數默認為0. 查缺OpenCV參考文獻來獲取更多信息.

通過畫出檢測到的直線來顯示結果.

`for ( size_t i = 0 ; i < lines.size(); i ++ ){ float rho = lines[i][ 0 ], theta = lines[i][ 1 ]; Point pt1, pt2; double a = cos(theta), b = sin(theta); double x0 = a * rho, y0 = b * rho; pt1.x = cvRound(x0 + 1000 * ( - b)); pt1.y = cvRound(y0 + 1000 * (a)); pt2.x = cvRound(x0 - 1000 * ( - b)); pt2.y = cvRound(y0 - 1000 * (a)); line( cdst, pt1, pt2, Scalar( 0 , 0 , 255 ), 3 , CV_AA);}`

統計概率霍夫線變換
1.首先, 你要執行變換:

`vector<Vec4i> lines; HoughLinesP(dst, lines, 1, CV_PI/180, 50, 50, 10 );`
帶有以下自變量:

dst: 邊緣檢測的輸出圖像. 它應該是個灰度圖 (但事實上是個二值化圖) * lines: 儲存着檢測到的直線的參數對 $(x_{start}, y_{start}, x_{end}, y_{end})$ 的容器

rho : 參數極徑 $r$ 以像素值為單位的分辨率. 我們使用 1 像素.

theta: 參數極角 $\theta$ 以弧度為單位的分辨率. 我們使用 1度 (即CV_PI/180)

threshold: 要”檢測” 一條直線所需最少的的曲線交點 * minLinLength: 能組成一條直線的最少點的數量. 點數量不足的直線將被拋棄.

maxLineGap: 能被認為在一條直線上的亮點的最大距離.

通過畫出檢測到的直線來顯示結果.

`for ( size_t i = 0 ; i < lines.size(); i ++ ){ Vec4i l = lines[i]; line( cdst, Point(l[ 0 ], l[ 1 ]), Point(l[ 2 ], l[ 3 ]), Scalar( 0 , 0 , 255 ), 3 , CV_AA);}`

顯示原始圖像和檢測到的直線:

`imshow( "source" , src);imshow( "detected lines" , cdst);`

等待用戶按鍵推出程序

`waitKey();`

結果

2.2霍夫圓變換

代碼：

#include <opencv2/imgproc.hpp> #include <opencv2/highgui.hpp> #include <math.h> using namespace cv; using namespace std; int main(int argc, char** argv) { Mat img, gray; if( argc != 2 || !(img=imread(argv[1], 1)).data) return -1; cvtColor(img, gray, COLOR_BGR2GRAY); // smooth it, otherwise a lot of false circles may be detected GaussianBlur( gray, gray, Size(9, 9), 2, 2 ); vector<Vec3f> circles; HoughCircles(gray, circles, HOUGH_GRADIENT,2, gray.rows/4, 200, 100 ); 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]); // draw the circle center circle( img, center, 3, Scalar(0,255,0), -1, 8, 0 ); // draw the circle outline circle( img, center, radius, Scalar(0,0,255), 3, 8, 0 ); } namedWindow( "circles", 1 ); imshow( "circles", img ); waitKey(0);   return 0; }

得益於一個新算法HOUGH_GRADIENT_ALT支持，比以前有很大程度的精度提試對比如下：

2.3相關函數具體調用
2.3.1標准霍夫變換HoughLines()函數

`void HoughLines(InputArray image,OutputArray lines, double rho, double theta, int threshold, double srn=0,double stn=0 )`
第一個參數，InputArray類型的image，輸入圖像，即源圖像，需為8位的單通道二進制圖像，可以將任意的源圖載入進來后由函數修改成此格式后，再填在這里。
第二個參數，InputArray類型的lines，經過調用HoughLines函數后儲存了霍夫線變換檢測到線條的輸出矢量。每一條線由具有兩個元素的矢量 $(p, θ)$ 表示，其中， $p$ 是離坐標原點((0,0)（也就是圖像的左上角）的距離。 $θ$ 是弧度線條旋轉角度（0_{垂直線，π/2}水平線）。
第三個參數，double類型的rho，以像素為單位的距離精度。另一種形容方式是直線搜索時的進步尺寸的單位半徑。PS:Latex中/rho就表示 $ρ$ 。
第四個參數，double類型的theta，以弧度為單位的角度精度。另一種形容方式是直線搜索時的進步尺寸的單位角度。
第五個參數，int類型的threshold，累加平面的閾值參數，即識別某部分為圖中的一條直線時它在累加平面中必須達到的值。大於閾值threshold的線段才可以被檢測通過並返回到結果中。
第六個參數，double類型的srn，有默認值0。對於多尺度的霍夫變換，這是第三個參數進步尺寸rho的除數距離。粗略的累加器進步尺寸直接是第三個參數rho，而精確的累加器進步尺寸為rho/srn。
第七個參數，double類型的stn，有默認值0，對於多尺度霍夫變換，srn表示第四個參數進步尺寸的單位角度theta的除數距離。且如果srn和stn同時為0，就表示使用經典的霍夫變換。否則，這兩個參數應該都為正數。
2.3.2 統計概率霍夫變換(HoughLinesP)

`C++: void HoughLinesP(InputArray image, OutputArray lines, double rho, double theta, int threshold, double minLineLength=0, double maxLineGap=0 )`
第一個參數，InputArray類型的image，輸入圖像，即源圖像，需為8位的單通道二進制圖像，可以將任意的源圖載入進來后由函數修改成此格式后，再填在這里。
第二個參數，InputArray類型的lines，經過調用HoughLinesP函數后后存儲了檢測到的線條的輸出矢量，每一條線由具有四個元素的矢量 $(x_{1}, y_{1}, x_{2}, y_{2} ）$ 表示，其中， $(x_{1}, y_{1})$ 和 $(x_{2}, y_{2})$ 是每個檢測到的線段的結束點。
第三個參數，double類型的rho，以像素為單位的距離精度。另一種形容方式是直線搜索時的進步尺寸的單位半徑。
第四個參數，double類型的theta，以弧度為單位的角度精度。另一種形容方式是直線搜索時的進步尺寸的單位角度。
第五個參數，int類型的threshold，累加平面的閾值參數，即識別某部分為圖中的一條直線時它在累加平面中必須達到的值。大於閾值threshold的線段才可以被檢測通過並返回到結果中。
第六個參數，double類型的minLineLength，有默認值0，表示最低線段的長度，比這個設定參數短的線段就不能被顯現出來。
第七個參數，double類型的maxLineGap，有默認值0，允許將同一行點與點之間連接起來的最大的距離。
2.3.3OpenCV圓變換函數 HoughCircles

`C++: void HoughCircles(InputArray image,OutputArray circles, int method, double dp, double minDist, double param1=100,double param2=100, int minRadius=0, int maxRadius=0)`
第一個參數，InputArray類型的image，輸入圖像，即源圖像，需為8位的灰度單通道圖像。
第二個參數，InputArray類型的circles，經過調用HoughCircles函數后此參數存儲了檢測到的圓的輸出矢量，每個矢量由包含了3個元素的浮點矢量(x, y, radius)表示。
第三個參數，int類型的method，即使用的檢測方法，目前OpenCV中就霍夫梯度法一種可以使用，它的標識符為CV_HOUGH_GRADIENT，在此參數處填這個標識符即可。
第四個參數，double類型的dp，用來檢測圓心的累加器圖像的分辨率於輸入圖像之比的倒數，且此參數允許創建一個比輸入圖像分辨率低的累加器。上述文字不好理解的話，來看例子吧。例如，如果dp= 1時，累加器和輸入圖像具有相同的分辨率。如果dp=2，累加器便有輸入圖像一半那么大的寬度和高度。
第五個參數，double類型的minDist，為霍夫變換檢測到的圓的圓心之間的最小距離，即讓我們的算法能明顯區分的兩個不同圓之間的最小距離。這個參數如果太小的話，多個相鄰的圓可能被錯誤地檢測成了一個重合的圓。反之，這個參數設置太大的話，某些圓就不能被檢測出來了。
第六個參數，double類型的param1，有默認值100。它是第三個參數method設置的檢測方法的對應的參數。對當前唯一的方法霍夫梯度法CV_HOUGH_GRADIENT，它表示傳遞給canny邊緣檢測算子的高閾值，而低閾值為高閾值的一半。
第七個參數，double類型的param2，也有默認值100。它是第三個參數method設置的檢測方法的對應的參數。對當前唯一的方法霍夫梯度法CV_HOUGH_GRADIENT，它表示在檢測階段圓心的累加器閾值。它越小的話，就可以檢測到更多根本不存在的圓，而它越大的話，能通過檢測的圓就更加接近完美的圓形了。
第八個參數，int類型的minRadius,有默認值0，表示圓半徑的最小值。
第九個參數，int類型的maxRadius,也有默認值0，表示圓半徑的最大值。

三、理論

在自動化分析數位圖片的問題里，其中一個常有的子問題是偵測某些簡單的直線、圓形、橢圓形。在多數情況下，邊緣偵測器(edge detector)會先用來做圖片前處理，將原本的圖片變成只含有邊緣的圖片。因為圖片的不完美或是邊緣偵測的不完美，導致有些點(point)或像素(pixel)缺漏，或是有噪聲使得邊緣偵測器所得的邊界偏離了實際的邊界。所以無法直觀的將檢測出的邊緣分成直線、圓形、橢圓形的集合，而霍夫變換解決上述問題，借由霍夫變換算法中的投票步驟，在復雜的參數空間中找到圖形的參數，電腦可以由參數得知該邊緣(edge)是哪種形狀。

3.1 霍夫線變換原理。

3.1.1 在平面直角坐標系(x-y)中，一條直線可以用方程式

表示，而可以視為參數空間中的一點。當直線垂直於軸時，斜率為無限大，若用電腦數值計算時會很不方便。

3.1.2 Duda 和 Hart 提出使用Hesse normal form來表示直線的參數，但是這個本身的原理比較繁瑣，我們可以通過圖像和極坐標的方法來進行簡化。

一方面，x=0的時候y的值，這個點是可以確定的；r/sin(Θ) f(0) = r/sin(Θ)

另一方面，這條直線的斜率（紅色）是tan(PI-Θ) = tan(PI/2 + Θ),根據誘導公式： tan（π/2+α）=－cotα ,也就是斜率為：－cot(Θ)

根據2.1.1中得到的定義，這條直線可以表示為：

進行化簡便可得到（簡單變換）:

這里還有另一種方法，繞過定義直接使用圖像作業，也應該掌握。

3.1.3 對於一個給定點我們在極坐標對極徑極角平面繪出所有通過它的直線, 將得到一條正弦曲線.

例如, 對於給定點(8,6)我們可以繪出下圖 (在平面，

有例如，對於給定點(3,4)情況時作圖（

)：會發現曲線的形狀類似一條正弦曲線，通常

p

這個變化非常有價值，是整個Hough變化的核心！

具體的，通過將霍夫參數空間量化為有限間隔或累加器單元來實現變換。隨着算法的運行，每個算法都把

(x_{i}, y_{i})

轉換為一個離散化的

(r, θ)

曲線，並且沿着這條曲線的累加器單元被遞增。累加器陣列中產生的峰值表示圖像中存在相應的直線的有力證據。換句話說，將每個交點看成一次投票，也就是說

(r, θ)

A (r, θ) = A (r, θ) + 1

，所有點都如此進行計算后，可以設置一個閾值，投票大於這個閾值的可以認為是找到的直線。

那么給定一個點，通過該點的所有直線的參數的集合，會在平面上形成一個三角函數，可由下方程式證明

因此，給定很多點，判斷這些點是否共線(concurrent lines)的問題，經由霍夫變換之后，變成判斷一堆曲線(每一個點在平面上代表一條曲線)是否在平面上相交於同一點的問題(concurrent curves)。繼續上面那張抽象的圖，如果兩個不同點進行上述操作后得到的曲線在平面相交, 這就意味着它們通過同一條直線. 例如,接上面的例子我們繼續繪圖, 得到下圖:

這三條曲線在平面相交於點 (0.925, 9.6), 坐標表示的是參數對

或者是說點

, 點

和點組成的平面內的的直線。

以上的說明表明，一般來說, 一條直線能夠通過在平面尋找交於一點的曲線數量來檢測。而越多曲線交於一點也就意味着這個交點表示的直線由更多的點組成. 一般來說我們可以通過設置直線上點的閾值來定義多少條曲線交於一點我們才認為檢測到了一條直線。

這就是霍夫線變換要做的. 它追蹤圖像中每個點對應曲線間的交點. 如果交於一點的曲線的數量超過了閾值, 那么可以認為這個交點所代表的參數對在原圖像中為一條直線。它的最大特點在於作業域的變化。

3.1.4范例：輸入的圖片中有兩條粗直線，經過霍夫變換后的結果得到accumaltor矩陣，右圖就是把accumaltor矩陣畫出來，越亮值越大，越黑值越小。在右圖中，有兩個很明顯的亮點，這兩個亮點分別代表兩條不同參數的直線，與輸入的圖片(左圖)吻合。然后讀取矩陣的兩個最大值就可以得出這兩條線距畫面中心距離以及角度。

霍夫變換的效率取決於輸入圖片的品質，邊緣必須要正確呈現才能讓霍夫變換有效率，當圖片有噪聲的時候，在霍夫變換的前一級要做抑制噪聲的動作。

3.2、霍夫圓變換原理

霍夫圓變換（CHT）是一種用於數字圖像處理中的基本特征提取技術，用於檢測不完整圖像中的圓。通過在Hough參數空間中進行“投票”，然后在累加器矩陣中選擇局部最大值來生成候選圓。

它是霍夫變換的一種特殊形式。

3.2.1理論

二維圖像中，圓可以這樣來表示：

其中（a，b）是圓的中心，r是半徑。如果二維點（x，y）是固定的，則可以根據（1）找到參數。參數空間將是三維（a，b，r）。並且所有滿足（x，y）的參數都將位於頂點為（x，y，0）的倒置直角錐的表面上。在3D空間中，圓弧參數可以通過2D圓弧上的點定義的許多圓錐曲面的交點來識別。這個過程可以分為兩個階段。第一個階段是固定半徑，然后在2D參數空間中找到最佳的圓心。第二階段是在一維參數空間中找到最佳半徑。

　　圓檢測是視覺處理的基礎應用問題目前主流的檢測算法是基於hough變換的方法其基本思想是將圖像的空間域變換到參數空間用大多數邊界點滿足的某種參數形式來描述圖像中的邊緣曲線通過累加投票求得峰值對應點即為有效圖元信息該方法具有可靠性高對噪聲變形部分區域殘缺邊緣不連續適應性強等特點但傳統hough變換有如下缺陷 1、一到多的參數映射引起的計算量大；2、內存占用空間大；3、參數量化間隔標准確定難。

查找已知半徑R的參數

如果半徑固定，則參數空間將減小為2D（圓心的位置）。對於原始圓上的每個點（x，y），它可以根據（1）定義一個以（x，y）為中心的半徑為R的圓。參數空間中所有此類圓的交點將對應於原始圓的中心點。

考慮原始圖像中一個圓上的4個點（左）。圓圈霍夫變換顯示在右側。注意，假定半徑是已知的。對於原始圖像中四個點（白點）中的每個（x，y），它可以在Hough參數空間中定義一個以半徑r為中心的（x，y）圓。累加器矩陣用於跟蹤相交點。在參數空間中，圓通過的點的投票數將增加一。然后可以找到局部最大值點（右圖中心的紅色點）。最大值的位置（a，b）將是原始圓的中心。

具有已知半徑R的多個圓

用相同的技術可以找到半徑相同的多個圓

累加器矩陣和投票

在實踐中，引入累加器矩陣以找到參數空間中的交點。首先，我們需要使用網格將參數空間划分為“存儲桶”，並根據網格生成累加器矩陣。累加器矩陣中的元素表示參數空間中穿過參數空間中相應網格單元的“圓”數。該號碼也稱為“投票號碼”。最初，矩陣中的每個元素均為零。然后，對於原始空間中的每個“邊緣”點，我們可以在參數空間中制定一個圓，並增加圓通過的網格單元的投票數。此過程稱為“投票”。

投票后，我們可以在累加器矩陣中找到局部最大值。局部最大值的位置與原始空間中的圓心相對應。【如果半徑已經知道的話，將極大程度降低這里的運算量】

查找未知半徑的圓參數

由於參數空間是3D，因此累加器矩陣也將是3D。我們可以遍歷可能的半徑；對於每個半徑，我們使用以前的技術。最后，在3D累加器矩陣中找到局部最大值。累加器數組在3D空間中應為A [x，y，r]。投票應針對每個像素，半徑和theta A [x，y，r] + = 1。

算法：

對於每個A [a，b，r] = 0;

在圖像高斯模糊上處理濾波算法，將圖像轉換為灰度（grayScaling），使Canny運算符，Canny運算符給出圖像的邊緣。

表決累加器中所有可能的圓圈。

累加器A的局部最大投票圈子給出了圈子Hough空間。

累加器的最大投票圈子給出了該圈子。

 
                
                 
                  For each pixel(x,y) 
                 
 
                     For each radius r  
                 = 
                 10 
                  to r  
                 = 
                 60 
                 // the possible radius 
                 
 
                       For each theta t  
                 = 
                 0 
                  to  
                 360 
                    
                 // the possible  theta 0 to 360 
                 
 
                          a  
                 = 
                  x – r  
                 * 
                  cos(t  
                 * 
                  PI  
                 / 
                 180 
                 );  
                 //polar coordinate for center 
                 
 
                          b  
                 = 
                  y – r  
                 * 
                  sin(t  
                 * 
                  PI  
                 / 
                 180 
                 );   
                 //polar coordinate for center 
                 
 
                          A[a,b,r]  
                 += 
                 1 
                 ;  
                 //voting 
                 
 
                       end 
                 
 
                     end 
                 
 
                   end

算法部分，最核心、最精彩的是“域的轉換”部分。這是天才的設想，其思考結果和方法，都非常值得借鑒。

四、代碼解析：

4.0詳細貼心的注釋

告訴你算法的使用方法和近期發展方向

/** @brief Finds circles in a grayscale image using the Hough transform.

The function finds circles in a grayscale image using a modification of the Hough transform.

Example: :

@include snippets/imgproc_HoughLinesCircles.cpp

@note Usually the function detects the centers of circles well. However, it may fail to find correct

radii. You can assist to the function by specifying the radius range ( minRadius and maxRadius ) if

you know it. Or, in the case of #HOUGH_GRADIENT method you may set maxRadius to a negative number

to return centers only without radius search, and find the correct radius using an additional procedure.

It also helps to smooth image a bit unless it's already soft. For example,

GaussianBlur() with 7x7 kernel and 1.5x1.5 sigma or similar blurring may help.

@param image 8-bit, single-channel, grayscale input image.

@param circles Output vector of found circles. Each vector is encoded as 3 or 4 element

floating-point vector \f$(x, y, radius)\f$ or \f$(x, y, radius, votes)\f$ .

@param method Detection method, see #HoughModes. The available methods are #HOUGH_GRADIENT and #HOUGH_GRADIENT_ALT.

@param dp Inverse ratio of the accumulator resolution to the image resolution. For example, if

dp=1 , the accumulator has the same resolution as the input image. If dp=2 , the accumulator has

half as big width and height. For #HOUGH_GRADIENT_ALT the recommended value is dp=1.5,

unless some small very circles need to be detected.

@param minDist Minimum distance between the centers of the detected circles. If the parameter is

too small, multiple neighbor circles may be falsely detected in addition to a true one. If it is

too large, some circles may be missed.

@param param1 First method-specific parameter. In case of #HOUGH_GRADIENT and #HOUGH_GRADIENT_ALT,

it is the higher threshold of the two passed to the Canny edge detector (the lower one is twice smaller).

Note that #HOUGH_GRADIENT_ALT uses #Scharr algorithm to compute image derivatives, so the threshold value

shough normally be higher, such as 300 or normally exposed and contrasty images.

@param param2 Second method-specific parameter. In case of #HOUGH_GRADIENT, it is the

accumulator threshold for the circle centers at the detection stage. The smaller it is, the more

false circles may be detected. Circles, corresponding to the larger accumulator values, will be

returned first. In the case of #HOUGH_GRADIENT_ALT algorithm, this is the circle "perfectness" measure.

The closer it to 1, the better shaped circles algorithm selects. In most cases 0.9 should be fine.

If you want get better detection of small circles, you may decrease it to 0.85, 0.8 or even less.

But then also try to limit the search range [minRadius, maxRadius] to avoid many false circles.

@param minRadius Minimum circle radius.

@param maxRadius Maximum circle radius. If <= 0, uses the maximum image dimension. If < 0, #HOUGH_GRADIENT returns

centers without finding the radius. #HOUGH_GRADIENT_ALT always computes circle radiuses.

@sa fitEllipse, minEnclosingCircle

4.1霍夫線變化

在OpenCV中，霍夫線變換是一種用來尋找直線的方法. 在使用霍夫線變換之前, 首先要對圖像進行邊緣檢測的處理，也即霍夫線變換的直接輸入只能是邊緣二值圖像.

OpenCV中的霍夫線變換有如下三種：

<1>標准霍夫變換（StandardHough Transform，SHT），由HoughLines函數調用。

<2>多尺度霍夫變換（Multi-ScaleHough Transform，MSHT），由HoughLines函數調用。

<3>累計概率霍夫變換（ProgressiveProbabilistic Hough Transform，PPHT），由HoughLinesP函數調用。

4.1 標准霍夫線變換算法流程

讀取原始圖像，並轉換成灰度圖，利用閾值分割或者邊緣檢測算子轉換成二值化邊緣圖像
初始化霍夫空間，令所有 $N u m (θ, p) = 0$
對於每一個像素點 $(x, y)$ ，在參數空間中找出所有滿足 $x c o s θ + y s i n θ = p$ 的 $(θ, p)$ 對,然后令 $N u m (θ, p) = N u m (θ, p) + 1$
統計所有 $N u m (θ, p)$ 的大小，取出 $N u m (θ, p) > τ$ 的參數（ $τ$ 是所設的閾值），從而得到一條直線。
將上述流程取出的直線，確定與其相關線段的起始點與終止點（有一些算法，如蝴蝶形狀寬度，峰值走廊之類）

static void HoughLinesStandard( InputArray src, OutputArray lines, int type, float rho, float theta, int threshold, int linesMax, double min_theta, double max_theta ) { CV_CheckType(type, type == CV_32FC2 || type == CV_32FC3, "Internal error"); Mat img = src.getMat(); int i, j; float irho = 1 / rho; CV_Assert( img.type() == CV_8UC1 ); CV_Assert( linesMax > 0 ); const uchar * image = img.ptr(); //得到圖像的指針 int step = ( int)img.step; //得到圖像的步長(通道) int width = img.cols; //得到圖像的寬 int height = img.rows; //得到圖像的高 int max_rho = width + height; int min_rho = -max_rho; CV_CheckGE(max_theta, min_theta, "max_theta must be greater than min_theta"); //由角度和距離的分辨率得到角度和距離的數量，即霍夫變換后角度和距離的個數 int numangle = cvRound((max_theta - min_theta) / theta); // 霍夫空間,角度方向的大小 int numrho = cvRound(((max_rho - min_rho) + 1) / rho); # if defined HAVE_IPP && IPP_VERSION_X100 > = 810 && !IPP_DISABLE_HOUGH if (type == CV_32FC2 && CV_IPP_CHECK_COND) { IppiSize srcSize = { width, height }; IppPointPolar delta = { rho, theta }; IppPointPolar dstRoi[ 2] = {{(Ipp32f) min_rho, (Ipp32f) min_theta},{(Ipp32f) max_rho, (Ipp32f) max_theta}}; int bufferSize; int nz = countNonZero(img); int ipp_linesMax = std : :min(linesMax, nz numangle /threshold); int linesCount = 0; std : :vector <Vec2f > _lines(ipp_linesMax); IppStatus ok = ippiHoughLineGetSize_8u_C1R(srcSize, delta, ipp_linesMax, &bufferSize); Ipp8u buffer = ippsMalloc_8u_L(bufferSize); if (ok > = 0) {ok = CV_INSTRUMENT_FUN_IPP(ippiHoughLine_Region_8u32f_C1R, image, step, srcSize, (IppPointPolar ) &_lines[ 0], dstRoi, ipp_linesMax, &linesCount, delta, threshold, buffer);}; ippsFree(buffer); if (ok > = 0) { lines.create(linesCount, 1, CV_32FC2); Mat(linesCount, 1, CV_32FC2, &_lines[ 0]).copyTo(lines); CV_IMPL_ADD(CV_IMPL_IPP); return; } setIppErrorStatus(); } # endif //為累加器數組分配內存空間,使用一維數組表示二維空間 Mat _accum = Mat : :zeros( (numangle + 2), (numrho + 2), CV_32SC1 ); //為排序數組分配內存空間 std : :vector < int > _sort_buf; //為正弦和余弦列表分配內存空間 AutoBuffer < float > _tabSin(numangle); AutoBuffer < float > _tabCos(numangle); //分別定義上述內存空間的地址指針 int accum = _accum.ptr < int >(); float tabSin = _tabSin.data(), tabCos = _tabCos.data(); // create sin and cos table createTrigTable( numangle, min_theta, theta, irho, tabSin, tabCos); // stage 1. fill accumulator //執行步驟1，逐點進行霍夫空間變換，並把結果放入累加器數組內 for( i = 0; i < height; i ++ ) for( j = 0; j < width; j ++ ) { //只對圖像的非零值處理，即只對圖像的邊緣像素進行霍夫變換 if( image[i * step + j] != 0 ) for( int n = 0; n < numangle; n ++ ) { //根據公式: ρ = xcosθ + ysinθ //cvRound()函數：四舍五入 int r = cvRound( j * tabCos[n] + i * tabSin[n] ); //numrho是ρ的最大值，或者說最大取值范圍 r += (numrho - 1) / 2; //過界預防 accum[(n + 1) * (numrho + 2) + r + 1] ++; //霍夫空間內的位置 } } // stage 2. find local maximums // 執行步驟2，找到局部極大值，即非極大值抑制 findLocalMaximums( numrho, numangle, threshold, accum, _sort_buf ); // stage 3. sort the detected lines by accumulator value //執行步驟3，對存儲在sort_buf數組內的累加器的數據按由大到小的順序進行排序 std : :sort(_sort_buf.begin(), _sort_buf.end(), hough_cmp_gt(accum)); // stage 4. store the first min(total,linesMax) lines to the output buffer linesMax = std : :min(linesMax, ( int)_sort_buf.size()); double scale = 1. /(numrho + 2); //定義一個尺度 lines.create(linesMax, 1, type); Mat _lines = lines.getMat(); for( i = 0; i < linesMax; i ++ ) // 依據霍夫空間分辨率,計算直線的實際r,theta參數 { //CvLinePolar 直線的數據結構 //CvLinePolar結構在該文件的前面被定義 LinePolar line; //idx為極大值在累加器數組的位置 int idx = _sort_buf[i]; //分離出該極大值在霍夫空間中的位置 //因為n是從0開始的，而之前為了防止越界，所以將所有的n+1了，因此下面要-1，同理r int n = cvFloor(idx scale) - 1; int r = idx - (n + 1) (numrho + 2) - 1; //最終得到極大值所對應的角度和距離 line.rho = (r - (numrho - 1) * 0. 5f) * rho; line.angle = static_cast < float >(min_theta) + n * theta; if (type == CV_32FC2) { _lines.at <Vec2f >(i) = Vec2f(line.rho, line.angle); } else { //存儲到序列內 CV_DbgAssert(type == CV_32FC3); _lines.at <Vec3f >(i) = Vec3f(line.rho, line.angle, ( float)accum[idx]); } } }

 
                        
                         
                         // 霍夫空間,局部最大點,采用四領域判斷,比較。(也可以使8鄰域或者更大的方式),如果不判斷局部最大值,同時選用次大值與最大值,就可能會是兩個相鄰的直線,但實際是一條直線。 
                         
 
                         // 選用最大值,也是去除離散的近似計算帶來的誤差,或合並近似曲線。 
                         
 
                         static  
                         void 
                         
findLocalMaximums(  
                         int numrho,  
                         int numangle,  
                         int threshold, 
                         
                    
                         const  
                         int  
                         *accum, std 
                         : 
                         :vector 
                         < 
                         int 
                         > 
                         & sort_buf ) 
                         
{ 
                         
     
                         for( 
                         int r  
                         =  
                         0; r  
                         < numrho; r 
                         ++ ) 
                         
         
                         for( 
                         int n  
                         =  
                         0; n  
                         < numangle; n 
                         ++ ) 
                         
        { 
                         
             
                         //得到當前值在累加器數組的位置 
                         
             
                         int base  
                         = (n 
                         + 
                         1)  
                         * (numrho 
                         + 
                         2)  
                         + r 
                         + 
                         1; 
                         
             
                         //得到計數值，並以它為基准，看看它是不是局部極大值 
                         
             
                         if( accum[base]  
                         > threshold  
                         && 
                         
                accum[base]  
                         > accum[base  
                         -  
                         1]  
                         && accum[base]  
                         > 
                         = accum[base  
                         +  
                         1]  
                         && 
                         
                accum[base]  
                         > accum[base  
                         - numrho  
                         -  
                         2]  
                         && accum[base]  
                         > 
                         = accum[base  
                         + numrho  
                         +  
                         2] ) 
                         
                 
                         //把極大值位置存入排序數組內——sort_buf 
                         
                sort_buf.push_back(base); 
                         
        } 
                         
}

4.2 統計概率霍夫變換算法流程

標准霍夫變換本質上是把圖像映射到它的參數空間上，它需要計算所有的M個邊緣點，這樣它的運算量和所需內存空間都會很大。如果在輸入圖像中只是處理 $m (m < M)$ 個邊緣點，則這m個邊緣點的選取是具有一定概率性的，因此該方法被稱為概率霍夫變換（Probabilistic Hough Transform）。該方法還有一個重要的特點就是能夠檢測出線端，即能夠檢測出圖像中直線的兩個端點，確切地定位圖像中的直線。
HoughLinesP函數就是利用概率霍夫變換來檢測直線的。它的一般步驟為：

隨機抽取圖像中的一個特征點，即邊緣點，如果該點已經被標定為是某一條直線上的點，則繼續在剩下的邊緣點中隨機抽取一個邊緣點，直到所有邊緣點都抽取完了為止；
對該點進行霍夫變換，並進行累加和計算；
選取在霍夫空間內值最大的點，如果該點大於閾值的，則進行步驟4，否則回到步驟1；
根據霍夫變換得到的最大值，從該點出發，沿着直線的方向位移，從而找到直線的兩個端點；
計算直線的長度，如果大於某個閾值，則被認為是好的直線輸出，回到步驟1。

 
                  
                   
                   static 
                   void 
                   
icvHoughLinesProbabilistic( CvMat 
                   * image, 
                   
                             
                   float rho, 
                   float theta, 
                   int threshold, 
                   
                             
                   int lineLength, 
                   int lineGap, 
                   
                            CvSeq 
                   *lines, 
                   int linesMax ) 
                   
{ 
                   
     
                   //accum為累加器矩陣，mask為掩碼矩陣 
                   
    cv 
                   : 
                   :Mat accum, mask; 
                   
    cv 
                   : 
                   :vector 
                   < 
                   float 
                   > trigtab;     
                   //用於存儲事先計算好的正弦和余弦值 
                   
     
                   //開辟一段內存空間 
                   
    cv 
                   : 
                   :MemStorage storage(cvCreateMemStorage( 
                   0)); 
                   
     
                   //用於存儲特征點坐標，即邊緣像素的位置 
                   
    CvSeq 
                   * seq; 
                   
    CvSeqWriter writer; 
                   
     
                   int width, height;     
                   //圖像的寬和高 
                   
     
                   int numangle, numrho;     
                   //角度和距離的離散數量 
                   
     
                   float ang; 
                   
     
                   int r, n, count; 
                   
    CvPoint pt; 
                   
     
                   float irho 
                   = 
                   1 
                   / rho;     
                   //距離分辨率的倒數 
                   
    CvRNG rng 
                   = cvRNG( 
                   - 
                   1);     
                   //隨機數 
                   
     
                   const 
                   float 
                   * ttab;     
                   //向量trigtab的地址指針 
                   
    uchar 
                   * mdata0;     
                   //矩陣mask的地址指針 
                   
     
                   //確保輸入圖像的正確性 
                   
    CV_Assert( CV_IS_MAT(image) 
                   && CV_MAT_TYPE(image 
                   - 
                   >type) 
                   == CV_8UC1 ); 
                   
 
                   
    width 
                   = image 
                   - 
                   >cols;     
                   //提取出輸入圖像的寬 
                   
    height 
                   = image 
                   - 
                   >rows;     
                   //提取出輸入圖像的高 
                   
     
                   //由角度和距離分辨率，得到角度和距離的離散數量 
                   
    numangle 
                   = cvRound(CV_PI 
                   / theta); 
                   
    numrho 
                   = cvRound(((width 
                   + height) 
                   * 
                   2 
                   + 
                   1) 
                   / rho); 
                   
     
                   //創建累加器矩陣，即霍夫空間 
                   
    accum.create( numangle, numrho, CV_32SC1 ); 
                   
     
                   //創建掩碼矩陣，大小與輸入圖像相同 
                   
    mask.create( height, width, CV_8UC1 ); 
                   
     
                   //定義trigtab的大小，因為要存儲正弦和余弦值，所以長度為角度離散數的2倍 
                   
    trigtab.resize(numangle 
                   * 
                   2); 
                   
     
                   //累加器矩陣清零 
                   
    accum 
                   = cv 
                   : 
                   :Scalar( 
                   0); 
                   
     
                   //避免重復計算，事先計算好所需的所有正弦和余弦值 
                   
     
                   for( ang 
                   = 
                   0, n 
                   = 
                   0; n 
                   < numangle; ang 
                   += theta, n 
                   ++ ) 
                   
    { 
                   
        trigtab[n 
                   * 
                   2] 
                   = ( 
                   float)(cos(ang) 
                   * irho); 
                   
        trigtab[n 
                   * 
                   2 
                   + 
                   1] 
                   = ( 
                   float)(sin(ang) 
                   * irho); 
                   
    } 
                   
     
                   //賦值首地址 
                   
    ttab 
                   = 
                   &trigtab[ 
                   0]; 
                   
    mdata0 
                   = mask.data; 
                   
     
                   //開始寫入序列 
                   
    cvStartWriteSeq( CV_32SC2, 
                   sizeof(CvSeq), 
                   sizeof(CvPoint), storage, 
                   &writer ); 
                   
 
                   
     
                   // stage 1. collect non-zero image points 
                   
     
                   //收集圖像中的所有非零點，因為輸入圖像是邊緣圖像，所以非零點就是邊緣點 
                   
     
                   for( pt.y 
                   = 
                   0, count 
                   = 
                   0; pt.y 
                   < height; pt.y 
                   ++ ) 
                   
    { 
                   
         
                   //提取出輸入圖像和掩碼矩陣的每行地址指針 
                   
         
                   const uchar 
                   * data 
                   = image 
                   - 
                   >data.ptr 
                   + pt.y 
                   *image 
                   - 
                   >step;   
                   // step是每一行的字節大小 此行代碼就轉到當前遍歷的這一行 
                   
        uchar 
                   * mdata 
                   = mdata0 
                   + pt.y 
                   *width; 
                   
         
                   for( pt.x 
                   = 
                   0; pt.x 
                   < width; pt.x 
                   ++ ) 
                   
        { 
                   
             
                   if( data[pt.x] )     
                   //是邊緣點 
                   
            { 
                   
                mdata[pt.x] 
                   = (uchar) 
                   1;     
                   //掩碼的相應位置置1 
                   
                CV_WRITE_SEQ_ELEM( pt, writer );    把該坐標位置寫入序列 
                   
            } 
                   
             
                   else     
                   //不是邊緣點 
                   
                mdata[pt.x] 
                   = 
                   0;     
                   //掩碼的相應位置清0 
                   
        } 
                   
    } 
                   
     
                   //終止寫序列，seq為所有邊緣點坐標位置的序列 
                   
    seq 
                   = cvEndWriteSeq( 
                   &writer ); 
                   
    count 
                   = seq 
                   - 
                   >total;     
                   //得到邊緣點的數量 
                   
 
                   
     
                   // stage 2. process all the points in random order 
                   
     
                   //隨機處理所有的邊緣點 
                   
     
                   for( ; count 
                   > 
                   0; count 
                   -- ) 
                   
    { 
                   
         
                   // choose random point out of the remaining ones 
                   
         
                   //步驟1，在剩下的邊緣點中隨機選擇一個點，idx為不大於count的隨機數 
                   
         
                   int idx 
                   = cvRandInt( 
                   &rng) 
                   % count; 
                   
         
                   //max_val為累加器的最大值，max_n為最大值所對應的角度 
                   
         
                   int max_val 
                   = threshold 
                   - 
                   1, max_n 
                   = 
                   0; 
                   
         
                   //由隨機數idx在序列中提取出所對應的坐標點 
                   
        CvPoint 
                   * point 
                   = (CvPoint 
                   *)cvGetSeqElem( seq, idx ); 
                   
         
                   //定義直線的兩個端點 
                   
        CvPoint line_end[ 
                   2] 
                   = {{ 
                   0, 
                   0}, { 
                   0, 
                   0}}; 
                   
         
                   float a, b; 
                   
         
                   //累加器的地址指針，也就是霍夫空間的地址指針 
                   
         
                   int 
                   * adata 
                   = ( 
                   int 
                   *)accum.data; 
                   
         
                   int i, j, k, x0, y0, dx0, dy0, xflag; 
                   
         
                   int good_line; 
                   
         
                   const 
                   int shift 
                   = 
                   16; 
                   // 
                   
         
                   //提取出坐標點的橫、縱坐標 
                   
        i 
                   = point 
                   - 
                   >y; 
                   
        j 
                   = point 
                   - 
                   >x; 
                   
 
                   
         
                   // "remove" it by overriding it with the last element 
                   
         
                   //用序列中的最后一個元素覆蓋掉剛才提取出來的隨機坐標點 
                   
         
                   *point 
                   = 
                   *(CvPoint 
                   *)cvGetSeqElem( seq, count 
                   - 
                   1 ); 
                   
 
                   
         
                   // check if it has been excluded already (i.e. belongs to some other line) 
                   
         
                   //檢測這個坐標點是否已經計算過，也就是它已經屬於其他直線 
                   
         
                   //因為計算過的坐標點會在掩碼矩陣mask的相對應位置清零 
                   
         
                   if( 
                   !mdata0[i 
                   *width 
                   + j] )     
                   //該坐標點被處理過 
                   
             
                   continue;     
                   //不做任何處理，繼續主循環 
                   
 
                   
         
                   // update accumulator, find the most probable line 
                   
         
                   //步驟2，更新累加器矩陣，找到最有可能的直線 
                   
         
                   for( n 
                   = 
                   0; n 
                   < numangle; n 
                   ++, adata 
                   += numrho ) 
                   
        { 
                   
             
                   //由角度計算距離 
                   
            r 
                   = cvRound( j 
                   * ttab[n 
                   * 
                   2] 
                   + i 
                   * ttab[n 
                   * 
                   2 
                   + 
                   1] ); 
                   
            r 
                   += (numrho 
                   - 
                   1) 
                   / 
                   2; 
                   //防止r為負數 
                   
             
                   //在累加器矩陣的相應位置上數值加1，並賦值給val 
                   
             
                   int val 
                   = 
                   ++adata[r]; 
                   
             
                   //更新最大值，並得到它的角度 
                   
             
                   if( max_val 
                   < val ) 
                   
            { 
                   
                max_val 
                   = val; 
                   
                max_n 
                   = n; 
                   
            } 
                   
        } 
                   
 
                   
         
                   // if it is too "weak" candidate, continue with another point 
                   
         
                   //步驟3，如果上面得到的最大值小於閾值，則放棄該點，繼續下一個點的計算 
                   
         
                   if( max_val 
                   < threshold ) 
                   
             
                   continue; 
                   
 
                   
         
                   // from the current point walk in each direction 
                   
         
                   // along the found line and extract the line segment 
                   
         
                   //步驟4，從當前點出發，沿着它所在直線的方向前進，直到達到端點為止 
                   
        a 
                   = 
                   -ttab[max_n 
                   * 
                   2 
                   + 
                   1];     
                   //a=-sinθ 
                   
        b 
                   = ttab[max_n 
                   * 
                   2];     
                   //b=cosθ 
                   
         
                   //當前點的橫、縱坐標值 
                   
        x0 
                   = j; 
                   
        y0 
                   = i; 
                   
         
                   //確定當前點所在直線的角度是在45度～135度之間，還是在0～45或135度～180度之間 
                   
         
                   if( fabs(a) 
                   > fabs(b) )     
                   //在45度～135度之間 
                   
        { 
                   
            xflag 
                   = 
                   1;     
                   //置標識位，標識直線的粗略方向 
                   
             
                   //確定橫、縱坐標的位移量 
                   
            dx0 
                   = a 
                   > 
                   0 
                   ? 
                   1 
                   : 
                   - 
                   1;   
                   
            dy0 
                   = cvRound( b 
                   *( 
                   1 
                   << shift) 
                   /fabs(a) ); 
                   
             
                   //確定縱坐標 
                   
            y0 
                   = (y0 
                   << shift) 
                   + ( 
                   1 
                   << (shift 
                   - 
                   1)); 
                   
        } 
                   
         
                   else     
                   //在0～45或135度～180度之間 
                   
        { 
                   
            xflag 
                   = 
                   0;    
                   //清標識位 
                   
             
                   //確定橫、縱坐標的位移量 
                   
            dy0 
                   = b 
                   > 
                   0 
                   ? 
                   1 
                   : 
                   - 
                   1; 
                   
            dx0 
                   = cvRound( a 
                   *( 
                   1 
                   << shift) 
                   /fabs(b) ); 
                   
             
                   //確定橫坐標 
                   
            x0 
                   = (x0 
                   << shift) 
                   + ( 
                   1 
                   << (shift 
                   - 
                   1)); 
                   
        } 
                   
         
                   //搜索直線的兩個端點 
                   
         
                   for( k 
                   = 
                   0; k 
                   < 
                   2; k 
                   ++ ) 
                   
        { 
                   
             
                   //gap表示兩條直線的間隙，x和y為搜索位置，dx和dy為位移量 
                   
             
                   int gap 
                   = 
                   0, x 
                   = x0, y 
                   = y0, dx 
                   = dx0, dy 
                   = dy0; 
                   
             
                   //搜索第二個端點的時候，反方向位移 
                   
             
                   if( k 
                   > 
                   0 ) 
                   
                dx 
                   = 
                   -dx, dy 
                   = 
                   -dy; 
                   
 
                   
             
                   // walk along the line using fixed-point arithmetics, 
                   
             
                   // stop at the image border or in case of too big gap 
                   
             
                   //沿着直線的方向位移，直到到達圖像的邊界或大的間隙為止 
                   
             
                   for( ;; x 
                   += dx, y 
                   += dy ) 
                   
            { 
                   
                uchar 
                   * mdata; 
                   
                 
                   int i1, j1; 
                   
                 
                   //確定新的位移后的坐標位置 
                   
                 
                   if( xflag ) 
                   
                { 
                   
                    j1 
                   = x; 
                   
                    i1 
                   = y 
                   >> shift; 
                   
                } 
                   
                 
                   else 
                   
                { 
                   
                    j1 
                   = x 
                   >> shift; 
                   
                    i1 
                   = y; 
                   
                } 
                   
                 
                   //如果到達了圖像的邊界，停止位移，退出循環 
                   
                 
                   if( j1 
                   < 
                   0 
                   || j1 
                   > 
                   = width 
                   || i1 
                   < 
                   0 
                   || i1 
                   > 
                   = height ) 
                   
                     
                   break; 
                   
                 
                   //定位位移后掩碼矩陣位置 
                   
                mdata 
                   = mdata0 
                   + i1 
                   *width 
                   + j1; 
                   
 
                   
                 
                   // for each non-zero point: 
                   
                 
                   //    update line end, 
                   
                 
                   //    clear the mask element 
                   
                 
                   //    reset the gap 
                   
                 
                   //該掩碼不為0，說明該點可能是在直線上 
                   
                 
                   if( 
                   *mdata ) 
                   
                { 
                   
                    gap 
                   = 
                   0;     
                   //設置間隙為0 
                   
                     
                   //更新直線的端點位置 
                   
                    line_end[k].y 
                   = i1; 
                   
                    line_end[k].x 
                   = j1; 
                   
                } 
                   
                     
                   //掩碼為0，說明不是直線，但仍繼續位移，直到間隙大於所設置的閾值為止 
                   
                 
                   else 
                   if( 
                   ++gap 
                   > lineGap )     
                   //間隙加1 
                   
                     
                   break; 
                   
            } 
                   
        } 
                   
         
                   //步驟5，由檢測到的直線的兩個端點粗略計算直線的長度 
                   
         
                   //當直線長度大於所設置的閾值時，good_line為1，否則為0 
                   
        good_line 
                   = abs(line_end[ 
                   1].x 
                   - line_end[ 
                   0].x) 
                   > 
                   = lineLength 
                   || 
                   
                    abs(line_end[ 
                   1].y 
                   - line_end[ 
                   0].y) 
                   > 
                   = lineLength; 
                   
         
                   //再次搜索端點，目的是更新累加器矩陣和更新掩碼矩陣，以備下一次循環使用 
                   
         
                   for( k 
                   = 
                   0; k 
                   < 
                   2; k 
                   ++ ) 
                   
        { 
                   
             
                   int x 
                   = x0, y 
                   = y0, dx 
                   = dx0, dy 
                   = dy0; 
                   
 
                   
             
                   if( k 
                   > 
                   0 ) 
                   
                dx 
                   = 
                   -dx, dy 
                   = 
                   -dy; 
                   
 
                   
             
                   // walk along the line using fixed-point arithmetics, 
                   
             
                   // stop at the image border or in case of too big gap 
                   
             
                   for( ;; x 
                   += dx, y 
                   += dy ) 
                   
            { 
                   
                uchar 
                   * mdata; 
                   
                 
                   int i1, j1; 
                   
 
                   
                 
                   if( xflag ) 
                   
                { 
                   
                    j1 
                   = x; 
                   
                    i1 
                   = y 
                   >> shift; 
                   
                } 
                   
                 
                   else 
                   
                { 
                   
                    j1 
                   = x 
                   >> shift; 
                   
                    i1 
                   = y; 
                   
                } 
                   
 
                   
                mdata 
                   = mdata0 
                   + i1 
                   *width 
                   + j1; 
                   
 
                   
                 
                   // for each non-zero point: 
                   
                 
                   //    update line end, 
                   
                 
                   //    clear the mask element 
                   
                 
                   //    reset the gap 
                   
                 
                   if( 
                   *mdata ) 
                   
                { 
                   
                     
                   //if語句的作用是清除那些已經判定是好的直線上的點對應的累加器的值，避免再次利用這些累加值 
                   
                     
                   if( good_line )     
                   //在第一次搜索中已經確定是好的直線 
                   
                    { 
                   
                         
                   //得到累加器矩陣地址指針 
                   
                        adata 
                   = ( 
                   int 
                   *)accum.data; 
                   
                         
                   for( n 
                   = 
                   0; n 
                   < numangle; n 
                   ++, adata 
                   += numrho ) 
                   
                        { 
                   
                            r 
                   = cvRound( j1 
                   * ttab[n 
                   * 
                   2] 
                   + i1 
                   * ttab[n 
                   * 
                   2 
                   + 
                   1] ); 
                   
                            r 
                   += (numrho 
                   - 
                   1) 
                   / 
                   2; 
                   
                            adata[r] 
                   --;     
                   //相應的累加器減1 
                   
                        } 
                   
                    } 
                   
                     
                   //搜索過的位置，不管是好的直線，還是壞的直線，掩碼相應位置都清0，這樣下次就不會再重復搜索這些位置了，從而達到減小計算邊緣點的目的 
                   
                     
                   *mdata 
                   = 
                   0; 
                   
                } 
                   
                 
                   //如果已經到達了直線的端點，則退出循環 
                   
                 
                   if( i1 
                   == line_end[k].y 
                   && j1 
                   == line_end[k].x ) 
                   
                     
                   break; 
                   
            } 
                   
        } 
                   
         
                   //如果是好的直線 
                   
         
                   if( good_line ) 
                   
        { 
                   
            CvRect lr 
                   = { line_end[ 
                   0].x, line_end[ 
                   0].y, line_end[ 
                   1].x, line_end[ 
                   1].y }; 
                   
             
                   //把兩個端點壓入序列中 
                   
            cvSeqPush( lines, 
                   &lr ); 
                   
             
                   //如果檢測到的直線數量大於閾值，則退出該函數 
                   
             
                   if( lines 
                   - 
                   >total 
                   > 
                   = linesMax ) 
                   
                 
                   return; 
                   
        } 
                   
    } 
                   
}

4.3霍夫圓變化

HoughCircles函數實現了圓形檢測，它使用的算法也是改進的霍夫變換——2-1霍夫變換（21HT）。也就是把霍夫變換分為兩個階段，從而減小了霍夫空間的維數。第一階段用於檢測圓心，第二階段從圓心推導出圓半徑。檢測圓心的原理是圓心是它所在圓周所有法線的交匯處，因此只要找到這個交點，即可確定圓心，該方法所用的霍夫空間與圖像空間的性質相同，因此它僅僅是二維空間。檢測圓半徑的方法是從圓心到圓周上的任意一點的距離（即半徑）是相同，只要確定一個閾值，只要相同距離的數量大於該閾值，我們就認為該距離就是該圓心所對應的圓半徑，該方法只需要計算半徑直方圖，不使用霍夫空間。圓心和圓半徑都得到了，那么通過公式1一個圓形就得到了。從上面的分析可以看出，2-1霍夫變換把標准霍夫變換的三維霍夫空間縮小為二維霍夫空間，因此無論在內存的使用上還是在運行效率上，2-1霍夫變換都遠遠優於標准霍夫變換。但該算法有一個不足之處就是由於圓半徑的檢測完全取決於圓心的檢測，因此如果圓心檢測出現偏差，那么圓半徑的檢測肯定也是錯誤的。

version 1:

2-1霍夫變換的具體步驟為：

1）首先對圖像進行邊緣檢測，調用opencv自帶的cvCanny()函數，將圖像二值化，得到邊緣圖像。
2）對邊緣圖像上的每一個非零點。采用cvSobel()函數，計算x方向導數和y方向的導數，從而得到梯度。從邊緣點，沿着梯度和梯度的反方向，對參數指定的min_radius到max_radium的每一個像素，在累加器中被累加。同時記下邊緣圖像中每一個非0點的位置。
3）從（二維）累加器中這些點中選擇候選中心。這些中心都大於給定的閾值和其相鄰的四個鄰域點的累加值。
4）對於這些候選中心按照累加值降序排序，以便於最支持的像素的中心首次出現。
5）對於每一個中心，考慮到所有的非0像素（非0，梯度不為0），這些像素按照與其中心的距離排序，從最大支持的中心的最小距離算起，選擇非零像素最支持的一條半徑。
6）如果一個中心受到邊緣圖像非0像素的充分支持，並且到前期被選擇的中心有足夠的距離。則將圓心和半徑壓入到序列中，得以保留。

 
                
                 
                 /****************************************************************************************\
*                                     Circle Detection                                   *
\****************************************************************************************/ 
                 
 
                 /*------------------------------------霍夫梯度法------------------------------------------*/ 
                 
 
                 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, R_THRESH 
                 = 
                 30;   
                 //One=1024,1左移10位2*10，R_THRESH是起始值，賦給max_count,后續會被覆蓋。 
                 
    cv 
                 : 
                 :Ptr 
                 <CvMat 
                 > dx, dy;   
                 //Ptr是智能指針模板，將CvMat對象封裝成指針 
                 
    cv 
                 : 
                 :Ptr 
                 <CvMat 
                 > edges, accum, dist_buf; 
                 //edges邊緣二值圖像,accum為累加器圖像,dist_buf存放候選圓心到滿足條件的邊緣點的半徑 
                 
    std 
                 : 
                 :vector 
                 < 
                 int 
                 > sort_buf; 
                 //用來進行排序的中間對象。在adata累加器排序中，其存放的是offset即偏移位置,int型。在ddata距離排序中，其存儲的和下標是一樣的值。 
                 
    cv 
                 : 
                 :Ptr 
                 <CvMemStorage 
                 > storage; 
                 //內存存儲器。創建的序列用來向其申請內存空間。 
                 
  
                 
    
                 int x, y, i, j, k, center_count, nz_count; 
                 //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; 
                 //rows,cols邊緣圖像的行數和列數,arows,acols是累加器圖像的行數和列數 
                 
    
                 int astep, 
                 *adata; 
                 //adata指向累加器數據域的首地址，用位置作為下標,astep為累加器每行的大小，以字節為單位 
                 
    
                 float 
                 * ddata; 
                 //ddata即dist_data，距離數據 
                 
    CvSeq 
                 *nz, 
                 *centers; 
                 //nz為非0，即邊界，centers為存放的候選中心的位置。 
                 
    
                 float idp, dr; 
                 //idp即inv_dp,dp的倒數 
                 
    CvSeqReader reader; 
                 //順序讀取序列中的每個值 
                 
  
                 
    edges 
                 = cvCreateMat( img 
                 - 
                 >rows, img 
                 - 
                 >cols, CV_8UC1 ); 
                 //邊緣圖像 
                 
    cvCanny( img, edges, MAX(canny_threshold 
                 / 
                 2, 
                 1), canny_threshold, 
                 3 ); 
                 //調用canny，變為二值圖像，0和非0即0和255 
                 
  
                 
    dx 
                 = cvCreateMat( img 
                 - 
                 >rows, img 
                 - 
                 >cols, CV_16SC1 ); 
                 //16位單通道圖像，用來存儲二值邊緣圖像的x方向的一階導數 
                 
    dy 
                 = cvCreateMat( img 
                 - 
                 >rows, img 
                 - 
                 >cols, CV_16SC1 ); 
                 //y方向的 
                 
    cvSobel( img, dx, 
                 1, 
                 0, 
                 3 ); 
                 //計算x方向的一階導數 
                 
    cvSobel( img, dy, 
                 0, 
                 1, 
                 3 ); 
                 //計算y方向的一階導數 
                 
  
                 
    
                 if( dp 
                 < 
                 1.f ) 
                 //控制dp不能比1小 
                 
        dp 
                 = 
                 1.f; 
                 
    idp 
                 = 
                 1.f 
                 /dp; 
                 
    accum 
                 = cvCreateMat( cvCeil(img 
                 - 
                 >rows 
                 *idp) 
                 + 
                 2, cvCeil(img 
                 - 
                 >cols 
                 *idp) 
                 + 
                 2, CV_32SC1 ); 
                 //cvCeil返回不小於參數的最小整數。32為單通道 
                 
    cvZero(accum); 
                 //初始化累加器為0 
                 
  
                 
    storage 
                 = cvCreateMemStorage(); 
                 //創建內存存儲器，使用默認參數0.默認大小為64KB 
                 
    nz 
                 = cvCreateSeq( CV_32SC2, 
                 sizeof(CvSeq), 
                 sizeof(CvPoint), storage ); 
                 //創建序列，用來存放非0點 
                 
    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; 
                 //cvMat對象的union對象的i成員成員 
                 
    
                 //step是矩陣中行的長度，單位為字節。我們使用到的矩陣是accum它的深度是CV_32SC1即32位的int 型。 
                 
    
                 //如果我們知道一個指針如int* p;指向數組中的一個元素， 則可以通過p+accum->step/adata[0]來使指針移動到p指針所指元素的，正對的下一行元素 
                 
    astep 
                 = accum 
                 - 
                 >step 
                 / 
                 sizeof(adata[ 
                 0]); 
                 
  
                 
    
                 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); 
                 //存儲 x方向sobel一階導數的矩陣的每一行的指向第一個元素的指針 
                 
        
                 const 
                 short 
                 * dy_row 
                 = ( 
                 const 
                 short 
                 *)(dy 
                 - 
                 >data.ptr 
                 + y 
                 *dy 
                 - 
                 >step); 
                 //y 
                 
  
                 
        
                 //遍歷邊緣的二值圖像和偏導數的圖像 
                 
        
                 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) ) 
                 //如果在邊緣圖像（存儲邊緣的二值圖像）某一點如A(x0,y0)==0則對一下點進行操作。vx和vy同時為0，則下一個 
                 
                
                 continue; 
                 
  
                 
            
                 float mag 
                 = sqrt(vx 
                 *vx 
                 +vy 
                 *vy); 
                 //求梯度圖像 
                 
            assert( mag 
                 > 
                 = 
                 1 ); 
                 //如果mag為0,說明沒有邊緣點，則stop。這里用了assert宏定義 
                 
            sx 
                 = cvRound((vx 
                 *idp) 
                 *ONE 
                 /mag); 
                 //  vx為該點的水平梯度（梯度幅值已經歸一化）；ONE為為了用整數運算代替浮點數引入的一個因子，為2^10 
                 
            sy 
                 = cvRound((vy 
                 *idp) 
                 *ONE 
                 /mag); 
                 
  
                 
            x0 
                 = cvRound((x 
                 *idp) 
                 *ONE); 
                 
            y0 
                 = cvRound((y 
                 *idp) 
                 *ONE); 
                 
  
                 
            
                 for( k 
                 = 
                 0; k 
                 < 
                 2; k 
                 ++ ) 
                 //k=0在梯度方向，k=1在梯度反方向對累加器累加。這里之所以要反向，因為對於一個圓上一個點，從這個點沿着斜率的方向的，最小半徑到最大半徑。在圓的另一邊與其相對應的點，有對應的效果。 
                 
            { 
                 
                x1 
                 = x0 
                 + min_radius 
                 * sx;  
                 
                y1 
                 = y0 
                 + min_radius 
                 * sy; 
                 
  
                 
                
                 for( r 
                 = min_radius; r 
                 < 
                 = max_radius; x1 
                 += sx, y1 
                 += sy, r 
                 ++ ) 
                 //x1=x1+sx即,x1=x0+min_radius*sx+sx=x0+(min_radius+1)*sx求得下一個點。sx為斜率 
                 
                { 
                 
                    
                 int x2 
                 = x1 
                 >> SHIFT, y2 
                 = y1 
                 >> SHIFT; 
                 //變回真實的坐標 
                 
                    
                 if( ( 
                 unsigned)x2 
                 > 
                 = ( 
                 unsigned)acols 
                 ||    
                 //如果x2大於累加器的行 
                 
                        ( 
                 unsigned)y2 
                 > 
                 = ( 
                 unsigned)arows ) 
                 
                        
                 break; 
                 
                    adata[y2 
                 *astep 
                 + x2] 
                 ++; 
                 //由於c語言是按行存儲的。即等價於對accum數組進行了操作。 
                 
                } 
                 
  
                 
                sx 
                 = 
                 -sx; sy 
                 = 
                 -sy; 
                 
            } 
                 
  
                 
            pt.x 
                 = x; pt.y 
                 = y; 
                 
            cvSeqPush( nz, 
                 &pt ); 
                 //把非零邊緣並且梯度不為0的點壓入到堆棧 
                 
        } 
                 
    } 
                 
  
                 
    nz_count 
                 = nz 
                 - 
                 >total; 
                 
    
                 if( 
                 !nz_count ) 
                 //如果nz_count==0則返回 
                 
        
                 return; 
                 
  
                 
    
                 for( y 
                 = 
                 1; y 
                 < arows 
                 - 
                 1; y 
                 ++ )      
                 //這里是從1到arows-1，因為如果是圓的話，那么圓的半徑至少為1，即圓心至少在內層里面 
                 
    { 
                 
        
                 for( x 
                 = 
                 1; x 
                 < acols 
                 - 
                 1; x 
                 ++ ) 
                 
        { 
                 
            
                 int base 
                 = y 
                 *(acols 
                 + 
                 2) 
                 + x; 
                 //計算位置，在accum圖像中 
                 
            
                 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] );   
                 //把序列轉換成數組，即把序列centers中的數據放入到sort_buf的容器中。 
                 
  
                 
    icvHoughSortDescent32s( 
                 &sort_buf[ 
                 0], center_count, adata ); 
                 //快速排序，根據sort_buf中的值作為下標，依照adata中對應的值進行排序,將累加值大的下標排到前面 
                 
    cvClearSeq( centers ); 
                 //清空序列 
                 
    cvSeqPushMulti( centers, 
                 &sort_buf[ 
                 0], center_count ); 
                 //重新將中心的下標存入centers 
                 
  
                 
  
                 
    dist_buf 
                 = cvCreateMat( 
                 1, nz_count, CV_32FC1 ); 
                 //創建一個32為浮點型的一個行向量 
                 
    ddata 
                 = dist_buf 
                 - 
                 >data.fl; 
                 //使ddata執行這個行向量的首地址 
                 
  
                 
    dr 
                 = dp; 
                 
    min_dist 
                 = MAX( min_dist, dp ); 
                 //如果輸入的最小距離小於dp，則設在為dp 
                 
    min_dist 
                 *= min_dist; 
                 
  
                 
  
                 
    
                 for( i 
                 = 
                 0; i 
                 < centers 
                 - 
                 >total; i 
                 ++ )    
                 //對於每一個中心點 
                 
    { 
                 
  
                 
  
                 
        
                 int ofs 
                 = 
                 *( 
                 int 
                 *)cvGetSeqElem( centers, i ); 
                 //獲取排序的中心位置，adata值最大的元素，排在首位  ，offset偏移位置 
                 
        y 
                 = ofs 
                 /(acols 
                 + 
                 2) 
                 - 
                 1; 
                 //這里因為edge圖像比accum圖像小兩個邊。 
                 
        x 
                 = ofs 
                 - (y 
                 + 
                 1) 
                 *(acols 
                 + 
                 2) 
                 - 
                 1; 
                 //求得y坐標 
                 
        
                 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 
                 ++ ) 
                 //中存儲已經找到的圓；若當前候選圓心與其中的一個圓心距離<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 ) 
                 //當前候選圓心與任意已檢測的圓心距離不小於min_dist時，才有j==circles->total 
                 
            
                 continue; 
                 
        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; 
                 //sort_buf同上，但是這里的sort_buf的下標值和元素值是相同的,重新利用 
                 
                k 
                 ++; 
                 //k和j是兩個游標 
                 
            } 
                 
        } 
                 
  
                 
        
                 int nz_count1 
                 = k, start_idx 
                 = nz_count1 
                 - 
                 1; 
                 
        
                 if( nz_count1 
                 == 
                 0 ) 
                 
            
                 continue;   
                 //如果一個候選中心到（非零邊界且梯度>0）確定的點的距離中，沒有滿足條件的，則從下一個中心點開始。 
                 
        dist_buf 
                 - 
                 >cols 
                 = nz_count1; 
                 //用來存放真是的滿足條件的非零元素(三個約束：非零點，梯度不為0，到圓心的距離在min_radius和max_radius中間) 
                 
        cvPow( dist_buf, dist_buf, 
                 0. 
                 5 ); 
                 //對dist_buf中的元素開根號.求得半徑 
                 
        icvHoughSortDescent32s( 
                 &sort_buf[ 
                 0], nz_count1, ( 
                 int 
                 *)ddata ); 
                 ////對與圓心的距離按降序排列，索引值在sort_buf中 
                 
  
                 
  
                 
        dist_sum 
                 = start_dist 
                 = ddata[sort_buf[nz_count1 
                 - 
                 1]]; 
                 //dist距離，選取半徑最小的作為起始值 
                 
  
                 
  
                 
        
                 //下邊for循環里面是一個算法。它定義了兩個游標（指針）start_idx和j，j是外層循環的控制變量。而start_idx為控制當兩個相鄰的數組ddata的數據發生變化時，即d-start_dist>dr時，的步進。 
                 
        
                 for( j 
                 = nz_count1 
                 - 
                 2; j 
                 > 
                 = 
                 0; j 
                 -- ) 
                 //從小到大。選出半徑支持點數最多的半徑 
                 
        { 
                 
            
                 float d 
                 = ddata[sort_buf[j]]; 
                 
  
                 
            
                 if( d 
                 > max_radius ) 
                 //如果求得的候選圓點到邊界的距離大於參數max_radius，則停止，因為d是第一個出現的最小的（按照從大到小的順序排列的） 
                 
                
                 break; 
                 
  
                 
            
                 if( d 
                 - start_dist 
                 > dr ) 
                 //如果當前的距離減去最小的>dr(==dp) 
                 
            { 
                 
                
                 float r_cur 
                 = ddata[sort_buf[(j 
                 + start_idx) 
                 / 
                 2]]; 
                 //當前半徑設為符合該半徑的中值，j和start_idx相當於兩個游標 
                 
                
                 if( (start_idx 
                 - j) 
                 *r_best 
                 > 
                 = max_count 
                 *r_cur 
                 ||   
                 //如果數目相等時，它會找半徑較小的那個。這里是判斷支持度的算法 
                 
                    (r_best 
                 < FLT_EPSILON 
                 && start_idx 
                 - j 
                 > 
                 = max_count) ) 
                 //程序這個部分告訴我們，無法找到同心圓，它會被外層最大，支持度最多（支持的點最多）所覆蓋。 
                 
                { 
                 
                    r_best 
                 = r_cur; 
                 //如果 符合當前半徑的點數(start_idx - j)/ 當前半徑>= 符合之前最優半徑的點數/之前的最優半徑 || 還沒有最優半徑時,且點數>30時；其實直接把r_best初始值置為1即可省去第二個條件 
                 
                    max_count 
                 = start_idx 
                 - j; 
                 //maxcount變為符合當前半徑的點數，更新max_count值，后續的支持度大的半徑將會覆蓋前面的值。 
                 
                } 
                 
                start_dist 
                 = d; 
                 
                start_idx 
                 = j; 
                 
                dist_sum 
                 = 
                 0; 
                 //如果距離改變較大，則重置distsum為0，再在下面的式子中置為當前值 
                 
            } 
                 
            dist_sum 
                 += d; 
                 //如果距離改變較小，則加上當前值（dist_sum）在這里好像沒有用處。 
                 
        } 
                 
  
                 
        
                 if( max_count 
                 > R_THRESH ) 
                 //符合條件的圓周點大於閾值30，則將圓心、半徑壓棧 
                 
        { 
                 
            c[ 
                 0] 
                 = cx; 
                 
            c[ 
                 1] 
                 = cy; 
                 
            c[ 
                 2] 
                 = ( 
                 float)r_best; 
                 
            cvSeqPush( circles, c ); 
                 
            
                 if( circles 
                 - 
                 >total 
                 > circles_max ) 
                 //circles_max是個很大的數,其值為INT_MAX 
                 
                
                 return; 
                 
        } 
                 
    } 
                 
} 
                 
  
                 
CV_IMPL CvSeq 
                 * 
                 
cvHoughCircles1( 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; 
                 
    
                 int canny_threshold 
                 = cvRound(param1); 
                 //cvRound返回和參數最接近的整數值，對一個double類型進行四舍五入 
                 
    
                 int acc_threshold 
                 = cvRound(param2); 
                 
  
                 
    img 
                 = cvGetMat( img, 
                 &stub ); 
                 //將img轉化成為CvMat對象 
                 
  
                 
    
                 if( 
                 !CV_IS_MASK_ARR(img))   
                 //圖像必須為8位，單通道圖像 
                 
        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" ); 
                 
  
                 
    min_radius 
                 = MAX( min_radius, 
                 0 ); 
                 
    
                 if( max_radius 
                 < 
                 = 
                 0 ) 
                 //用來控制當使用默認參數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; 
                 
  
                 
        
                 //數組應該是CV_32FC3類型的單列數組。 
                 
        
                 if( 
                 !CV_IS_MAT_CONT( mat 
                 - 
                 >type ) 
                 || (mat 
                 - 
                 >rows 
                 != 
                 1 
                 && mat 
                 - 
                 >cols 
                 != 
                 1) 
                 || 
                 //連續,單列,CV_32FC3類型 
                 
            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 ); 
                 //由於是單列，故elem_size為mat->rows+mat->cols-1 
                 
        circles_max 
                 = circles 
                 - 
                 >total; 
                 
        cvClearSeq( circles ); 
                 //清空序列的內容（如果傳入的有數據的話） 
                 
    } 
                 
    
                 else 
                 
        CV_Error( CV_StsBadArg, 
                 "Destination is not CvMemStorage* nor CvMat*" ); 
                 
  
                 
    
                 switch( method ) 
                 
    { 
                 
    
                 case CV_HOUGH_GRADIENT 
                 : 
                 
        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 ) 
                 //給定一個指向圓存儲的數組指針值，則返回0，即NULL 
                 
    { 
                 
        
                 if( mat 
                 - 
                 >cols 
                 > mat 
                 - 
                 >rows ) 
                 //因為不知道傳入的是列向量還是行向量。 
                 
            mat 
                 - 
                 >cols 
                 = circles 
                 - 
                 >total; 
                 
        
                 else 
                 
            mat 
                 - 
                 >rows 
                 = circles 
                 - 
                 >total; 
                 
    } 
                 
    
                 else 
                 //如果是傳入的是內存存儲器，則返回一個指向一個序列的指針。 
                 
        result 
                 = circles; 
                 
  
                 
    
                 return result; 
                 
}

version 2:

第一階段：檢測圓心

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步驟，得到所有的圓半徑

 
                 
                  
                  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; 
                  
        } 
                  
    } 
                  
}

參考文章：

https://www.cnblogs.com/kk17/p/9693132.html

來自為知筆記(Wiz)

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 Hough transform(霍夫變換) 概率霍夫變換（Progressive Probabilistic Hough Transform）原理詳解 Matlab 霍夫變換 ( Hough Transform）直線檢測 Hough Transform直線檢測第三章霍夫變換(Hough Transform) Hough變換原理（ZZ） hough變換檢測直線和圓 Hough變換原理 opencv —— HoughCircles 霍夫圓變換原理及圓檢測 hough變換是如何檢測出直線和圓的？