關於統計概率的霍夫線變換原始碼在下篇補上（我還沒來得及去看），這次直接按照流程把霍夫圓變換擼一遍。

主要參考部落格：

原理分析: （是一篇文章，帶有原始碼分析的）

1.經典霍夫圓變換的原理

霍夫圓變換和霍夫線變換的原理類似。霍夫線變換是兩個引數(r,θ)，霍夫圓需要三個引數，圓心的x,y座標和圓的半徑.如下對應的三個引數c1,c2,c3。

例如:

其形狀和

類似，該函式是由z=x沿z軸旋轉而成的圓錐曲面。(原諒我不會畫圖，不夠深動形象)

對於xy平面的一個點x0,y0(上述對應的點為（1,1）），則對應的由c1,c2,c3組成三維空間的空間曲面。對於c1,c2,c3平面的一個點，則對應的在xy平面它是一個圓。

   對於在x,y平面上的三個點(x0,y0)，(x1,y1)，(x2,y2)，在c1,c2,c3三維空間是對應的三個空間曲面(此時c1,c2,c3相當於常量）

求解這三個方程，我們可以得到c1,c2,c3的值。這說明(x0,y0)，(x1,y1)，(x2,y2)這三個點在由c1,c2,c3所確定的圓上（即c1,c2,c3分別表示圓的圓心x座標、圓心y座標以及圓的半徑），且三個點對應由c1,c2,c3確定的空間的三個空間曲面。進一步說明，在xy平面，三個點在同一個圓上，則它們對應的空間曲面相交於一點（即點(c1,c2,c3))。
故我們如果知道一個邊界上的點的數目，足夠多，且這些點與之對應的空間曲面相交於一點。則這些點構成的邊界，就接近一個圓形。
上述描述的是標準霍夫圓變換的原理，由於三維空間的計算量大大增大的原因, 標準霍夫圓變化很難被應用到實際中。

OpenCV實現的是一個比標準霍夫圓變換更為靈活的檢測方法: 霍夫梯度法, 也叫2-1霍夫變換(21HT),
它的原理依據是圓心一定是在圓上的每個點的模向量上, 這些圓上點模向量的交點就是圓心, 霍夫梯度法的第一步就是找到這些圓心, （圓心包含了圓心處的x和y座標）這樣三維的累加平面就又轉化為二維累加平面. 第二步根據所有候選中心的邊緣非0畫素對其的支援程度來確定半徑。（備註：降維的思想很重要）

2、霍夫梯度法的原理（兩個版本的解釋，對應著看，一定會懂）

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畫素的充分支援，並且到前期被選擇的中心有足夠的距離。則將圓心和半徑壓入到序列中，得以保留。

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

3.Samples:

先給出OpenCV的霍夫圓檢測函式以及引數詳解。

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，表示圓半徑的最大值。

#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<iostream>
using namespace cv;
using namespace std;
Mat g_scrImage,g_scrImagecopy, g_dstImage, g_midImage;
vector<Vec3f>circles;
int g_nthreshold = 15;
static void on_HoughCircles(int, void*);
void showText();

int main(){
	g_scrImage = imread("circle.jpg");
	g_scrImagecopy = imread("circle.jpg");
	imshow("SourceImage", g_scrImage);
	

	namedWindow("Design Sketch");
	createTrackbar("yuzhi", "Design Sketch", &g_nthreshold, 30, on_HoughCircles);
	on_HoughCircles(1.5, 0);
	while (char(waitKey()) != 'q'){}
	return 0;
}
static void on_HoughCircles(int, void*){
	cvtColor(g_scrImage, g_midImage, CV_BGR2GRAY);
	GaussianBlur(g_midImage, g_midImage, Size(9,9), 2, 2);//必須要用高斯blur。
	double g_nthresholdcopy = (double)g_nthreshold / 10;
	HoughCircles( g_midImage, circles, CV_HOUGH_GRADIENT, g_nthresholdcopy, g_midImage.rows/20, 100, 60, 0, 0 );//輸入影象為灰度圖
	for (size_t i = 0; i < circles.size(); i++){
		//提取出圓心座標
		Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
		//提取出圓半徑
		int radius = cvRound(circles[i][2]);
		//圓心
		circle( g_scrImage, center, 3, Scalar(0,0,0), -1, 8, 0 );
		//圓
		circle( g_scrImage, center, radius, Scalar(0,0,255), 3, 8, 0 );
	}
	imshow("Design Sketch", g_scrImage);
	g_scrImage = g_scrImagecopy.clone();
}
void showText(){
	cout << "\n\n\n\t請調整滾動條觀察影象效果~\n\n"
		"\n\n\t\t\t\t\t\t\t\t by淺墨";

}

備註：

1.我加了一個引數來控制霍夫空間解析度，只是為了看一下引數對檢測效果的影響，後面還有幾個引數，而且每個引數的改變對結果影響都很大，即漏檢和錯檢的機率很大。（感興趣可以試一下）

2.畫圓函式詳解

cvCircle(CvArr* img, CvPoint center, int radius, CvScalar color, int thickness=1, int lineType=8, int shift=0)
img為源影象指標
center為畫圓的圓心座標
radius為圓的半徑
color為設定圓的顏色，規則根據B（藍）G（綠）R（紅）
thickness 如果是正數，表示組成圓的線條的粗細程度。否則，表示圓是否被填充
line_type 線條的型別。預設是8
shift 圓心座標點和半徑值的小數點位數

淺墨總結的不錯：截下來，如圖所示

效果如下：

原圖：

dp = 1.5

dp = 1.9

4.原始碼分析

HoughCircles函式在sources/modules/imgproc/src/hough.cpp檔案內被定義：

void cv::HoughCircles( InputArray _image, OutputArray _circles,
                       int method, double dp, double min_dist,
                       double param1, double param2,
                       int minRadius, int maxRadius )
{
    //定義一段記憶體
    Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
    Mat image = _image.getMat();    //提取輸入影象矩陣
    CvMat c_image = image;    //矩陣轉換
    //呼叫cvHoughCircles函式
    CvSeq* seq = cvHoughCircles( &c_image, storage, method,
                    dp, min_dist, param1, param2, minRadius, maxRadius );
    //把序列轉換為矩陣
    seqToMat(seq, _circles);
}

cvHoughCircles函式為：

CV_IMPL CvSeq*
cvHoughCircles( CvArr* src_image, void* circle_storage,
                int method, double dp, double min_dist,
                double param1, double param2,
                int min_radius, int max_radius )
{
    CvSeq* result = 0;

    CvMat stub, *img = (CvMat*)src_image;
    CvMat* mat = 0;
    CvSeq* circles = 0;
    CvSeq circles_header;
    CvSeqBlock circles_block;
    int circles_max = INT_MAX;    //輸出最多圓形的數量，設為無窮多
    //canny邊緣檢測中雙閾值中的高閾值
    int canny_threshold = cvRound(param1);
    //累加器閾值
    int acc_threshold = cvRound(param2);

    img = cvGetMat( img, &stub );
    //確保輸入影象是灰度影象
    if( !CV_IS_MASK_ARR(img))
        CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
    //記憶體空間是否存在
    if( !circle_storage )
        CV_Error( CV_StsNullPtr, "NULL destination" );
    //確保引數的正確性
    if( dp <= 0 || min_dist <= 0 || canny_threshold <= 0 || acc_threshold <= 0 )
        CV_Error( CV_StsOutOfRange, "dp, min_dist, canny_threshold and acc_threshold must be all positive numbers" );
    //圓的最小半徑要大於0
    min_radius = MAX( min_radius, 0 );
    //圓的最大半徑如果小於0，則設最大半徑為影象寬和長度的最大值，
    //如果最大半徑小於最小半徑，則設最大半徑為最小半徑加兩個畫素的寬度
    if( max_radius <= 0 )
        max_radius = MAX( img->rows, img->cols );
    else if( max_radius <= min_radius )
        max_radius = min_radius + 2;

    if( CV_IS_STORAGE( circle_storage ))
    {
        circles = cvCreateSeq( CV_32FC3, sizeof(CvSeq),
            sizeof(float)*3, (CvMemStorage*)circle_storage );
    }
    else if( CV_IS_MAT( circle_storage ))
    {
        mat = (CvMat*)circle_storage;

        if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) ||
            CV_MAT_TYPE(mat->type) != CV_32FC3 )
            CV_Error( CV_StsBadArg,
            "The destination matrix should be continuous and have a single row or a single column" );

        circles = cvMakeSeqHeaderForArray( CV_32FC3, sizeof(CvSeq), sizeof(float)*3,
                mat->data.ptr, mat->rows + mat->cols - 1, &circles_header, &circles_block );
        circles_max = circles->total;
        cvClearSeq( circles );
    }
    else
        CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );
    //選擇哪種演算法檢測圓，目前只有2-1霍夫變換
    switch( method )
    {
    case CV_HOUGH_GRADIENT:
        //呼叫icvHoughCirclesGradient函式
        icvHoughCirclesGradient( img, (float)dp, (float)min_dist,
                                min_radius, max_radius, canny_threshold,
                                acc_threshold, circles, circles_max );
          break;
    default:
        CV_Error( CV_StsBadArg, "Unrecognized method id" );
    }

    if( mat )
    {
        if( mat->cols > mat->rows )
            mat->cols = circles->total;
        else
            mat->rows = circles->total;
    }
    else
        result = circles;
    //輸出圓
    return result;
}

也有兩個版本的註釋，一樣的，結合來看

version1：論文中的

/****************************************************************************************\
*                                     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;
}

兩個結合起來看更加容易理解，本想綜合在一起，但是弄在一起覺得很彆扭。還是分開吧！

version2：（趙老師）

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

霍夫空間是用於統計座標點，解析度沒有必要與原影象一致，所以用dp來調整霍夫空間的大小，從而提高運算速度。而ONE的作用是為了減小運算過程中的累積誤差，提高精度，如：1/3=0.333333，如果保留小數點後5位，則為0.33333，但如果向左位移3位，則為333.33333，這個數的精度要比前一個高，只要我們把所有的數都向左位移3位，最終的結果再向右位移3位，最後的數值精度會提高不少的。基本思想就是這個。

其他應該沒啥了，不過霍夫也算是複雜了，建議先看完線變換，再來看霍夫圓變換，思想以及原始碼都會好理解很多。

折騰了幾天，總算是弄完了，很多也米有解釋很好，歡迎留言討論，any problem is ok!