浅谈霍夫变换

最新推荐文章于 2024-07-21 23:48:51 发布

Lucas_Wang_Robot

最新推荐文章于 2024-07-21 23:48:51 发布

阅读量821

点赞数

分类专栏： Computer Vision 文章标签：图像处理霍夫变换 C++ 应用

本文链接：https://blog.csdn.net/lump_night5/article/details/38612453

版权

Computer Vision 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

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

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

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

这个任务对一般人来说太容易了，然而要让电脑顺利完成检测直线的任务，却是个头疼的问题。我们可以分析一下，图中点只有一部分在直线上，其他的点离直线很远，这些无关的点可以看成噪声数据。分成两种情况：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 specified
		m_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 accumulator
	float 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_MULTILINEHT
	IplImage * 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);
#endif
	vvLinePoints.clear();
	m_vk.clear();
	m_vb.clear();
	vector<SPointXY> vOtherPoints, vLinePoints;
	while ( true )
	{
		HoughTransform( vLinePoints, vOtherPoints );
		if ( vLinePoints.empty() )
		{
			break;	
		}
#ifdef TEST_MULTILINEHT
		DrawLine( pImage, m_k, m_b, GREEN );
		cvShowImage( WinName, pImage );
		cvWaitKey(0);
#endif
		vvLinePoints.push_back( vLinePoints );
		m_vk.push_back( m_k );
		m_vb.push_back( m_b );
		m_vOrigPoints = vOtherPoints;
	}
#ifdef TEST_MULTILINEHT
	cvShowImage( WinName, pImage );
	cvWaitKey(0);
#endif
	return 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;
}

检测的结果：