OpenCV霍夫变换系列（前篇）-经典霍夫线变换

前言：最近新来的的我的大学室友（现在也是我的学弟）在研究霍夫线变换，我之前只是知道这玩意可以拿来做直线检测，并没有深入研究，那既然提到了，还是按照我们的老规矩，原理，示例以及OpenCV这一套流程走下来。

菜鸟一枚，好多写的不好，有点啰嗦，见谅

主要参考博客：

原理部分：http://www.cnblogs.com/php-rearch/p/6760683.html（相当清楚，不解释，）

源码分析1:http://blog.csdn.net/zhaocj/article/details/50281537(赵春江老师，很多源码都给出了详解，尤其是那篇sift，看的热血沸腾)

源码分析2:http://blog.csdn.net/traumland/article/details/51319644

源码分析3:http://blog.csdn.net/sunshine_in_moon/article/details/45271647

Samples1:http://blog.csdn.net/poem_qianmo/article/details/26977557/   (还是浅墨)

Samples2:http://www.cnblogs.com/skyfsm/p/6881686.html

还有一些博主未给出链接，也衷心表示感谢！

霍夫变换（Hough）

1.基本原理：

一条直线可由两个点A=(X1,Y1)和B=(X2,Y2)确定(笛卡尔坐标)

另一方面，也可以写成关于（k,q）的函数表达式（霍夫空间）：

对应的变换可以通过图形直观表示：

变换后的空间成为霍夫空间。即：笛卡尔坐标系中一条直线，对应霍夫空间的一个点。


反过来同样成立（霍夫空间的一条直线，对应笛卡尔坐标系的一个点）：

再来看看A、B两个点，对应霍夫空间的情形：

一步步来，再看一下三个点共线的情况：

可以看出如果笛卡尔坐标系的点共线，这些点在霍夫空间对应的直线交于一点：这也是必然，共线只有一种取值可能。

如果不止一条直线呢？再看看多个点的情况（有两条直线）：

其实（3，2）与（4，1）也可以组成直线，只不过它有两个点确定，而图中A、B两点是由三条直线汇成，这也是霍夫变换的后处理的基本方式：选择由尽可能多直线汇成的点。

看看，霍夫空间：选择由三条交汇直线确定的点（中间图），对应的笛卡尔坐标系的直线（右图）。

到这里问题似乎解决了，已经完成了霍夫变换的求解，但是如果像下图这种情况呢？

k=∞是不方便表示的，而且q怎么取值呢，这样不是办法。因此考虑将笛卡尔坐标系换为：极坐标表示。

备注：很多文章一来就是极坐标，感觉没有这样看着顺畅，舒服。

在极坐标系下，其实是一样的：极坐标的点--->霍夫空间的直线，这个地方要注意，这个必须注意，只有垂直的时候才可能是

一一对应关系。

只不过霍夫空间不再是[k,q]的参数，而是[ρ,θ]的参数，给出对比图：

是不是就一目了然了？

有一个离散化的过程，本质上就是这样（这张图太深动形象了，怪不得后面源码分析就有博主总是在说格子数）

交点怎么求解呢？细化成坐标形式，取整后将交点对应的坐标进行累加，最后找到数值最大的点就是求解的[ρ,θ]，也就求解出了直线。(OpenCV不是这样做的，但是也用到了累加的思想)

下面给出OpenCV中HoughLines的霍夫变换的直线检测步骤:

1、对边缘图像进行霍夫空间变换；

2、在4邻域内找到霍夫空间变换的极大值；

3、对这些极大值安装由大到小顺序进行排序，极大值越大，越有可能是直线；

这里说明一下极大值的含义，指的是霍夫空间中相交于某一点的曲线的数目的极大值，不要理解为数学上某一条曲线的极大值
4、输出直线。

2. Samples

OpenCV中的霍夫直线检测的函数为HoughLines

还有一个改进版本的HoughLinesP函数（统计概论霍夫直线检测）

函数原型分别如下：

