一、霍夫（圆）变换的广泛使用和简要历史
霍夫变换是一种特征提取方法，被广泛应用在图像处理和计算机视觉应用中。霍夫变换是用来辨别找出物件中的特征，例如：线条。他的算法流程大致如下，给定一个物件、要辨别的形状的种类，算法会在参数空间中执行投票来决定物体的形状，而这是由累加空间（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实现了以下两种霍夫线变换:
- 标准霍夫线变换
- 原理在上面的部分已经说明了. 它能给我们提供一组参数对
的集合来表示检测到的直线 - 在OpenCV 中通过函数 HoughLines 来实现
- 统计概率霍夫线变换
- 这是执行起来效率更高的霍夫线变换. 它输出检测到的直线的端点
- 在OpenCV 中它通过函数 HoughLinesP 来实现
代码
- 这个程序是用来做什么的?
- 加载一幅图片
- 对图片进行 标准霍夫线变换 或是 统计概率霍夫线变换.
- 分别在两个窗口显示原图像和绘出检测到直线的图像.
- 我们将要说明的例程能从 这里 下载。 一个更高级的版本 (能同时演示标准霍夫线变换和统计概率霍夫线变换并带有活动条来改变变换的阈值) 能从 这里 下载。
#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;using namespace std;
void help(){
cout << "
This program demonstrates line finding with the Hough transform.
"
"Usage:
"
"./houghlines <image_name>, Default is pic1.jpg
" << 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 = 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);
}
#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);
}
#endifimshow("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: 储存着检测到的直线的参数对
的容器 * rho : 参数极径 
以像素值为单位的分辨率. 我们使用 1 像素. - 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: 储存着检测到的直线的参数对
的容器 - rho : 参数极径 以像素值为单位的分辨率. 我们使用 1 像素.
- 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.1标准霍夫线变换和统计概率霍夫线变换
OpenCV实现了以下两种霍夫线变换:
- 标准霍夫线变换
- 原理在上面的部分已经说明了. 它能给我们提供一组参数对
- 在OpenCV 中通过函数 HoughLines 来实现
- 统计概率霍夫线变换
- 这是执行起来效率更高的霍夫线变换. 它输出检测到的直线的端点
- 在OpenCV 中它通过函数 HoughLinesP 来实现
代码
- 这个程序是用来做什么的?
- 加载一幅图片
- 对图片进行 标准霍夫线变换 或是 统计概率霍夫线变换.
- 分别在两个窗口显示原图像和绘出检测到直线的图像.
- 我们将要说明的例程能从 这里 下载。 一个更高级的版本 (能同时演示标准霍夫线变换和统计概率霍夫线变换并带有活动条来改变变换的阈值) 能从 这里 下载。
#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;using namespace std;
void help(){
cout << "
This program demonstrates line finding with the Hough transform.
"
"Usage:
"
"./houghlines <image_name>, Default is pic1.jpg
" << 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 = 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);
}
#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);
}
#endifimshow("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: 储存着检测到的直线的参数对 的容器 * rho : 参数极径 以像素值为单位的分辨率. 我们使用 1 像素.
- 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: 储存着检测到的直线的参数对 的容器
- rho : 参数极径 以像素值为单位的分辨率. 我们使用 1 像素.
- theta: 参数极角 以弧度为单位的分辨率. 我们使用 1度 (即CV_PI/180)
- threshold: 要”检测” 一条直线所需最少的的曲线交点 * minLinLength: 能组成一条直线的最少点的数量. 点数量不足的直线将被抛弃.
- maxLineGap: 能被认为在一条直线上的亮点的最大距离.
- dst: 边缘检测的输出图像. 它应该是个灰度图 (但事实上是个二值化图) * lines: 储存着检测到的直线的参数对
通过画出检测到的直线来显示结果.
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函数后储存了霍夫线变换检测到线条的输出矢量。每一条线由具有两个元素的矢量表示，其中，是离坐标原点((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函数后后存储了检测到的线条的输出矢量，每一条线由具有四个元素的矢量表示，其中，和是每个检测到的线段的结束点。
第三个参数，double类型的rho，以像素为单位的距离精度。另一种形容方式是直线搜索时的进步尺寸的单位半径。
第四个参数，double类型的theta，以弧度为单位的角度精度。另一种形容方式是直线搜索时的进步尺寸的单位角度。
第五个参数，int类型的threshold，累加平面的阈值参数，即识别某部分为图中的一条直线时它在累加平面中必须达到的值。大于阈值threshold的线段才可以被检测通过并返回到结果中。
第六个参数，double类型的minLineLength，有默认值0，表示最低线段的长度，比这个设定参数短的线段就不能被显现出来。
第七个参数，double类型的maxLineGap，有默认值0，允许将同一行点与点之间连接起来的最大的距离。
2.3.3OpenCV圆变换函数 HoughCircles
void HoughLines(InputArray image,OutputArray lines, double rho, double theta, int threshold, double srn=0,double stn=0 )C++: void HoughLinesP(InputArray image, OutputArray lines, double rho, double theta, int threshold, double minLineLength=0, double maxLineGap=0 )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)情况时作图（)：会发现曲线的形状类似一条正弦曲线，通常
)：会发现曲线的形状类似一条正弦曲线，通常
这个变化非常有价值，是整个Hough变化的核心！
具体的，通过将霍夫参数空间量化为有限间隔或累加器单元来实现变换。随着算法的运行，每个算法都把转换为一个离散化的 曲线，并且沿着这条曲线的累加器单元被递增。累加器阵列中产生的峰值表示图像中存在相应的直线的有力证据。 换句话说，将每个交点看成一次投票，也就是说，所有点都如此进行计算后，可以设置一个阈值，投票大于这个阈值的可以认为是找到的直线。
那么给定一个点
，通过该点的所有直线的参数
的集合，会在
平面上形成一个三角函数，可由下方程式证明
因此，给定很多点，判断这些点是否共线(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 标准霍夫线变换算法流程- 读取原始图像，并转换成灰度图，利用阈值分割或者边缘检测算子转换成二值化边缘图像
- 初始化霍夫空间， 令所有
- 对于每一个像素点，在参数空间中找出所有满足的对,然后令
- 统计所有的大小，取出的参数（是所设的阈值），从而得到一条直线。
- 将上述流程取出的直线，确定与其相关线段的起始点与终止点（有一些算法，如蝴蝶形状宽度，峰值走廊之类）
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 统计概率霍夫变换算法流程
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);
}
}标准霍夫变换本质上是把图像映射到它的参数空间上，它需要计算所有的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;
}
}
}参考文章：