opencv 霍夫曼变换

霍夫变换不仅可以找出图片中的直线,也可以找出圆,椭圆,三角形等等,只要你能定义出直线方程,圆形的方程等等.

不得不说,现在网上的各种博客质量真的不行,网上一堆文章,乱TM瞎写,误人子弟.本身自己就没有理解的很清楚,又不去读算法实现的源码,写的云山雾罩的,越看越懵逼.

霍夫变换本身的思路是很简明的.这篇文章我们就以霍夫直线变换说明算法的思想.

霍夫变换
思考一下,二维平面里我们怎么表达直线.

有两种表达方式:

直角坐标系(也叫笛卡尔坐标系)
极坐标系(也叫球坐标系)

第一种就是最常见的直角坐标系下的表达:y=ax+b的形式.
第二种就是极坐标系下的表达:
我们把直角坐标系下的直线方程用r,theta去表达直线方程的斜率和截距.

则得到极坐标下的表达: r=xcosθ+ysinθ

假设图像中某像素点坐标为(x,y).在直角坐标系下穿过这一点我们可以画出无数条直线.
转化到一个r-θ坐标系下,我们就可以绘制出一条曲线.也就是r=xcosθ+ysinθ中的x,y是已知数,θ和r是未知数

这条曲线上每一个θ对应一个r,代表了一条直线.这些直线的共同点是他们都穿过了坐标为(x,y)的像素点.

针对图像中的每一个像素点,我都可以绘制出一条曲线来表达穿过该点的无数条直线. 那曲线的交点代表什么呢? 很显然,代表着交点处的(θ,r)所代表的直线即穿过了像素点A,又穿过了像素点B,像素点C…

怎么样叫做"找到图中的一条直线"
回到我们的问题,我们想找出图像中的一条线.意味着什么?
很多博客说了,意味着找出一条直线,尽可能多地穿过各个像素点.

我TM随便在图像上画直线,不都能穿过很多像素点吗?
实际上,应该是找出一条直线尽可能多地穿过"有效像素点".这也是为什么霍夫变换前一定要先做边缘检测的原因.经过canny检测以后(不知道的参考上一篇文章),得到的图像矩阵,只有在边缘处其像素灰度值才是比较大的,反映在图像上就是白色亮点,在非边缘处,其灰度值是0,反映在图像上就是黑色.这些代表了边缘的像素点就是有效像素点.

即:假如我能找到这么一条直线,穿过了很多个有效像素点(这个就是我们需要调参的阈值),那我就说我在图像中找到了一条直线. . 同理,找圆,找三角形还是找任意形状都是一个道理.

比方说,下面这个图

你就找不到一条直线,穿过很多个白点.所以图中是不存在直线的.

霍夫变换的过程
canny边缘检测提取出边缘
对边缘图像中的每个像素点,
伪代码如下

for (every pixel)
{
    if(pixel is effective edge pixel)
    {
        for(int theta = 0; theta < 360; theta++)
        {
            r=xcosθ+ysinθ;//x,y为pixel坐标
            accum(theta,r) += 1; //(theta,r)所代表的直线经过的像素点数量加1
        }
    }
}

for(every element in accum)
{
    if (count of (theta,r) > thershold)
    {
        find line (theta,r)
    }
}

opencv示例
houghlines api

其中, double rho, double theta,决定了最终有多少种(theta,r)的组合.决定了过每个像素点的线的可能情况.这个值越小,粒度就越细,需要的计算量也越大. 一般取rho=1,即1像素.theta取1度.
下面是一个提取车位图片中直线的示例

import sys
import math
import cv2 as cv
import numpy as np
def test():
    src = cv.imread("/home/sc/disk/keepgoing/opencv_test/houghtest.jpg")
    src = cv.GaussianBlur(src, (3, 3), 0)
    gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
    
    dst = cv.Canny(src, 150, 300, None, 3)
    lines = cv.HoughLines(dst, 1, np.pi / 180, 150, None, 0, 0)
    
    # Copy edges to the images that will display the results in BGR
    cdst = cv.cvtColor(dst, cv.COLOR_GRAY2BGR)
    cdstP = np.copy(cdst)
    
    lines = cv.HoughLines(dst, 1, np.pi / 180, 200, None, 0, 0)
    
    if lines is not None:
        for i in range(0, len(lines)):
            rho = lines[i][0][0]
            theta = lines[i][0][1]
            a = math.cos(theta)
            b = math.sin(theta)
            x0 = a * rho
            y0 = b * rho
            pt1 = (int(x0 + 1000*(-b)), int(y0 + 1000*(a)))
            pt2 = (int(x0 - 1000*(-b)), int(y0 - 1000*(a)))
            cv.line(cdst, pt1, pt2, (0,0,255), 3, cv.LINE_AA)
    
    
    cv.imshow("origin",src)
    cv.imshow("dst1",dst)
    cv.imshow("dst2",cdst)
    if 27 == cv.waitKey():
        cv.destroyAllWindows()

test()

opencv源码解读
opencv 官方实现