//函数HoughLines的原型为：
void HoughLines(InputArray image,OutputArray lines, double rho, double theta, int threshold, double srn=0,double stn=0 )
/*
image为输入图像，要求是单通道的二值图像
lines为输出直线向量，两个元素的向量（ρ，θ）代表一条直线，ρ是从原点（图像的左上角）的距离，θ是直线的角度（单位是弧度），0表示垂直线，π/2表示水平线
rho为距离分辨率
theta为角度分辨率
threshold为阈值，在步骤2中，只有大于该值的点才有可能被当作极大值，即至少有多少条正弦曲线交于一点才被认为是直线
srn和stn在多尺度霍夫变换的时候才会使用，在这里我们只研究经典霍夫变换的源码
 */

//函数HoughLinesP的原型为：
void HoughLinesP(InputArray image, OutputArray lines, double rho, double theta, int threshold, double minLineLength=0, double maxLineGap=0 )
/*
第一个参数，InputArray类型的image，输入图像，即源图像，需为8位的单通道二进制图像，可以将任意的源图载入进来后由函数修改成此格式后，再填在这里。
第二个参数，InputArray类型的lines，经过调用HoughLinesP函数后后存储了检测到的线条的输出矢量，每一条线由具有四个元素的矢量(x_1,y_1, x_2, y_2）  表示，其中，(x_1, y_1)和(x_2, y_2) 是是每个检测到的线段的结束点。
第三个参数，double类型的rho，以像素为单位的距离精度。另一种形容方式是直线搜索时的进步尺寸的单位半径。
第四个参数，double类型的theta，以弧度为单位的角度精度。另一种形容方式是直线搜索时的进步尺寸的单位角度。
第五个参数，int类型的threshold，累加平面的阈值参数，即识别某部分为图中的一条直线时它在累加平面中必须达到的值。大于阈值threshold的线段才可以被检测通过并返回到结果中。
第六个参数，double类型的minLineLength，有默认值0，表示最低线段的长度，比这个设定参数短的线段就不能被显现出来。
第七个参数，double类型的maxLineGap，有默认值0，允许将同一行点与点之间连接起来的最大的距离。
*/

由于函数HoughLines只输出直线所对应的角度和距离，所以在进行直线检测的时候还要把其转换为直角坐标系下的数据， 另外输入图像还必须是边缘图像，下面就是具体的实例：

#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>
#include <iostream>
#include <Windows.h>

using namespace cv;
using namespace std;

Mat src, dst, cdst;
const char* filename = "1.jpg";
const char *winname="hough";
const char *trackbarname="1.houghlines\n2.houghlinep\n3.houghcircle";
int savevalue=100,houghtype;
const int maxthreshold=150;

void help()
{
	cout << "\nThis program demonstrates line finding with the Hough transform.\n";
}
void choisehoughlines()
{

	vector<Vec2f> lines;
	Canny(src, dst,50, 200, 3);
	cvtColor(dst, cdst, CV_GRAY2BGR);
	HoughLines(dst, lines, 1, CV_PI/180, savevalue+10, 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), 1, CV_AA);
		float angle = theta / CV_PI *180;
		//cout<<"line "<< i << ": "<<"rho: "<<rho<<"theta: "<<theta
			//<<"angle: "<< angle<<endl;
	}
	imshow(winname, cdst); 
}
void choisehoughlinep()
{
	vector<Vec4i> lines;
	Canny(src, dst,50, 200, 3);
	cvtColor(dst, cdst, CV_GRAY2BGR);
	HoughLinesP(dst, lines, 1, CV_PI/180, savevalue+10, savevalue+10, 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), 1, CV_AA);
	}
	imshow(winname, cdst);

}
void choisehoughcircle()
{
	vector<Vec3f> circles;
	cvtColor(src, cdst, CV_GRAY2BGR);
	/// Apply the Hough Transform to find the circles
	HoughCircles( src, circles, CV_HOUGH_GRADIENT, 1, dst.rows/10, 200, savevalue+10, 0, 0 );
	/// Draw the circles detected
	printf("%d",circles.size());
	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 center
		circle( cdst, center, 3, Scalar(0,255,0), -1, 8, 0 );
		// circle outline
		circle( cdst, center, radius, Scalar(0,0,255), 3, 8, 0 );
	}
	imshow(winname, cdst);
}
void choice(int,void *)
{
	switch (houghtype)
	{
	case 0:choisehoughlines();break;
	case 1:choisehoughlinep();break;
	//case 2:choisehoughcircle();break;
	}
}
int main(int argc, char** argv)
{
	SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN |FOREGROUND_INTENSITY);
	//system("color 3E");
	help();
	src = imread(filename, 0);
	if(src.empty())
	{
		help();
		cout << "can not open " << filename << endl;
		return -1;
	}
	namedWindow(winname);
	createTrackbar(trackbarname,winname,&houghtype,2,choice);
	createTrackbar("thresholdvalue",winname,&savevalue,maxthreshold,choice);
	imshow("source", src);
	choice(0,0);
	while(char(waitKey())!= 'q');

	return 0;
}

