霍夫变换不仅可以找出图片中的直线,也可以找出圆,椭圆,三角形等等,只要你能定义出直线方程,圆形的方程等等.
不得不说,现在网上的各种博客质量真的不行,网上一堆文章,乱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源码解读
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.