static void
HoughLinesStandard( InputArray src, OutputArray lines, int type,
                    float rho, float theta,
                    int threshold, int linesMax,
                    double min_theta, double max_theta )
{
    CV_CheckType(type, type == CV_32FC2 || type == CV_32FC3, "Internal error");

    Mat img = src.getMat();

    int i, j;
    float irho = 1 / rho;

    CV_Assert( img.type() == CV_8UC1 );
    CV_Assert( linesMax > 0 );

    const uchar* image = img.ptr();
    int step = (int)img.step;
    int width = img.cols;
    int height = img.rows;

    int max_rho = width + height;
    int min_rho = -max_rho;

    CV_CheckGE(max_theta, min_theta, "max_theta must be greater than min_theta");

    int numangle = cvRound((max_theta - min_theta) / theta);
    int numrho = cvRound(((max_rho - min_rho) + 1) / rho);

#if defined HAVE_IPP && IPP_VERSION_X100 >= 810 && !IPP_DISABLE_HOUGH
    if (type == CV_32FC2 && CV_IPP_CHECK_COND)
    {
        IppiSize srcSize = { width, height };
        IppPointPolar delta = { rho, theta };
        IppPointPolar dstRoi[2] = {{(Ipp32f) min_rho, (Ipp32f) min_theta},{(Ipp32f) max_rho, (Ipp32f) max_theta}};
        int bufferSize;
        int nz = countNonZero(img);
        int ipp_linesMax = std::min(linesMax, nz*numangle/threshold);
        int linesCount = 0;
        std::vector<Vec2f> _lines(ipp_linesMax);
        IppStatus ok = ippiHoughLineGetSize_8u_C1R(srcSize, delta, ipp_linesMax, &bufferSize);
        Ipp8u* buffer = ippsMalloc_8u_L(bufferSize);
        if (ok >= 0) {ok = CV_INSTRUMENT_FUN_IPP(ippiHoughLine_Region_8u32f_C1R, image, step, srcSize, (IppPointPolar*) &_lines[0], dstRoi, ipp_linesMax, &linesCount, delta, threshold, buffer);};
        ippsFree(buffer);
        if (ok >= 0)
        {
            lines.create(linesCount, 1, CV_32FC2);
            Mat(linesCount, 1, CV_32FC2, &_lines[0]).copyTo(lines);
            CV_IMPL_ADD(CV_IMPL_IPP);
            return;
        }
        setIppErrorStatus();
    }
#endif


Mat _accum = Mat::zeros( (numangle+2), (numrho+2), CV_32SC1 );
std::vector<int> _sort_buf;
AutoBuffer<float> _tabSin(numangle);
AutoBuffer<float> _tabCos(numangle);
int *accum = _accum.ptr<int>();
float *tabSin = _tabSin.data(), *tabCos = _tabCos.data();

// create sin and cos table
createTrigTable( numangle, min_theta, theta,
                 irho, tabSin, tabCos);

// stage 1. fill accumulator
for( i = 0; i < height; i++ )
    for( j = 0; j < width; j++ )
    {
        if( image[i * step + j] != 0 )
            for(int n = 0; n < numangle; n++ )
            {
                int r = cvRound( j * tabCos[n] + i * tabSin[n] );
                r += (numrho - 1) / 2;
                accum[(n+1) * (numrho+2) + r+1]++;
            }
    }

// stage 2. find local maximums
findLocalMaximums( numrho, numangle, threshold, accum, _sort_buf );

// stage 3. sort the detected lines by accumulator value
std::sort(_sort_buf.begin(), _sort_buf.end(), hough_cmp_gt(accum));

// stage 4. store the first min(total,linesMax) lines to the output buffer
linesMax = std::min(linesMax, (int)_sort_buf.size());
double scale = 1./(numrho+2);

lines.create(linesMax, 1, type);
Mat _lines = lines.getMat();
for( i = 0; i < linesMax; i++ )
{
    LinePolar line;
    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 = static_cast<float>(min_theta) + n * theta;
    if (type == CV_32FC2)
    {
        _lines.at<Vec2f>(i) = Vec2f(line.rho, line.angle);
    }
    else
    {
        CV_DbgAssert(type == CV_32FC3);
        _lines.at<Vec3f>(i) = Vec3f(line.rho, line.angle, (float)accum[idx]);
    }
}

}
stage1即核心逻辑,挨个遍历有效像素,统计出各种(theta,r)代表的直线穿过的像素点点的数量

Mat _accum = Mat::zeros( (numangle+2), (numrho+2), CV_32SC1 );
可以看到统计直线穿过的点数量的矩阵的个数是 (2 + numangle) x (numrho+2),即与我们传入的double rho, double theta有关.这个值越小,相应的我们搜索的直线数量就越多.

opencv的实现里有一些可能是出于工程上的考虑,这点不太确定,比如这里为什么要(2 + numangle) x (numrho+2) 而不是 numangle x numrho

int max_rho = width + height;
int min_rho = -max_rho;
为什么是w + h,而没有用开平方根求对角线长度.
希望知道的朋友可以留言告诉我.

// stage 2. find local maximums

static void
findLocalMaximums( int numrho, int numangle, int threshold,
                   const int *accum, std::vector<int>& sort_buf )
{
    for(int r = 0; r < numrho; r++ )
        for(int 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.push_back(base);
        }
}

寻找计数的局部最大值.类似于非极大值抑制.进一步细化检测到的直线,把局部的很相似的直线只取最精准的.

// stage 3. sort the detected lines by accumulator value
按accum数量大小排序

// stage 4. store the first min(total,linesMax) lines to the output buffer
保存前n条lines到输出Buffer.