效果如下：

这是一个综合示例，我把霍夫圆变换注释掉了，那个放在下一篇单独讲解，其中输出经典霍夫线变换的角度的代码我也注释掉了，

需要的可以自己加上。

对于经典的线变换计算得到的两点的坐标为（ρcosθ-1000sinθ，ρsinθ+1000cosθ），（ρcosθ+1000sinθ，ρsinθ-1000cosθ），我个人觉得

计算过程如下：

结论：

houghlines的计算效率比较低O(n*n*m)，耗时较长，而且没有检测出直线的端点。
统计概论霍夫直线检测houghlinesP是一个改进，不仅执行效率较高，而且能检测到直线的两个端点。
思想：
先随机检测出一部分直线，然后将直线上点的排查掉，再进行其他直线的检测
1. 首先仅统计图像中非零点的个数，对于已经确认是某条直线上的点就不再变换了。
2. 对所以有非零点逐个变换到霍夫空间
a. 并累加到霍夫统计表(图像)中，并统计最大值
b. 最大值与阈值比较，小于阈值，则继续下一个点的变换
c. 若大于阈值，则有一个新的直线段要产生了
d. 计算直线上线段的端点、长度，如果符合条件，则保存此线段，并mark这个线段上的点不参与其他线段检测的变换。

3.源码分析

HoughLines函数是在sources/modules/imgproc/src/hough.cpp文件中定义的：

void cv::HoughLines( InputArray _image, OutputArray _lines,
                     double rho, double theta, int threshold,
                     double srn, double stn )
{
    //申请一段内存，用于存储霍夫变换后所检测到的直线
    Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
    //提取输入图像矩阵
    Mat image = _image.getMat();
    CvMat c_image = image;
    //调用由C语音写出的霍夫变换的函数
    CvSeq* seq = cvHoughLines2( &c_image, storage, srn == 0 && stn == 0 ?
                    CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,
                    rho, theta, threshold, srn, stn );
    //把由cvHoughLines2函数得到的霍夫变换直线序列转换为数组
    seqToMat(seq, _lines);
}

cvHoughLines2函数是霍夫变换直线检测的关键，函数的输出为所检测到的直线序列，它的第3个形参“srn == 0 && stn == 0 ? CV_HOUGH_STANDARD :CV_HOUGH_MULTI_SCALE”表示的是该霍夫变换是经典霍夫变换还是多尺度霍夫变换，它是由变量srn和stn决定的，只要这两个变量有一个不为0，就进行多尺度霍夫变换，否则为经典霍夫变换。另外cvHoughLines2函数不仅可以用于经典霍夫变换和多尺度霍夫变换，还可以用于概率霍夫变换。

