opencv之霍夫曼变换

  • 2019 年 10 月 3 日
  • 笔记

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

不得不说,现在网上的各种博客质量真的不行,网上一堆文章,乱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.