浅谈霍夫变换

霍夫变换是计算机视觉和图像处理领域中一种重要的特征提取算法，它最常见的应用是在噪声中找到诸如直线、圆、椭圆等特征。我第一次接触霍夫变换是在研一的图像处理课上，当时不是很明白，直到后来做项目时需要定制该算法，动手写起来才算是真正的理解该算法。

作为低端码农，我首先关心的是：输入和输出是什么，下面用实例讲一下：

如下图所示，我们的输入是这些灰色的点，很明显的可以看出，这些点中隐藏着一条直线。所以我们输出就是在这些点中找到那条最显著的直线。

这个任务对一般人来说太容易了，然而要让电脑顺利完成检测直线的任务，却是个头疼的问题。我们可以分析一下，图中点只有一部分在直线上，其他的点离直线很远，这些无关的点可以看成噪声数据。分成两种情况：1) 噪声数据很少时，我们可以用最小二乘法(  Least Square )求出直线系数：y = ax + b; 2)噪声数据较多时，可以采用非常暴力的方法：任意取出两个点，求出直线方程，计算其他点到这条线的距离，当距离小于某个阈值时认为点在直线上，并记录下一共有多少个点对应于这条直线。假设有100个点，一共可以计算出50 * 49 = 2450条直线，在这些直线中选出点数最多那条就是所要的结果。这是穷举法，当然耗时耗力，但绝对正确。事实上，由该穷举法衍生出了另外一种重要的参数估计方法：RANSAC，我会在以后的博客中介绍。

霍夫变换主要用于第二种情况，当噪声很多时，我们如何找到直线？我们可以在数据上穷举(每两个点算一次)，最终的目标是算出直线的参数a和b，反其道行之，为什么不穷举所有a,b呢？对于每个a和b，可以算出在直线上的点数，在众多的ab中选出点数最多的就行。

这种想法在实现时碰到一个问题：如何穷举a和b？看起来ab的范围似乎是无穷的。我们换一种直线的表示方法：

如上图所示，红色的直线由两个参数表示：原点到直线的距离r和直线与坐标轴的夹角sita，对任意一条直线来说，夹角sita在[-180, 180]在，距离r的长度也容易计算，所以穷举直线参数成为可能。

当给定(r, sita)时，直线方程就确定了。然而穷举策略中"计算每个点到直线的距离"这一步实在是太浪费资源了。我们得接着想办法。这时候就得Hough大神出现了，大神说，换个角度看世界，两个点确定一条直线，也就是确定一对(r, sita)，那么一个点呢？当然是很多很多对(r, sita)啦！每个点( x, y )相当于公民，为自己支持的直线( r, sita )投一票，总统就确定了！当给定一个点，比如说( 1, 1 )，r和sita的关系是： r = cos(sita) + sin(sita)，就是正弦曲线，当三个点在一条直线上时，三条正弦曲线就会相交于一点：

下图是很多点的情况，颜色越亮表示是两条正弦曲线交点的次数越多。

知道r和sita的范围后，我们可以对这个二维空间进行栅格化，当一个点确定的正弦曲线经过某栅格时，就把这个栅格的票数加1。

好了，现在我们可以写出霍夫变换的算法了：

Step1:：确定r和sita的范围。

Step2：对每个点计算( r, sita ), r = x * cos(sita) + y * sin(sita)。

Step3：统计，找出得票最多的点对( r, sita )。

每一步都有些注意点，Step1中，sita的范围是取( -180, 180 ), ( 0, 180 )还是( -90, 90 )？这要取决于r！当sita在(0, 180)时，r可正可负，sita在(0, 360)时r只能大于0;Step3选择时，要选那些”峰值“点，即该点的得票数要大于周围的点。

我把Hough变换封装成类(根据OpenCV1.0中的HoughTransform() )，如果有遗漏或Bug希望大家指正。

参考文献：

[1] http://en.wikipedia.org/wiki/Hough_Transform

[2] http://docs.opencv.org/doc/tutorials/imgproc/imgtrans/hough_lines/hough_lines.html