CV_IMPL CvSeq*
cvHoughLines2( CvArr* src_image, void* lineStorage, int method,
               double rho, double theta, int threshold,
               double param1, double param2 )
{
    CvSeq* result = 0;

    CvMat stub, *img = (CvMat*)src_image;
    CvMat* mat = 0;
    CvSeq* lines = 0;
    CvSeq lines_header;
    CvSeqBlock lines_block;
    int lineType, elemSize;
    int linesMax = INT_MAX;    //输出最多直线的数量，设为无穷多
    int iparam1, iparam2;

    img = cvGetMat( img, &stub );
    //确保输入图像是8位单通道
    if( !CV_IS_MASK_ARR(img))
        CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
    //内存空间申请成功，以备输出之用
    if( !lineStorage )
        CV_Error( CV_StsNullPtr, "NULL destination" );
    //确保rho，theta，threshold这三个参数大于0
    if( rho <= 0 || theta <= 0 || threshold <= 0 )
        CV_Error( CV_StsOutOfRange, "rho, theta and threshold must be positive" );

    if( method != CV_HOUGH_PROBABILISTIC )    //经典霍夫变换和多尺度霍夫变换
    {
        //输出直线的类型为32位浮点双通道，即ρ和θ两个浮点型变量
        lineType = CV_32FC2;
        elemSize = sizeof(float)*2;
    }
    else    //概率霍夫变换
    {
        //输出直线的类型为32位有符号整型4通道，即两个像素点的4个坐标
        lineType = CV_32SC4;
        elemSize = sizeof(int)*4;
    }
    //判断lineStorage的类型，经分析lineStorage只可能是STORAGE类型
    if( CV_IS_STORAGE( lineStorage ))
    {
        //在lineStorage内存中创建一个序列，用于存储霍夫变换的直线，lines为该序列的指针
        lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage );
    }
    else if( CV_IS_MAT( lineStorage ))
    {
        mat = (CvMat*)lineStorage;

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

        if( CV_MAT_TYPE( mat->type ) != lineType )
            CV_Error( CV_StsBadArg,
            "The destination matrix data type is inappropriate, see the manual" );

        lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,
                                         mat->rows + mat->cols - 1, &lines_header, &lines_block );
        linesMax = lines->total;
        cvClearSeq( lines );
    }
    else
        CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );

    iparam1 = cvRound(param1);
    iparam2 = cvRound(param2);

    switch( method )
    {
    case CV_HOUGH_STANDARD:    //经典霍夫变换
          icvHoughLinesStandard( img, (float)rho,
                (float)theta, threshold, lines, linesMax );
          break;
    case CV_HOUGH_MULTI_SCALE:    //多尺度霍夫变换
          icvHoughLinesSDiv( img, (float)rho, (float)theta,
                threshold, iparam1, iparam2, lines, linesMax );
          break;
    case CV_HOUGH_PROBABILISTIC:    //概率霍夫变换
          icvHoughLinesProbabilistic( img, (float)rho, (float)theta,
                threshold, iparam1, iparam2, lines, linesMax );
          break;
    default:
        CV_Error( CV_StsBadArg, "Unrecognized method id" );
    }
    //在前面判断lineStorage类型时，已经确定它是STORAGE类型，因此没有进入else if内，也就是没有对mat赋值，所以在这里进入的是else
    if( mat )
    {
        if( mat->cols > mat->rows )
            mat->cols = lines->total;
        else
            mat->rows = lines->total;
    }
    else
        result = lines;
    //返回lines序列的指针，该序列存储有霍夫变换的直线
    return result;
}

上面这两段代码我没怎么研究，我只是分析了经典霍夫变换的icvHoughLinesStandard函数：

并且其中有一个地方还是不明白,已经标注出来，希望有人解释。

