浅析霍夫变换检测直线和圆

原文：http://blog.csdn.net/shanchuan2012/article/details/74010561

1.本文目的

结合实例和一些演示来理解霍夫变换中的精髓：空间转换。

我想我把我对空间转换的理解写完了，这篇文章也就结束了，所以文章的内容不会那么完整，不会完整的讲霍夫变换（比如一般形状的检测，说实话我是初学也不会），也不会讲在opencv中的实现（我只能抄代码，自己也写不出来）。

为什么要写这篇文章？我是自学+初学图像处理，没有人可以问，只能上网搜索，网上的教程总有地方我看不明白，有时候我会怀疑是不是自己理解力不行？其实更多时候我觉得不是，问题其实在于两点，一是教程讲的不好，不能引领新手逐步深入，二是没有找到合适的教程，不同阶段看的东西是不一样的。所以我把自己的理解写出来，希望可以帮助跟我有同样困惑的人，也希望自己不要误导了别人，让人产生更多的困惑。

2.简介

这篇文章不是教程，不讲基本概念，要了解霍夫变换可以先参考以下文章：

【OpenCV入门教程之十四】OpenCV霍夫变换：霍夫线变换，霍夫圆变换合辑-毛星云

以一般化视角串联霍夫变换(hough transform)，从直线到圆再到广义霍夫变换

Hough transform(霍夫变换)

看了以上资料我有很多地方是想不明白的，特别是看到检测圆的时候，更别说一般形状了。

刚看到霍夫变换的时候感觉挺神奇的，它来源于图像处理中的几何形状检测问题，用什么办法能够检测直线、圆或者别的形状呢？最简单的是直线，如果让我来解决，要用什么办法才能把直线检测出来呢？一时也想不到什么办法（最小二乘可以做拟合，但是一张图上有多条直线怎么办？）。看到霍夫变换的时候想到了傅里叶变换，在这个空间不好解决的问题或者不明显的特征，给他来个变换，在另外一个空间上看问题，说不定就清楚了。

3.霍夫变换检测直线

3.1笛卡尔坐标下的例子

假如我们有如下几个点：

它们完全在一条直线上，可以先告诉你们直线的方程是：$y=x+4$。在我们不知情的情况下要如何检测？

直线的一般表示是：$y=kx+b$，$k$和$b$是两个参数，现在反过来，将两个参数当做自变量和变量，将方程写成：$b=-xk+y$，这个时候坐标轴就发生了变化，也就是空间发生了转换，原来坐标轴是$xy$，现在是$kb$。

以点$(1,5)$ （该点所在空间此处称为图像空间，因为我们是在做图像处理，这些点都对应于图像的像素）为例，代入直线方程$y=kx+b$ ，确定一条直线至少要两个点，只有一个点是不能确定一条直线的，但可以知道的是，直线一定过这个点，过一个点的直线有无数条，像这样（这个时候已经可以思考了，每个点上都有无数条直线经过，那么怎么找出同时经过这些点的直线？特征是什么？没错，用条件：“斜率和截距一样”就可以把那条直线找出来，这也是做空间变换的原因）：

由此得到：$b=-k+5$ ，这是一条直线的方程：

所以在$xy$空间中的一个点，对应于$kb$空间（此处称为参数空间）上的一条线。这个比较好理解，所谓的$kb$空间，代表了直线的斜率和截距，既然过一个点的线有无数条，那么斜率和截距也有无数个，$kb$空间上的线就描述了斜率和截距的无数种组合，所以形成了一条线。将图像空间上所有对应的点在参数空间上的直线都画出来就是这样：

可以看到它们相交于一个点了，这个点是$(k=1,b=4)$ ，没错，它就是我们要找的特征点（前面说过要怎么从无数条直线中找到经过这些点的直线），这正是图像空间上所有点连成的直线的斜率和截距。这样就把图像空间上的那条直线检测出来了，方程为$y=1\ast x+4$ 。

但是实际上用极坐标，为什么要用极坐标而不用笛卡尔坐标？上面给了回答，这里不解释。

3.2极坐标下的例子

如何用一个角度和一段距离来表示过一个点的直线呢？方法是从原点作该直线的垂线，原点到该直线的距离用$\rho$表示，垂线与x轴的夹角用$\theta$表示，这样就可以用$(\theta ,\rho )$来表示过一个点的任意直线了。