#pragma once#ifndef H_HOUGH_H#define H_HOUGH_H#include "Commons.h"class CHough{public:float m_AngResolution;float m_RhoResolution;float m_MinAng;float m_MaxAng;float m_MinRho;float m_MaxRho;int m_PtsNumThr;float m_DistThrSquared;double m_k;double m_b;vector<double> m_vk;vector<double> m_vb;vector<SPointXY> m_vCosSinTab;vector<SPointXY> m_vOrigPoints;CHough(){m_AngResolution  = 2.0f;m_RhoResolution  = 1.0f;m_MinAng         = 0;m_MaxAng         = 180;m_PtsNumThr      = 3;m_DistThrSquared = 3.0f * 3.0f;m_k              = 0.0;m_b              = 0.0;////Rho Range need to specifiedm_MinRho         = 0;m_MaxRho         = 0;m_vCosSinTab.clear();SPointXY tempPt;const float factor = CV_PI * m_AngResolution / 180.0f;for ( int i = 0; i < 180.0f / m_AngResolution; i++ ){tempPt.x = cos( i * factor );tempPt.y = sin( i * factor );m_vCosSinTab.push_back( tempPt );}m_vk.clear();m_vb.clear();}int Initial( float MinRho, float MaxRho, const vector<SPointXY> & vOrigPoints );int HoughTransform( vector<SPointXY> & vLinePoints, vector<SPointXY> & vOtherPoints );int MultiLineHT( vector< vector<SPointXY> > & vvLinePoints );int Test();protected:private:};#endif

