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 circlesHoughCircles( src, circles, CV_HOUGH_GRADIENT, 1, dst.rows/10, 200, savevalue+10, 0, 0 );/// Draw the circles detectedprintf("%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 centercircle( cdst, center, 3, Scalar(0,255,0), -1, 8, 0 );// circle outlinecircle( 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 );  //用序列存放多条直线    }}

恍恍惚惚，结束了