【计算机视觉】霍夫线/圆变换从原理到源码详解

【计算机视觉】霍夫线/圆变换从原理到源码详解

---恢复内容开始---

@

1 简述

　　霍夫变换是一个经典的特征提取技术，本文主要说的是霍夫线/圆变换，即从图像中获取直线与圆，同时需要对图像进行二值化操作，效果如下。

　　霍夫变换目的通过投票程序在特定类型的形状内找到对象的不完美示例，这个投票程序是在一个参数空间中进行的，在这个参数空间中，候选对象被当作所谓的累加器空间中的局部最大值来获得，所述累加器空间由用于计算霍夫变换的算法明确的构建，霍夫变换不仅能够识别出图像中有无需要检测的图形，而且能够定位到该图像（包括位置、角度等），有标准霍夫变换，多尺度霍夫变换，统计概率霍夫变换3种较为经典的方法，主要优点是能容忍特征边界描述中的间隙，并且相对不受图像噪声的影响。

2 标准霍夫线变换原理

2.1 霍夫变换直线的方程

　　直线的方程形式有截距式，斜截式，一般式，点斜式等，但是这些形式都有“奇异”情况，其中一个截距为0，斜率为0， c=0（通常把c除到等式左面，减少一个参数的搜索）,另外某些形式的参数空间不是闭的，比如斜截式的斜率k，取值范围从0到无穷大，给量化搜索带来了困难。所以在霍夫线变换中采用了如下方程的形式：
\[p = xcos\theta + ysin\theta\]
其中\(p\)与\(\theta\)都是极坐标下点的极径与极角（也可以理解\(p\)是原点离直线最近的距离），但并不是极坐标方程，因为还在笛卡尔坐标下表示，两个参数的意义如下图：

　　这个方程的巧妙之处是角度与距离这两个参数都是有界的，而且没有奇异的情况，因此可以在参数空间中进行搜索。同时，直线霍夫变换有2个参数通过这一个等式相关联，所以程序在投票阶段只需要遍历其中一个，通过等式获取另一个参数从而进行投票与搜索峰值。

2.2 霍夫空间

　　对于直线来说，可以将图像的每一条直线与一对参数\((p, \theta)\)相关联，这个参数\((p, \theta)\)平面有时被称为霍夫空间，用于二维直线的集合，因此，对于某一指定的点，不同的\((p, \theta)\)代表了不同的直线，我们把\((p, \theta)\)之间的关系通过图像表示出来——固定一个点(3,4),角度取[0,2\(\pi\)]时，\(r\)的范围取值情况：

会发现曲线的形状类似一条正弦曲线，通常\(p\)越大，正弦曲线的振幅越大，反而就会越小。

2.3 检测直线的方法

　　可以得到一个结论，给定平面中的单个点，那么通过该点的所有直线的集合对应于(r,θ)平面中的正弦曲线，这对于该点是独特的。一组两个或更多点形成一条直线将产生在该线的(r,θ)处交叉的一条或多条正弦曲线。因此，检测共线点的问题可以转化为找到并发曲线的问题。
　　具体的，通过将霍夫参数空间量化为有限间隔或累加器单元来实现变换。随着算法的运行，每个算法都把\((x_{i}, y_{i})\)转换为一个离散化的 \((r,θ)\)曲线，并且沿着这条曲线的累加器单元被递增。累加器阵列中产生的峰值表示图像中存在相应的直线的有力证据。 换句话说，将每个交点看成一次投票，也就是说\(A(r,θ)=A(r,θ)+1\)，所有点都如此进行计算后，可以设置一个阈值，投票大于这个阈值的可以认为是找到的直线。

2.4 一个例子

　　举一个具体的例子，来体现霍夫变换的流程同时更深的理解一下该算法（此例子摘取某博客，下文会给出所有参考博客的链接）

　　霍夫变换可用于识别最适合一组给定边缘点的曲线的参数。该边缘描述通常从诸如Roberts Cross，Sobel或 Canny边缘检测器的特征检测算子获得，并且可能是嘈杂的，即其可能包含对应于单个整体特征的多个边缘片段。此外，由于边缘检测器的输出仅限定图像中的特征的位置，所以霍夫变换进一步是确定两个特征是什么（即检测其具有参数（或其他）的特征描述）以及 它们有多少个存在于图像中。
　　为了详细说明霍夫变换，我们从两个闭合矩形的简单图像开始。

使用 Canny边缘检测器可产生一组边界描述的这个部分，如下图：

　　这里我们看到了图像中的整体边界，但是这个结果并没有告诉我们边界描述中的特征的身份（和数量）。在这种情况下，我们可以使用Hough（线检测）变换来检测该图像的八个单独的直线段，从而确定对象的真实几何结构。
　　如果我们使用这些边缘/边界点作为Hough变换的输入，每一个点通过遍历不同的\((r, \theta)\)都会在笛卡尔空间中生成一条曲线，当被视为强度图像时， 累加器阵列看起来像如下所示：