还是在笛卡尔坐标下，$(\theta ,\rho )$的关系用$(x,y)$来表示，可以写为（如果不明白，可以先看看参考资料，这里面的示意图讲的很清楚）：
$$\rho =x\ast cos\theta +y\ast sin\theta$$
还是一样的思想，对于图像空间上的一个点，过这个点的直线有无数条，所以$(\theta ,\rho )$也有无数个组合，那么在参数空间就对应于一条线（这时是一条曲线，管他是什么线呢，反正它们都是描述那无数条直线的参数组合）。还是以点$(1,5)$ 为例，它对应于参数空间的线：

将图像空间上所有对应的点在参数空间上的直线都画出来就是这样：

它们也相交于一点，把交点找到（我这里没有计算交点，就当它是$({\theta }_0,{\rho }_0)$吧，大概是$(2.3562,2.7266)$ , 其实对应的是$(135^{\circ },4\ast sin(135/180))$ , 这是我反推的）后，代入方程：$\rho =x\ast cos\theta +y\ast sin\theta$ 就可以得到过所有点的直线了：
$$2.7266=x\ast cos(2.3562)+y\ast sin(2.3562)$$
整理后就是：$y=x+4$ .到此直线检测结束。

4.霍夫变换检测圆

都说检测圆和检测直线大体上是类似的，确实，从空间转换的角度来说是一样的，椭圆、任意形状的检测都是一样的，只要找到合适的转换方法。

圆的一般方程是这样的：
$$(x-a{)}^2+(y-b{)}^2=r^2$$
$(a,b)$为圆心，$r$为半径，还是先给一个例子：

图上面有一系列的点，它们是在同一个圆上，这里已经把圆的轮廓给出来了，这里只画了圆的一半，但已经足够说明问题了。要检测这个图像中的圆应该怎么办？还是用同样的思考方式，来个空间转换。将上面圆的方程重新写一下：
$$(a-x{)}^2+(b-y{)}^2=r^2$$
这时有三个参数，圆心$(a,b)$，半径$r$。

过一点的圆有多少个？无数个。以上图的第1个点$(-2,1)$为例，把它代入方程$(a-x{)}^2+(b-y{)}^2=r^2$，得到：$(a+2{)}^2+(b-1{)}^2=r^2$，这个方程描述的圆都经过点$(-2,1)$，它们在参数空间上看起来是这样的：

这里$r$取的是$[1,2,3,4,5]$，没取成无限的，要是取成无限了，图就没有办法看了。每一个参数空间上的圈圈上的点都对应于图像空间上经过点$(-2,1)$的一个圆，如下面这个点：

它代表了什么意思呢？它代表了一个经过点$(-2,1)$的一个圆，圆心为$(-2,0)$，半径为1（图上的半径），这样是不是就清楚了一些？从三维空间来看就是这样（为了能看得清楚，只画了半径为1，3，5的三个圆）：

上面的图是图像空间中一个点对应在参数空间的曲线，将所有点对应在参数空间的曲线都画出来就是这样的（同样为了看的清楚，只取了3个半径，分别为1，3，5）：

其实我们的目的就是在参数空间找聚集点，思想和检测直线是一样的，找参数空间中曲线交点最多的位置，可以看到在$r=3$的时候所有圆都相交于同一个点，这个点对应的圆就是图像空间中经过所有点的那个圆，这样就把圆给检测出来了。还可以看出，当r离真实值（$r=3$）越远时，圆的相交的就越不明显。

我想这时候对霍夫变换检测直线和圆的原理应该说清楚了，从直线到圆，参数从2个变成了3个，对应的参数空间从2维变成了3维。如果是别的形状，那么它的参数可能更多，对应的参数空间维数就越多，这时就很难画图了。不过通过上面的理解，相信我们已经不需要再借助图像来理解更高维的变换了，这时我们应该可以进行归纳了，不管多高维，我们都是去找在参数空间中曲线相交最多的位置。

5.数据不是那么完美的情况

上面给出来的直线和圆都完全在同一条直线和同一个圆上，但实际上可能没有那精确像下面这样：

那么参数空间中对应的相交点也不是完全在一个点上

圆的情况也是类似的，就不再给出。所以一般的做法是将参数空间分成很多个小格子，落在同一个格子里面的就算同一个交点。

6.小结