static void icvHoughLinesStandard( const CvMat* img, float rho, float theta,
                                    int threshold, CvSeq *lines, int lineMax)
{
    cv::AutoBuffer<int> _accum, _sort_buf;
    cv::AutoBuffer<float> _tabSin, _tabCos;

    const uchar* image;
    int step, width, height;
    int numangle, numrho;
    int total=0;
    float ang;
    int r, n;
    int i, j;
    float irho = 1 / rho;
    double scale;
    //再次确保输入图像的正确性
    CV_Assert( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1);

    image = img->data.ptr;    //得到图像的指针  
    step = img->step;    //得到图像的步长  
    width = img->cols;    //得到图像的宽  
    height = img->rows;    //得到图像的高  

    //由角度和距离的分辨率得到角度和距离的数量，即霍夫变换后角度和距离的个数
    numangle = cvRound(CV_PI / theta);     // 霍夫空间,角度方向的大小
    numrho = cvRound(((width + height)*2 + 1) / rho);  //r的空间范围

/*
allocator类是一个模板类，定义在头文件memory中，用于内存的分配、释放、管理，它帮助我们将内存分配和对象构造分离开来。具体地说，allocator类将内存的分配和对象的构造解耦，分别用allocate和construct两个函数完成，同样将内存的释放和对象的析构销毁解耦，分别用deallocate和destroy函数完成。
 */
    //为累加器数组分配内存空间  
    //该累加器数组其实就是霍夫空间，它是用一维数组表示二维空间
    _accum.allocate((numangle+2) * (numrho + 2)); 
    //为排序数组分配内存空间
    _sort_buf.allocate(numangle * numrho);
    //为正弦和余弦列表分配内存空间
    _tabSin.allocate(numangle);
    _tabCos.allocate(numangle);
    //分别定义上述内存空间的地址指针  
    int *accum = _accum, *sort_buf = _sort_buf;
    int *tabSin = _tabSin, *tabCos = _tabCos;
    //累加器数组清零
    memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );

    //为避免重复运算，事先计算好sinθi/ρ和cosθi/ρ
    for( ang = 0, n = 0; n < numangle; ang += theta, n++ ) //计算正弦曲线的准备工作,查表
    {
        tabSin[n] = (float)(sin(ang) * irho);
        tabCos[n] = (float)(cos(ang) * irho);
    }

    //stage 1. fill accumulator
    //执行步骤1，逐点进行霍夫空间变换，并把结果放入累加器数组内 
    for( i = 0 ; i < height; i++)
        for( j = 0; j < width; j++)
        {
            //只对图像的非零值处理，即只对图像的边缘像素进行霍夫变换
            if( image[i * step + j] != 0 )  //将每个非零点,转换为霍夫空间的离散正弦曲线,并统计。
                for( n = 0; n < numangle; n++ )
                {   
                    //根据公式: ρ = xcosθ + ysinθ
                    //cvRound()函数：四舍五入
                    r = cvRound( j * tabCos[n] + i * tabSin[n])
                    //numrho是ρ的最大值，或者说最大取值范围
                    r += (numrho - 1) / 2;    //这一步是真的看不太明白？？？
                    //(另一位博主的解释)哈哈,这里让我想了好久,为什么这样安排呢?  可能是因为θ、ρ正负问题  ，但我感觉解释不通
                    
                    
                    //r表示的是距离，n表示的是角点，在累加器内找到它们所对应的位置（即霍夫空间内的位置），其值加1                                                         
                    accum[(n+1) * (numrho+2)+ r + 1]++;
                    /* 
                    最初我也是下面这样理解的，觉得比较好理解，直观，但是是一维数组
                    哪来的行与列额!
                    n+1是为了第一行空出来
                    // numrho+2 是总共的列数,这里实际意思应该是半径的所有可能取值，加2是为了防止越界，但是说成列数我觉得也没错，并且我觉得这样比较好理解
                    //r+1 是为了把第一列空出来
                    //因为程序后面要比较前一行与前一列的数据,  这样省的越界
                    因此空出首行和首列*/
                }
        }
    // stage 2. find local maximums 
    // 执行步骤2，找到局部极大值，即非极大值抑制 
    // 霍夫空间,局部最大点,采用四领域判断,比较。(也可以使8邻域或者更大的方式),如果不判断局部最大值,同时选用次大值与最大值,就可能会是两个相邻的直线,但实际是一条直线。选用最大值,也是去除离散的近似计算带来的误差,或合并近似曲线。
    for( r = 0 ; r < numrho; r++ )
        for( 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[total++] = base;
        }

    //stage 3. sort the detected lines by accumulator value  
    //执行步骤3，对存储在sort_buf数组内的累加器的数据按由大到小的顺序进行排序
    icvHoughSortDescent32s( sort_buf, total, accum );
    /*OpenCV中自带了一个排序函数，名称为：
    void icvHoughSortDescent32s(int *sequence , int sum , int*data)，参数解释：
    第三个参数：数组的首地址
    第二个参数：要排序的数据总数目
    第一个参数：此数组存放data中要参与排序的数据的序号
    而且这个排序算法改变的只是sequence[]数组中的元素，源数据data[]未发生丝毫改变。
    */

    // stage 4. store the first min(total, linesMax ) lines to the output buffer,输出直线
    lineMax = MIN(lineMax, total);  //linesMax是输入参数，表示最多输出几条直线
    //事先定义一个尺度
    scale = 1./(numrho+2);
    for(i=0; i<linesMax; i++)    // 依据霍夫空间分辨率,计算直线的实际r,theta参数
    {
        //CvLinePolar 直线的数据结构
        //CvLinePolar结构在该文件的前面被定义
        CvLinePolar line;
        //idx为极大值在累加器数组的位置    
        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 );  //用序列存放多条直线
    }
}

恍恍惚惚，结束了