可以利用直方图均衡技术使得图像可以让我们看到包含在低强度像素值中的信息模式，如下图：

　　注意，虽然r​和θ​是概念上的极坐标，但是累加器空间以横坐标θ​和纵坐标r​的矩形绘制 。请注意，累加器空间环绕图像的垂直边缘，因此实际上只有8个实际峰值。
　　由梯度图像中的共线点生成的曲线(r,θ) 在霍夫变换空间中的峰中相交。这些交点表征原始图像的直线段。有许多方法可用于从累加器阵列中提取这些亮点或局部最大值。例如，一个简单的方法涉及阈值处理，然后 对累加器阵列图像中孤立的亮点集群进行一些细化处理。这里我们使用相对阈值来提取(r,θ) ，对应于原始图像中的每条直线边缘的点。（换句话说，我们只采用累加器数组中的那些局部最大值，其值等于或大于全局最大值的某个固定百分比。）
  从Hough变换空间（即，反霍夫变换）映射回笛卡尔空间产生一组图像主题的线描述。通过将该图像叠加在原始的反转版本上，我们可以确认霍夫变换找到两个矩形的8个真实边的结果，并且因此展示出了遮挡场景的基础几何形状。

　　请注意，在这个简单的例子中，检测到的和原始图像线的对齐精度显然不完美，这取决于累加器阵列的量化。（还要注意许多图像边缘有几条检测线，这是因为有几个附近的霍夫空间峰值具有相似的线参数值，存在控制这种效应的技术，但这里没有用来说明标准霍夫变换）
　　还要注意，由霍夫变换产生的线条的长度是无限的。如果我们希望识别产生变换参数的实际线段，则需要进一步的图像分析以便查看这些无限长线的哪些部分实际上具有点。
　　为了说明Hough技术对噪声的鲁棒性，Canny边缘描述已被破坏1％（由椒盐噪声引起）， 在Hough变换之前，如图12所示。

霍夫空间如下图所示：

将这个结果反霍夫变换（相对阈值设置为40%，并将它覆盖在原来的结果上）

可以通过变换图像来研究Hough变换对特征边界中的间隙的敏感性，如下图：

然后我们再将其变换到霍夫变换空间，如下图

然后使用40%的相对阈值去对图像反霍夫变换（同样也是叠加在原图上）

  在这种情况下，因为累加器空间没有接收到前面例子中的条目数量，只能找到7个峰值，但这些都是结构相关的线。

3 霍夫线变换的算法流程

3.1 标准霍夫线变换算法流程

1. 读取原始图像，并转换成灰度图，利用阈值分割或者边缘检测算子转换成二值化边缘图像
2. 初始化霍夫空间， 令所有\(Num(\theta, p) = 0\)
3. 对于每一个像素点\((x, y)\)，在参数空间中找出所有满足\(xcos\theta + ysin\theta = p\)的\((\theta, p)\)对,然后令\(Num(\theta, p) = Num(\theta, p) + 1\)
4. 统计所有\(Num(\theta, p)\)的大小，取出\(Num(\theta, p) > τ\)的参数（\(τ\)是所设的阈值），从而得到一条直线。

5. 将上述流程取出的直线，确定与其相关线段的起始点与终止点（有一些算法，如蝴蝶形状宽度，峰值走廊之类）

   3.2 统计概率霍夫变换算法流程

   标准霍夫变换本质上是把图像映射到它的参数空间上，它需要计算所有的M个边缘点，这样它的运算量和所需内存空间都会很大。如果在输入图像中只是处理\(m(m<M)\)个边缘点，则这m个边缘点的选取是具有一定概率性的，因此该方法被称为概率霍夫变换（Probabilistic Hough Transform）。该方法还有一个重要的特点就是能够检测出线端，即能够检测出图像中直线的两个端点，确切地定位图像中的直线。
   HoughLinesP函数就是利用概率霍夫变换来检测直线的。它的一般步骤为：

6. 随机抽取图像中的一个特征点，即边缘点，如果该点已经被标定为是某一条直线上的点，则继续在剩下的边缘点中随机抽取一个边缘点，直到所有边缘点都抽取完了为止；
7. 对该点进行霍夫变换，并进行累加和计算；
8. 选取在霍夫空间内值最大的点，如果该点大于阈值的，则进行步骤4，否则回到步骤1；
9. 根据霍夫变换得到的最大值，从该点出发，沿着直线的方向位移，从而找到直线的两个端点；

10. 计算直线的长度，如果大于某个阈值，则被认为是好的直线输出，回到步骤1。

4 OpenCV中的霍夫变换函数

4.1标准霍夫变换HoughLines()函数

void HoughLines(InputArray image,OutputArray lines, double rho, 
double theta, int threshold, double srn=0,double stn=0 )