#include "CHough.h"int CHough::Initial( float MinRho, float MaxRho, const vector<SPointXY> & vOrigPoints ){m_MinRho      = MinRho;m_MaxRho      = MaxRho;m_vOrigPoints = vOrigPoints;return 0;}int CHough::HoughTransform( vector<SPointXY> & vLinePoints, vector<SPointXY> & vOtherPoints ){vLinePoints.clear();vOtherPoints.clear();int AngNum = cvRound( ( m_MaxAng - m_MinAng ) / m_AngResolution );int RhoNum = cvRound( ( m_MaxRho - m_MinRho ) / m_RhoResolution );////注意：这里vAccum实质是二维数组，行代表Angle，列代表radius，////在m*n矩阵外再包上一层0元素，是为了下一步寻找局部极小值时比较操作更方便vector<int> vAccum( ( AngNum + 2 ) * ( RhoNum + 2 ), 0 );vector<int> vSort_buf;  //// stage 1. fill accumulatorfloat CosA = 0.0f, SinA = 0.0f, tempAng = 0.0f;int AngIdx = 0, rho = 0, ang = 0 ;for ( vector<SPointXY>::size_type i = 0; i != m_vOrigPoints.size(); ++i ){for( ang = 0; ang < AngNum; ang++ ){tempAng = m_MinAng + ang * m_AngResolution;// AngIdx = int(tempAng / 0.1f );AngIdx = int( tempAng / m_AngResolution );CosA    = m_vCosSinTab[AngIdx].x;SinA    = m_vCosSinTab[AngIdx].y;rho = cvRound( ( ( m_vOrigPoints[i].x * CosA + m_vOrigPoints[i].y * SinA ) - m_MinRho ) / m_RhoResolution ) ;//r += (RhoNum - 1) / 2;             //r有可能是负的，所以要加上(numrho - 1) / 2;vAccum[ (ang + 1) * (RhoNum + 2) + rho + 1 ]++;}}//// stage 2. find local maximums，四邻域比较for( int rho = 0; rho < RhoNum; rho++ ){for( int ang = 0; ang < AngNum; ang++ ){int base = (ang + 1) * (RhoNum + 2) + rho + 1;int Counter = vAccum[base];if( vAccum[base] > m_PtsNumThr &&vAccum[base] > vAccum[base - 1] && vAccum[base] >= vAccum[base + 1] &&vAccum[base] > vAccum[base - RhoNum - 2] && vAccum[base] >= vAccum[base + RhoNum + 2] ){vSort_buf.push_back(base);}}}if ( vSort_buf.size() == 0 )   ////没有检测到直线{return -1;}int MaxIdx = 0;for ( int i = 1; i < vSort_buf.size(); i++ ){if ( vAccum[ vSort_buf[i] ] > vAccum[ vSort_buf[MaxIdx] ] ){MaxIdx = i;}}//// stage 4. 反解直线, x * cos(sita) + y * sin(sita) = d;//int base = (n + 1) * (RhoNum + 2) + rho + 1;int idx = vSort_buf[MaxIdx];ang = idx / ( RhoNum + 2 ) - 1;rho = idx % ( RhoNum + 2 ) - 1;float d = m_MinRho + rho * m_RhoResolution;float sita = m_MinAng + ang * m_AngResolution;sita = sita * CV_PI /180;float eps = 1e-6;////if sita is very small, we will treat sin(sita) = eps;if ( fabs( sin(sita) ) < eps ){m_k = - 1.0 / eps;m_b = d / eps;}else{m_k = - cos(sita) / sin(sita);m_b = d / sin(sita);}////找出与 y = kx + b距离较近的点，认为这些点在同一条直线上float dist = 0.0f;for ( int i = 0; i != m_vOrigPoints.size(); ++i ){dist = dist_p2l_Square( m_vOrigPoints[i], m_k, m_b );if ( dist <= m_DistThrSquared ){vLinePoints.push_back( m_vOrigPoints[i] );}else{vOtherPoints.push_back( m_vOrigPoints[i] );}}return 0;}//// This Hough Transform can calculate multiple lines( k and b is different).int CHough::MultiLineHT( vector< vector<SPointXY> > & vvLinePoints ){#define TEST_MULTILINEHT#ifdef TEST_MULTILINEHTIplImage * pImage = cvCreateImage( cvSize( 270, 270 ), IPL_DEPTH_8U, 3 );cvSetZero( pImage );char WinName[200] = " Test Multi-line Hough Transform ";cvNamedWindow( WinName, CV_WINDOW_AUTOSIZE );SPointXY tempPt;for ( int i = 0; i != m_vOrigPoints.size(); i++ ){tempPt = m_vOrigPoints[i];cvCircle( pImage, cvPoint( tempPt.x, tempPt.y ), 2, RED, -1 );}cvShowImage( WinName, pImage );cvWaitKey(0);#endifvvLinePoints.clear();m_vk.clear();m_vb.clear();vector<SPointXY> vOtherPoints, vLinePoints;while ( true ){HoughTransform( vLinePoints, vOtherPoints );if ( vLinePoints.empty() ){break;}#ifdef TEST_MULTILINEHTDrawLine( pImage, m_k, m_b, GREEN );cvShowImage( WinName, pImage );cvWaitKey(0);#endifvvLinePoints.push_back( vLinePoints );m_vk.push_back( m_k );m_vb.push_back( m_b );m_vOrigPoints = vOtherPoints;}#ifdef TEST_MULTILINEHTcvShowImage( WinName, pImage );cvWaitKey(0);#endifreturn 0;}int CHough::Test(){vector<SPointXY> vOrigPoints;for ( int i = 0; i < 30; i++ ){vOrigPoints.push_back( SPointXY( 10 * i + 3, 35 * i + 10 ) );}for ( int i = 0; i < 70; i++ ){vOrigPoints.push_back( SPointXY(  i, i + 50 ) );}for ( int i = 0; i != 100; i++ ){SPointXY tempPt;tempPt.x = rand() * 400.0 / ( RAND_MAX + 1 );tempPt.y = rand() * 400.0 / ( RAND_MAX + 1 );vOrigPoints.push_back( tempPt );}Initial( -400, 400, vOrigPoints );vector< vector<SPointXY> > vvLinePoints;MultiLineHT( vvLinePoints );cout<<"There are "<<vvLinePoints.size()<<" Lines:"<<endl<<"Parameter(k, b ): "<<endl;for ( int i = 0; i != m_vk.size(); i++ ){cout<<m_vk[i]<<"  "<<m_vb[i]<<endl;}return 0;}

检测的结果：