我对霍夫变换的理解依然很肤浅，觉得奇妙的地方就是空间转换。

这里再来看以一般化视角串联霍夫变换(hough transform)，从直线到圆再到广义霍夫变换对作者简练的总结更能体会吧。

当然在实际的算法实现中，远比上面的图示来的复杂，不然后效率是个问题。

7.作图代码

这里附上作图用的matlab代码：

% 霍夫变换理解 圆 2维
close('all');

% 给出一个圆，圆心（1，1），半径3
r = 3;
x = -2 : 0.5 : 4; % x定义域
y = sqrt(r.^2 - (x-1).^2) + 1 + 0*(rand(1,length(x)) - 0.5); % 半个圆

% 图像空间上的圆
figure;
plot(x,y,'*'); % 实际点
hold on;
x0 = -2 : 0.01 : 4;
y0 = sqrt(r.^2 - (x0-1).^2) + 1;
plot(x0,y0); % 圆曲线
hold off;
axis equal; % 各个轴比例相同
xlabel('x');
ylabel('y');

% 参数空间
figure;
% set(gcf,'position',[100,100,500,500]);
theta = 0 : 0.01 : 2*pi;
for ir = 1 : 5
    for i = 1 : 1%length(x)
        a = x(i) - ir*cos(theta); % 极坐标
        b = y(i) - ir*sin(theta);
        plot(a,b); % 画圆
        hold on;
        plot(x(i),y(i), '*'); % 画圆心
    end
end

hold off;
grid on;
set(gca, 'xtick', -6:1:8);
axis equal;
xlabel('a');
ylabel('b');
print(gcf, '-djpeg', 'c2_2.jpg', '-r72');

% 霍夫变换(Hough Transform)
% 直线
clc;clear;close('all');

% 测试点
x = [0,1,2,3,4];
y = x + 4 + 1*(rand(1,length(x))-0.5);
plot(x,y,'*');
hold on;
y1 = x + 4;
plot(x,y1);

xlabel('x');
ylabel('y');
print(gcf, '-djpeg', 'line1.jpg', '-r72');

% 笛卡尔坐标
figure;
legd = {};
k = 0:3;
for i = 1 : length(x)
    b = -x(i)*k + y(i);
    plot(k,b)
    hold on;
    legd = [legd, num2str(i)];
end
hold off;
% legend(legd);
grid on;
xlabel('k');
ylabel('b');
print(gcf, '-djpeg', 'line2.jpg', '-r72');

% 极坐标
figure;
theta = 0/180*pi : 0.1 : 180/180*pi;
for i = 1 : length(x)
    ruo = x(i)*cos(theta) + y(i)*sin(theta);
    plot(theta, ruo);
    hold on;
end
hold off;
grid on;
xlabel('\theta');
ylabel('\rho');
print(gcf, '-djpeg', 'line3.jpg', '-r72');

% 霍夫变换理解 圆 3维
close('all');

% 给出一个圆，圆心（1，1），半径3
r = 3;
x = -2 : 0.5 : 4; % x定义域
y = sqrt(r.^2 - (x-1).^2) + 1 + 0*(rand(1,length(x)) - 0.5); % 半个圆

% 图像空间上的圆
figure;
plot(x,y,'*'); % 实际点
hold on;
x0 = -2 : 0.01 : 4;
y0 = sqrt(r.^2 - (x0-1).^2) + 1;
plot(x0,y0); % 圆曲线
hold off;
axis equal; % 各个轴比例相同
xlabel('x');
ylabel('y');
print(gcf, '-djpeg', 'c1.jpg', '-r72');

% 参数空间
% 3维图画的很慢，可能有地方弄错了
figure;
% set(gcf,'position',[100,100,500,500]);
theta = 0 : 0.01 : 2*pi;
for ir = 1 : 2 : 5
    for i = 1 : 1%length(x)
        a = x(i) - ir*cos(theta); % 极坐标
        b = y(i) - ir*sin(theta);
        r = ir*ones(1,length(a));
        plot3(r,a,b,'b'); % 画圆
        hold on;
        plot3(r,x(i),y(i), 'r*'); % 画圆心
    end
end

hold off;
grid on;
set(gca, 'xtick', -6:1:8);
% axis equal;
xlabel('r');
ylabel('a');
zlabel('b');
print(gcf, '-djpeg', 'c2.jpg', '-r72');