Hough变换：检测直线和圆
前言：Hough变换是一种在图像中寻找直线和圆的方法。我在实际的项目中，使用到了Hough变换检测圆，效果不错，所以写一篇文章，学习Hough变换的原理，并阅读Hough变换的源码，看看OpenCV是如何实现Hough变换的。Hough变换比较难理解，尤其是圆变换的部分，另外我写的也未必清楚，所以记录下来仅做学习、参考之用。
本文的结构如下：
- 1、Hough变换如何检测出直线
- 2、Hough变换如何检测出圆
- 3、Hough变换 OpenCV源码分析
- 4、实战：Python Hough变换的小例子
1、Hough变换如何检测出直线
我们要解决的问题
在进一步理解Hough变换之前，首先要明确，我们要解决的问题是什么。假设在一个平面上，或者说在一张图片内，存在很多很多的像素点，一些像素点构成了明显的直线，这样的直线可能有多条，我们要做的就是找到这些直线，并将这些直线画出来。比如下图中，房子的屋檐、门窗就包含特征明显的直线，Hough变换算法就是要检测出这些直线。首先来看看Hough变换的数学基础。
Hough变换的数学原理：
Hough变换检测直线的核心是直角坐标系到极坐标系的转换。在平面中，我们往往使用直角坐标系表示一条直线，比如直线：y = kx + b，就是在直角坐标系下给出的。
同样的一条直线，在极坐标系下也能被表示，如下图所示，直线a 也可以被(ρ,θ)表示出来。
当(ρ,θ)是一个确定的值时，直线是确定的，也就是说，当(ρ,θ)是一个确定的值时，比如(10, 30)，它所表示出的直线在平面内是唯一的。同样，当直线确定时，(ρ,θ)也是个确定值。
对于直线a上的点P1和P2，通过P1的直线有多条：【a、b、c、d】；通过P2的直线也有多条：【a、e、f】。我们把这些直线统统用(ρ,θ)来表示，看看能发现什么规律。
我们发现，存在一组相等的(ρ,θ)值：（10，30°）。由该值所确定的直线a，正是P1、P2所共线的直线a。如果两个点共线，那么一定它们会有相同的(ρ,θ)值。反过来讲，如果一个(ρ,θ)值被多个点拥有（或者说出现的次数多），那么这些点一定是共线的，并且这条直线一定是由这个(ρ,θ)值唯一确定的。
理解了上面的内容，就好办了。对于平面内的多个点，我们把通过它们的直线都求出来，然后计算出所有直线的(ρ,θ)，如果发现有一个(ρ,θ)出现的次数很多，那么由这个(ρ,θ)确定的直线，一定被多个点共线。通过统计(ρ,θ)出现的次数，就能找到平面内的直线。
如果在直角坐标系内找直线，我们需要做大量的计算，比如计算斜率，求直线方程等等。但是通过极坐标系，我们仅需要统计(ρ,θ)的数量，就能确定出直线了。因此可以说，Hough变换将大量的计算，转化为了计数问题，这正是它的巧妙之处。
2、Hough变换如何检测出圆
我们要解决的问题
与检测直线类似，在平面内，存在了大量的点。其中，一些点可能共圆，我们要尽可能多的找出这些圆。比如在下面的图像中，我们找到了一个大圆。
Hough圆变换的数学原理：
在平面直角坐标系中，圆的公式是：(x－a)²+(y－b)²=r²，有三个参数a、b、r，圆心坐标为（a，b），半径为r。
假设有一点N (xn, yn)，我们知道通过N的圆有无穷多个，比如下图中，我画出了4个通过N的圆。对于所有通过点N的圆来讲，势必满足公式：(xn－a)²+(yn－b)²=r²成立。其中，因为N是已知的，所以xn，yn是定值。而那些通过N的圆，有无穷多种，圆心、半径都是不一定的，所以a，b，r都是变量。
既然，a、b、r都是变量，我们可以用一个坐标系表示出a、b、r来，姑且叫它参数空间吧。那么，通过N点的所有圆，在这个参数空间中如何表示出来呢？首先来看：
- 当r = 0时，根据公式(xn－a)²+(yn－b)²=r² ，能够计算出 a = xn，b = yn，此时（a、b、r）就是参数空间中的一个点。
- 当r = 1时，我们又能计算出一组(a、b)值，这些(a、b)值均满足公式(xn－a)²+(yn－b)²=r²，在参数空间中刚好构成一个圆。
- 同理，当r等于2、3、4、… ∞时，我们都能求出一组(a、b)值。把r从1到∞得到的所有的（a、b），在参数空间出画出来，就是一个圆锥面。
这只是仅有一个点N的情况，在参数空间中构成了一个圆锥面。
假设现在又有多了一个点M（xm，ym），在参数空间中又构成了另一个圆锥面。如果这两个圆锥面之间有一个交点P(ap, bp，rp)。说明点M、N共圆(ap、bp、rp)。论证如下：
- P(ap,bp,rp)在N构成的圆锥面上，则有(xn－ap)²+(yn－bp)²=rp²;
- P(ap,bp,rp)在M构成的圆锥面上，则有(xm－ap)²+(ym－bp)²=rp²;
- 根据上面两个公式可知，存在一个圆，它的圆心是(ap，bp)，半径是rp，圆上刚好有两个点N、M，满足该圆的方程：(x－ap)²+(y－bp)²=rp²。N（xn，yn）、M（xm，ym）刚好是这个方程的两个解。
当有多个点时，在参数空间中就会有多个圆锥面。多个圆锥面之间会有多个交点，只要我们找到那些交点，就能找到多个点之间的共圆。
对比检测直线和检测圆，我们发现hough变换的实质是坐标系的转换，将计算问题转化为计数问题，这是它最巧妙，也是最有趣的地方。
。。两者对比。。
hough变换检测直线
| hough变换检测圆
|
直角坐标映射到极坐标
| 直角坐标映射到参数空间
|
统计(ρ,θ)出现的次数
| 统计(a，b，r)出现的次数
|
相同的(ρ,θ)说明两个点共线
| 相同的(a，b，r)说明两个点共圆
|
找到出现次数最多的(ρ,θ)，
说明这条直线上的点最多
| 找到出现次数最多的(a，b，r),
说明这个圆上的点最多
|
3、Hough变换 OpenCV源码分析
了解了Hough变换的原理，现在来看看OpenCV是如何实现这么复杂的算法的。
OpenCV2.4.9中，Hough变换的代码实现有1000多行，在这里不能一句一句分析。就把一些比较重要的函数拿出来读读。阅读大佬们的代码，也是一个我自提高的过程。
在OpenCV中，Hough变换的源码在opencv_imgproc/Src/hough.cpp中。在源码中，主要有三个函数：
- void cv::HoughLines( … ); // 标准Hough变换检测直线，调用icvHoughLinesStandard( … )函数
- void cv::HoughLinesP( … ); // Hough变换的优化版，通过分析点的子集，并估计这些点属于一条直线概率。
- void cv::HoughCircles( … ); // Hough变换检测圆
其中，有关HoughLinesP( … )的内容，本文没有涉及，所以在此不作分析。我们只看标准Hough变换检测直线，和检测圆的函数。
1、icvHoughLinesStandard( … )函数的核心部分
// stage 1. fill accumulator，计算(ρ,θ),这里面用(r,θ)来表示。
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] );
r += (numrho - 1) / 2; // numrho是ρ的最大值，或者说最大取值范围
accum[(n+1) * (numrho+2) + r+1]++; // 计数器累加，把(ρ,θ)二维空间用一维数组存储
}
}
}
// stage 2. find local maximums，找极大值
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[total++] = base;
}
}
// stage 3. sort the detected lines by accumulator value
icvHoughSortDescent32s( sort_buf, total, accum ); // 排序
// stage 4. store the first min(total,linesMax) lines to the output buffer，输出直线
linesMax = MIN(linesMax, total); // linesMax是输入参数，表示最多输出几条直线
scale = 1./(numrho+2);
for( i = 0; i < linesMax; i++ )
{
CvLinePolar line; // CvLinePolar 直线的数据结构
int idx = sort_buf[i]; // 找到索引（下标）
int n = cvFloor(idx*scale) - 1; // 向下取整
int r = idx - (n+1)*(numrho+2) - 1;
line.rho = (r - (numrho - 1)*0.5f) * rho;
line.angle = n * theta;
cvSeqPush( lines, &line ); // 用序列存放多条直线
2、icvHoughCirclesGradient( … )函数的核心部分
// Accumulate circle evidence for each edge pixel
// 为每个边缘像素（因为是边缘图，所以边缘像素就是非零点）累加值
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); // dx是对图像在x方向进行sobel求导的结果矩阵
const short* dy_row = (const short*)(dy->data.ptr + y*dy->step); // dy是对图像在y方向进行sobel求导的结果矩阵
for( x = 0; x < cols; x++ )
{
float vx, vy;
int sx, sy, x0, y0, x1, y1, r;
CvPoint pt;
vx = dx_row[x]; // x方向的梯度导数
vy = dy_row[x]; // y方向的梯度导数
// 如果不是边缘像素，或者x方向 y方向的导数为0
if( !edges_row[x] || (vx == 0 && vy == 0) )
continue;
float mag = sqrt(vx*vx+vy*vy); // 使用公式计算：G = sqrt（dx*dx + dy * dy）梯度导数
assert( mag >= 1 ); // 断言
sx = cvRound((vx*idp)*ONE/mag); // x方向的位移量
sy = cvRound((vy*idp)*ONE/mag); // y方向的位移量
x0 = cvRound((x*idp)*ONE); // dp分辨率，idp分辨率的倒数，ONE = 1 << 10
y0 = cvRound((y*idp)*ONE); // 把像素点的位置，转换到累加器的坐标系中
// Step from min_radius to max_radius in both directions of the gradient
// 对x方向和y方向，分别从最小的半径一步一步位移到最大半径，计算累加值
for(int k1 = 0; k1 < 2; k1++ )
{
x1 = x0 + min_radius * sx;
y1 = y0 + min_radius * sy;
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 ) // aclos、arows表示累加器的列、行
break;
adata[y2*astep + x2]++; // 累加器累加
}
sx = -sx; sy = -sy;
}
pt.x = x; pt.y = y;
cvSeqPush( nz, &pt ); // 存到圆周序列中
}
}
nz_count = nz->total;
if( !nz_count ) // 如果没有检测到圆，退出
return;
//Find possible circle centers，找到可能存在的圆心
for( y = 1; y < arows - 1; y++ )
{
for( x = 1; x < acols - 1; x++ )
{
int base = y*(acols+2) + x;// 得到计数值，并以它为基准，看看它是不是局部极大值
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) ); // 重定义sort_buf的大小，sort_buf是一个vector
cvCvtSeqToArray( centers, &sort_buf[0] ); // 把圆心存入sort_buf中
// 从大到小排序 : adata[sort_buf[0]] -> 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); // 累加器中圆心的y
x = ofs - (y)*(acols+2); // 累加器中圆心的x
//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循环执行break了，也就是说两个圆心太近了
if( j < circles->total )
continue;
// Estimate best radius 对当前的圆心，估计它的最佳半径
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中
k++; // k用于累计圆周上点的个数
}
}
int nz_count1 = k, start_idx = nz_count1 - 1;
if( nz_count1 == 0 ) // 相当于k等于0
continue;
dist_buf->cols = nz_count1; // 圆周上点的数量
cvPow( dist_buf, dist_buf, 0.5 ); // 求平方根
icvHoughSortDescent32s( &sort_buf[0], nz_count1, (int*)ddata ); // 对半径排序
dist_sum = start_dist = ddata[sort_buf[nz_count1-1]];
for( j = nz_count1 - 2; j >= 0; j-- )
{
float d = ddata[sort_buf[j]]; // 当前的半径（从小半径开始计算）
if( d > max_radius ) // 大于最大值，跳出
break;
if( d - start_dist > dr ) // 如果两次半径之差大于半径阈值（分辨率），说明不是同一个圆
{
float r_cur = ddata[sort_buf[(j + start_idx)/2]]; // 当前半径
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
if( max_count > acc_threshold )
{
float c[3];
c[0] = cx;
c[1] = cy;
c[2] = (float)r_best;
cvSeqPush( circles, c );
if( circles->total > circles_max )
return;
}
}
4、实战：Python Hough变换的小例子
说了这么多，其实在项目中，我们可以不用自己动手实现Hough变换，OpenCV已经帮我们写得很好了。
实战贴请移步：
参考：
http://blog.163.com/yuyang_tech/blog/static/21605008320130233343990/