第一个参数，InputArray类型的image，输入图像，即源图像，需为8位的单通道二进制图像，可以将任意的源图载入进来后由函数修改成此格式后，再填在这里。

第二个参数，InputArray类型的lines，经过调用HoughLines函数后储存了霍夫线变换检测到线条的输出矢量。每一条线由具有两个元素的矢量\((p, \theta)\)表示，其中，\(p\)是离坐标原点((0,0)（也就是图像的左上角）的距离。 \(\theta\)是弧度线条旋转角度（0~垂直线，π/2~水平线）。

第三个参数，double类型的rho，以像素为单位的距离精度。另一种形容方式是直线搜索时的进步尺寸的单位半径。PS:Latex中/rho就表示 \(\rho\)。

第四个参数，double类型的theta，以弧度为单位的角度精度。另一种形容方式是直线搜索时的进步尺寸的单位角度。

第五个参数，int类型的threshold，累加平面的阈值参数，即识别某部分为图中的一条直线时它在累加平面中必须达到的值。大于阈值threshold的线段才可以被检测通过并返回到结果中。

第六个参数，double类型的srn，有默认值0。对于多尺度的霍夫变换，这是第三个参数进步尺寸rho的除数距离。粗略的累加器进步尺寸直接是第三个参数rho，而精确的累加器进步尺寸为rho/srn。

第七个参数，double类型的stn，有默认值0，对于多尺度霍夫变换，srn表示第四个参数进步尺寸的单位角度theta的除数距离。且如果srn和stn同时为0，就表示使用经典的霍夫变换。否则，这两个参数应该都为正数。

4.2 统计概率霍夫变换(HoughLinesP)

C++: 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，允许将同一行点与点之间连接起来的最大的距离。

5 源码分析

此部分参考几位博主的博客，更改了里面一些我认为不太对的地方以及加上了一些自己的理解还有就是一些没有解释的地方

5.1 HoughLines()源码分析

阅读源码之后，发现源码确实细节方面做得比较好，无论是内存方面，还是效率方面，包括一些小细节，比如避免除法，等等。所以源码值得借鉴与反复参考，具体解释如下。

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的空间范围， 这里以周长作为rho的上界

/*
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;    //因为theta是从0到π的，所以cos(theta)是有负的，所以就所有的r += 最大值的一半，让极径都>0
                    //另一方面，这也预防了r小于0，让本属于这一行的点跑到了上一行，这样可能导致同一个位置却表示了多个不同的“点”


                    //r表示的是距离，n表示的是角点，在累加器内找到它们所对应的位置（即霍夫空间内的位置），其值加1
                    accum[(n+1) * (numrho+2)+ r + 1]++;
                    /*
                    将1维数组想象成2维数组
                    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];   //找到索引(下标)
        //分离出该极大值在霍夫空间中的位置
        //因为n是从0开始的，而之前为了防止越界，所以将所有的n+1了，因此下面要-1，同理r
        int n = cvFloor(idx*scale) - 1;   //向下取整
        int r = idx - (n+1)*(numrho+2) - 1;
        //最终得到极大值所对应的角度和距离
        line.rho = (r - (numrho - 1)*0.5f)*rho; //因为之前统一将r+= (numrho-1)/2， 因此需要还原以获得真实的rho
        line.angle = n * theta;
        //存储到序列内
        cvSeqPush( lines, &line );  //用序列存放多条直线
    }
}

5.2 HoughLinesP源码分析

这个代码与之前HoughLines风格完全不同，个人感觉没有之前HoughLines的优美，同时小优化也并不是特别注意，不过实现算法思想的大体思路还是值得我这个初学者借鉴的,但是里面那个shitf操作没看懂，不是很理解，如果有同学会的话，希望解释一下

static void
icvHoughLinesProbabilistic( CvMat* image,
                            float rho, float theta, int threshold,
                            int lineLength, int lineGap,
                            CvSeq *lines, int linesMax )
{
    //accum为累加器矩阵，mask为掩码矩阵
    cv::Mat accum, mask;
    cv::vector<float> trigtab;    //用于存储事先计算好的正弦和余弦值
    //开辟一段内存空间
    cv::MemStorage storage(cvCreateMemStorage(0));
    //用于存储特征点坐标，即边缘像素的位置
    CvSeq* seq;
    CvSeqWriter writer;
    int width, height;    //图像的宽和高
    int numangle, numrho;    //角度和距离的离散数量
    float ang;
    int r, n, count;
    CvPoint pt;
    float irho = 1 / rho;    //距离分辨率的倒数
    CvRNG rng = cvRNG(-1);    //随机数
    const float* ttab;    //向量trigtab的地址指针
    uchar* mdata0;    //矩阵mask的地址指针
    //确保输入图像的正确性
    CV_Assert( CV_IS_MAT(image) && CV_MAT_TYPE(image->type) == CV_8UC1 );

    width = image->cols;    //提取出输入图像的宽
    height = image->rows;    //提