霍夫变换 Hough Transform

by WZhang published 2026-02-07 views 67

Hough 变换原理与应用

前言： 详细介绍了 Hough 变换的基本思想、基本原理和应用等。

1. 基本概述

1.1 一些基本问题

Hough 变换简单概括就是原空间到参数空间的变换。

以一条直线为例，常用直线方程为 $y = kx+b$ ，这个是原空间。在这个直线方程里，x 和 y 是变量，k和b为直线方程参数。注意参数不是固定的，数量也不是固定的，比如这个直线不一定非要用 k,b （斜率和截距）这两个参数，下面会有举例说明。

何为到参数空间的变换呢？

直接了当地说就是，把参数当作变量（即坐标轴）。上面直线例子而言，就是把 k,b 当作横纵坐标构建的空间。注意，这一变换并非坐标系的变换（如：不是笛卡尔坐标向极坐标的变换；两者本质就不一样），Hough变换是始终以笛卡尔坐标系为基础，只是换了坐标轴（新空间中，坐标轴换成了原空间表达式中的参数，也是因此新空间叫参数空间。）

为何要变换到参数空间呢？

原空间中一些事物具有很强的关联性，但在原空间不好观察和把握，因此考虑从别的视角来表述它们关联性。参数空间就是一种，还有类似的傅里叶变换等等。所以，可以说，我们变换此空间，是为了更直观地揪出原坐标系下一系列点间的共同特点（或者说联系）。至于为什么选择参数空间，别问，问就是第一个想出这点的数学家的智慧和创新。并且，”很巧”的是，它在处理一些问题上确实很管用。

既然说参数选择不是固定的，如何选择最佳参数来构成参数空间？

这个不用操心，针对一些特定任务和问题，前人已经给出了答案了。

1.2 以例子说明

下面以三个例子说明上述自己理解得出的重点。

1.2.1 例子1：直线 $y = kx + b$ 到参数空间的变换（k，b为定值，如k=2，b=4）

事实上，使用这样例子并不好，容易给人带来误解，但又不得不从该例子说起。之所以说容易带来误解，原因在于：

我们遇到的问题，往往都是考察原坐标系下所给出的一系列点的关系问题。即基于点，去将这些点变换到参数空间；这个例子的表述，恰恰与此过程相反，该例子表述是已知结果（这些点是在一条直线上），去说明参数变换这回事。

两个是因果倒置的问题。这么说有些难以理解，下面例子会说明上面两句话要表达的含义。

我们考察直线到其参数空间的变换，这里参数就选择 k,b，变换如下图(a)(b)所示。

我们这样表述：左侧原空间下，对于直线 $y = kx + b$ ，将 $k,b$ 当作坐标轴的话，得到右侧参数空间，**该直线就对应于右侧空间一点 $(k,b)$ ** 。这样表述有些别扭，因为常规上，我们都是把问题落脚在点上，毕竟计算机存储和处理的都是离散点数据。

因此修改说： 左侧原空间下，直线 $y=kx+b$ 上的点对应右侧参数空间中一点$(k, b)$，如图(c)到(b)所示。 但这样说的话，同样别扭，因为这句话是错的。

Hough 变换不是常规的映射，不是点对点的对应关系，不是函数式的 f(x)。这点也可以说明它跟坐标系间变换有本质区别。

左侧原空间上的任意一点如 $(x_0,y_0)$，对应到右侧，不会是一个点，事实上是一条线。因为右侧表示的是直线的斜率和截距，而经过点 $(x_0,y_0)$ 的直线有无数条，也就对应着无数的 k 和 b。那是否意味着左侧这 $(x_0, y_0)$ 点，对应到右侧，就是充满坐标轴的所有点？当然不是，因为 k,b 虽然是随机的，取值可以 $[-\infty, \infty]$，但它俩不是各自随机，两者是有关系的。因为我们前面前提是经过 $(x_0,y_0)$ 点的直线，即

$$ y = k (x-x_0) + y_0 $$

说明：用这个 y = kx + b 这个式子不好，因为后面例子都是把 x、y 放到一侧，两者式子下，过点限制条件是由差别的。使用 $y-kx = b$或$y-kx-b=0$ 理解更好，这样的话，限制条件就是

$$ y - k x= y_0-k x_0\ \ \ \ 或 \ \ \ y-kx-b = y_0-kx_0 -b $$

也即

$$ y - kx = (y_0-kx) $$

这里会衍生出一个问题，即好像一个参数 k 就可以确定这一直线了。到这节最后会说到这个问题。先忽略之，或自己加以理解。

依据上式，也就是以 k， $y_0-kx_0$ 为参数，即 k 当横轴，$y_0-kx_0$ 当纵轴，或者写为我们更常见的

$$ b = -kx_0+y_0 $$

其中，b为纵轴，k为横轴。整个过程如下图(a)(b)所示

反过来，参数空间一点，就对应于原空间里斜率为 k，截距为 b 的直线。

总结一下，这两个空间的关系可以这样表述：左侧原空间的直线对应右侧参数空间一点；左侧原空间一点对应右侧空间一条直线。这是比较有意思的相互关系，圆形变换其实也有这个特点。参数即变量，变量即参数 的感觉。

基于此，我们就把刚才别扭的表述，说得准确些：对于左侧直线，我们用斜率 k 和截距 b 就可以描述它，并且唯一确定它，因此若建立一个斜率、截距空间，点(k，b)就可以表征一条原空间下斜率为 k 截距为 b 的直线。

那么是不是参数只能选择 k,b 呢？显然不是，比如你就可以选择参数空间里横坐标为 k+b, 纵坐标为 b，同等效果。但一般，参数选择要为计算、理解方便服务。

假设，我们发挥了自己的想象，采取了如下两个参数，即圆心到直线的距离 d 和这个垂线的角度 $\theta$，如下图所示。首先，我们明确的是，这两个参数是否可以唯一确定地表示这一直线？答案当时是肯定的，那么这个参数就能用。

那么，问题就来了，如果以 $d, \theta$ 为参数空间横纵坐标，那应该是怎样的。我们下面来考察一下：

首先，我们先把这条直线用 d 和 $\theta$ 表示出来，这个是高中数学题了，就不推导了（求解思路就是用 d 和$\theta$ 表征线上任意一点）。直接给出答案，这条线可以这么表示（这里没有常规地把 y 放在左侧，为了美观）

$$ d = xcos\theta + ysin\theta $$

同样的方法，对于无数 $d$ 和 $\theta$ ，我们要满足原空间里经过$(x_0,y_0)$，即

$$ x_0\cos\theta+y_0\sin\theta=x\cos\theta+y\sin\theta $$

也即有

$$ d = x_0\cos\theta + y_0\sin\theta $$

这里体现了 d 是有取值范围的。

根据三角函数的和差角公式，可以写成

$$ d = \sqrt {x_0+y_0}sin(\theta+\psi) $$

其中，$\psi$为常数，且$\psi=\arctan(\frac{x_0}{y_0})$。最终，我们就可以得到原空间直线上一点 $(x_0, y_0)$ ，在上述 $d$ 和 $\theta$ 参数空间下（以d为纵轴，$\theta$ 为横轴），所表示的东西。是什么呢？是一条正弦曲线。如下图 (b) 所示

也就是说，在原空间里的任意一点 $(x_0, y_0)$，如果变换到上述 $d,\theta$ 参数空间里，就是一条正弦曲线。同样是点对应线的关系。

上面是直接给出了使用 $d$ 和 $\theta$ 参数，一个是截距一个是角度，显得很突兀。但实际上，使用 d 和$\theta$ 参数，应用十分广泛，本质上，它是参考极了坐标的表达，另一方面，其应用广泛的原因在于 d 和 $\theta$ 都是有范围的，后面图像直线检测实例中会讲到。

上述只是直接给出了结果，是为了避免一些人将它与直角坐标系到极坐标系变换相混淆。下面给出推导：

使用极坐标系表示直角坐标里的点，有个对应关系，即

$$ x = \rho \cos \theta\ \ \ ,\ y= \rho \sin \theta $$

变换，即

$$ x\cos \theta = \rho \cos^2 \theta\ \ ,\ \ \ y\sin\theta = \rho\sin^2\theta $$

两等式累加，即

$$ x\cos\theta+y\sin\theta = \rho $$

也就是上面我们有提到的用 $d$ 和 $\theta$ 表示的直线公式。注意，这里的 $\rho$ 如果要是用于变换到参数空间，表示的不是直角坐标系里的一点到原点距离，而是过该点且垂线倾角为 $\theta$ 的直线距原点的距离，$\theta$ 同理。 $\rho,\theta$ 要是用于点变换到极坐标系下，就是我们常规表示，而且此时是点对点关系。再次强调两者不同点。

总结以上：原空间一点变换到参数空间的结果会随参数选择而变化，参数选择往往是基于计算方便、基于自己目的。

1.2.2 例子2：极坐标下的直线到参数空间的变换

问题来了，极坐标下直线一般方程是什么？同样是高中问题。

直角坐标下

$$ y = kx + b $$

使用常规推导，即$y=\rho\sin\theta,x=\rho\cos\theta$，带入

$$ \rho\sin\theta = k\rho\cos\theta + b $$

看着别扭，把 k 和 b 换下，即变成形如下式的形式

$$ \rho(\sin\theta-m\cos\theta)-n = 0 $$

其中，m,n 为定值。也即 m,n 知道后就可以确定极坐标系下唯一直线，因此可以用 m,n 构建该极坐标系下直线的参数空间，即对应于一点(m,n)。只是它的物理意义在坐标系（极坐标系）里不那么看得出来。m,n 的物理意义是，该直线（在极坐标系下）映射到直角坐标系下的斜率和截距。贴个极坐标系下直线图

上式也可以利用三角函数的和差公式变换下，即得

$$ \rho\sqrt{1+m^2} sin(\theta-\psi) -n=0 $$

也即

$$ \rho\sin(\theta-\psi)=\frac{n}{\sqrt{1+m^2}} $$

换下参数，即

$$ \rho\sin(\theta+\alpha) = b $$

其中，$\alpha,b$ 为常数，同样，$\alpha,b$ 知道后就可以确定极坐标系下唯一直线，同样可以用它构建参数空间。

三个小问题：

上面以$\alpha,b$ 为参数，变换到参数空间后，形状是什么？好奇的自行探究吧。
这里的$\alpha,b$ 可能不能随便取值，换句话说，并不是随意两个$\alpha,b$ 都可以表示一条直线。因为从公式演化来看，两者都是受 m,n 控制，自由度受控，可能有些 $\alpha,b$ 是无效的，但可以确定的是，极坐标内任意直线都可以找到一组 $\alpha,b$ 与之唯一对应。感觉是如此，不再去印证。
是否可以考虑参数 $\alpha，b$ 用极坐标系表示？即用极坐标系表示参数空间。我想应该是可以的。这是个有意思的问题。但常规上，都是在笛卡尔坐标系下。同样不再去深究。

1.2.3 例子3：圆到参数空间的变换

有了上面的经验，我们考察一下圆。前面已经说过，表达形式不同，可选参数也是多样的。这里只给出常用的参数选择。

圆的一般方程（只考虑笛卡尔坐标系下）：

$$ (x-a)^2+(x-b)^2 = r^2 $$

即圆心O坐标为 $(a,b)$，半径为 $r$。

选择 $a,b,r$ 为参数，注意，这里验证了前面提到的对于一些问题参数数量也是不一样。

按照前面对直线考察相同的步骤，参数空间里一点 (a,b,r) 就可以表示原空间里这一圆。不同之处在于，这里上升到了三维坐标。

同样考察对于原空间圆上一点 $(x_0, y_0)$ ，与直线相同的思路，经过该点有无数的圆，不同圆心O、半径r就对应着不同的过该点的圆。但同样(a,b,r)不是随意的，而是满足过$(x_0,y_0)$ 这一基本条件，即

$$ (x_0-a)^2+(y_0-b)^2=r^2 $$

这个也就是选择 a,b,r 作为坐标轴的参数空间下方程。写得顺眼点就是

$$ (a-x_0)^2 + (b-y_0)^2 = r^2 $$

更顺眼点，即

$$ (x-x_0)^2 + (y-y_0)^2 = z^2 $$

这个方程不就是圆锥方程吗。

参考直线时的表述，我们就可以说：过原空间一点的所有圆对应于参数空间里的一个圆锥。

如下图（a）(b)所示

【参考Fig.2自己想象吧，用python画太费劲了，懒。想象：左图(a)是一个点，然后有无数过该点的圆(用虚线表示)；右图(b)是一个圆锥。当然也可以有志之士帮忙画一下】

事实上，不如我们把两者综合起来说，即

对于原空间坐标中的一点，如果我们考察对象是过该点的直线，那么对应于参数空间里，就是一条直线；如果我们考察对象是过该点的圆，那么对应于参数空间里，就是一个圆锥

两个小问题：

维度相较于考察对象为直线时增加了一个，是不是很神奇？自行见解。
一个点对应了参数空间的一个圆锥，后面应用时(Hough圆检测)，计算量会非常非常大。因此一般不会直接应用，会有改进或别的思路。下面圆形检测部分会详解。

一个有意思的东西：

换个思路表示圆形，即

x = a + r$\cos\theta$

y = b + r$\sin\theta$

上面消去了$\theta$ ，会变成了圆一般方程 $(x-a)^2+(x-b)^2 = r^2$，这里我们消去r试试，即变成 $$ y-b = (x-a)\tan\theta $$ 如果以 $a,b,\theta$ 为参数会怎样呢？

一个更有意思的东西：

对于原空间中一点，考察经过该点的所有圆，事实上使用 (a,b) 就可以唯一表示一个圆了，即圆心为(a,b)，且经过该点，就已经确定唯一圆了，即只考虑参数 a,b 不就够了吗？

换个更简单的，对于原空间一点，考察经过该点的所有直线，事实上使用斜率 k 就可以唯一表示一条线了，即斜率为k，且经过该点，截距自然而然就确定了。那只考虑参数 k 不就够了吗？为什么参数空间里还要使用两个参数 k,b 来确定经过该点唯一直线。

其实这个问题有些诡辩的意味。稍微思考一下就可以反驳，如果仅用 k 表达，也即参数空间是个一维坐标轴，反推一下，参数空间里一个 k 能表达一条唯一的线吗？显然不能。加入b的原因，就可以简单理解为是为了把经过该点的信息囊括进去。

更简单直白点就是，参考例子 1 ，实际上 b 并不是单独参数，它只是 k的函数，即f(k)，如例子1中的纵轴其实是 $y_0-kx_0$。其目的就是为了表达它有个限制条件，即要经过$(x_0,y_0)$。

圆形Hough下，同样的原因。

至此，或许对 Hough 变换有了那么一点理解，或者说，有那么点印象了。

下面就直入主题，Hough变换的直线检测和圆检测，以及改进后的其它检测。

2. Hough 直线检测

将前面的结论再说一遍：

对于原空间里的一点，过该点有无数的直线，每条直线唯一对应着参数空间中的一点，而所有直线就对应参数空间里无数的点，这些无数点连起来就是一条直线$or$一条正弦线$or...$（取决于你选择的参数）。

先从解决一个简单问题入手。

2.1 问题

如何基于上述 Hough 变换的特点，检测三个点，即 $a(1, 2),b(3, 4),c(-1, 0)$ 是否在一条直线上？

2.2 思路

根据上面 Hough 直线变换的特点，这里使用 k,b 参数进行分析说明。

先考察 a(-1, 0)点，经过该点有无数的线，映射到 k,b 参数空间是一条直线。而 k,b 空间上该线一点也就对应于原空间这无数线中的其中一条。如下图(b)中红点，就对应于原空间红线。

下面给出参数空间曲线求解过程（可参考例子1 过程）：

经过 a(1,2) 点的所有直线可表示为： $$ y = k(x-1)+2 $$ 即 $$ y = kx-k+2 $$ 则参数空间中，横坐标为 k，表示原空间的直线斜率；纵坐标为 -k+2 ,表示原空间的直线截距，记为b。则参数空间中，线方程就是 $$ y = -x + 2 $$ 即 Fig.4 中图(b)所示。该线上每一点，对应左侧过点 a(1,2) 一条线。如红色所示。

同理方法，对于 b,c 点做一样的参数空间变换，最后结果如 Fig.5 所示

右侧参数空间里，有个特殊点 (1,1)，参数空间里三个线都经过这点。也即有

对于 a ，该特殊点含义是，原空间中，有个斜率为 k=1，截距 b=1 的直线经过它。
对于 b ，该特殊点含义是，原空间中，有个斜率为 k=1，截距 b=1 的直线经过它。
对于 c ，该特殊点含义是，原空间中，有个斜率为 k=1，截距 b=1 的直线经过它。

换句话说，原空间里有一条直线经过了这三点。即 Fig.5 (a) 中的红色线。

至此，我们达到了目的。即三点在一条线上。且该线的斜率为1，截距为1。也即图fig.5(b)中参数空间的红色交点。

2.3 总结

所以对于直线检测，我们只需对每个点进行参数空间的直线变换，然后查看参数空间的交点情况。

如果没有交点，就是不共线。如果有 N 个直线交于一点，对应原空间里就是 N 个点共线。我们可以把 N 称作叠加度，或者说参数空间中该点的亮度（名字乱起的）。亮度越高，共线的点越多，在原图中能直接观察到直线的可能性越大。

有意思的一点是，原空间两个点对应到参数空间会怎样？毕竟原空间中两点必共线，那参数空间是否必相交？参数空间里仅有两条线，且平行，对应原空间的两点是什么情况？

2.4 提升

2.4.1 实际问题

有了以上思路，就可以用 Hough 直线检测来解决实际问题了。

Hough 直线检测一般用于二值图中的直线检测，通常是一张图经过边缘算子后，检测边缘二值图是否存在直线。

如，利用 Hough 直线检测，检测 Lena 图中存在的直线。

2.4.2 求解思路

如果没有参考其他已有代码，而是基于上述一些知识，自己实现编程可能会遇到一些问题（比如我就是）。

即，如果我用 k,b 当参数空间，我可以很容易得到经过一点$(x_0,y_0)$，变换到参数空间里的直线为

$$ b = -x_0* k+ y_0 $$

但问题是，k 可以取无穷大，b 也是，毕竟是条直线。我怎么编程表示出这条直线，因为后续还要计算各点对应线的交点，叠加度。直线无法表示出，怎么求交点？

有一种思路是解方程，考虑所有点变换到参数空间后的线方程，即

$$ \left\{ \begin{matrix} \begin{aligned} b &= -x_0*k + y_0 \\ b &= -x_1*k + y_1\\ &\vdots \\ b &= -x_n*k + y_n \end{aligned} \end{matrix} \right. $$

问题转化为，得到一组组解（k,b），该解满足的方程数即为交点叠加度（或亮度）。具体解法，先将方程转化为参数矩阵，利用矩阵，然后慢慢折腾吧。（这是遇到问题后自己想到的，可行与否，未知，且不探究。感觉可行，实在不济，两两方程求解，把解带入其他方程验证。但计算复杂度可能会特别特别高。）

第二种思路就是使用之前提到的有取值范围的 $d,\theta$ 参数。 即原点到直线距离 d 和直线的法向量倾角 $\theta$ 。好处是，它俩都是有范围的。其中 $d\in[-\sqrt{x_0^2+y_0^2},\sqrt{x_0^2+y_0^2}]$ ,$\theta \in[0,2\pi]$ 。这似乎就可以用编程来实现了。

前面已经推导出了，过点 $(x_0,y_0)$ ，变换到参数空间 $d,\theta$ 后的曲线方程，即 $$ d = x_0\cos\theta + y_0\sin\theta $$ 或者 $$ d = \sqrt {x_0+y_0}sin(\theta+\psi)\ ,其中\psi=\arctan(\frac{x_0}{y_0}) $$ 编程实现如下：

def getLine(x0,y0,angs_resolution= 100):
    """
    x0: 原空间下点横坐标
    y0: 原空间下点纵坐标
    angs_resolution: \theta 划分精度
    功能：原空间(x0,y0)点，变换到 d, \theta 参数空间的曲线。
    说明：因为计算机存储是离散值，所以只是 \theta 取到一些值下的直线。当然，\theta 取值越多，越精细。
    """

    angs = np.linspace(0, 2*np.pi, angs_resolution) # 定义\theta 取到的离散值
    d = x0 * np.cos(angs) + y0 * np.sin(angs)
    return angs,d

测试一下，点 (x0, y0) 设为 (1, 1)，运行程序，可以得到如下结果

原空间中一点，对应到 $\theta,d$ 参数空间是一条正弦线。

求原空间所有点对应的正弦线，计算亮度

事实上，getLine()输出的是一系列离散点，点横坐标是 $\theta$ ，纵坐标是 d。我们查找所有输出值中 $(\theta,d)$ 重复次数就行了。重复次数，就是亮度。而且还有一点方便的是，对于每组getLine() 输出值，$\theta$ 都是相同的，因此我们只需检查每组对应位 d 是否相同（或小于某一阈值）即可。

2.4.3 自编程实现

完整代码：

import numpy as np
import matplotlib.pyplot as plt
import cv2

def getLine(x0,y0,angs_resolution= 100):
    """
    x0: 原空间下点横坐标
    y0: 原空间下点纵坐标
    angs_resolution: \theta 划分精度
    功能：原空间(x0,y0)点，变换到 d, \theta 参数空间的曲线。
    说明：因为计算机存储是离散值，所以只是 \theta 取到一些值下的直线。当然，\theta 取值越多，越精细。
    """
    
    angs = np.linspace(0, 2*np.pi, angs_resolution) # 定义\theta 取到的离散值
    d = x0 * np.cos(angs) + y0 * np.sin(angs)
    return angs,d

def Hough(edgeImg, angsDiv = 500, dDiv=1000):
    # 获取图像尺寸
    ySize, xSize = edgeImg.shape

    # 得到二值图中所有点坐标(x,y)
    y, x = np.where(edgeImg != 0)
    
    # 大致确定 d 范围
    dMax = np.sqrt(np.max(y)**2+np.max(x)**2)

    # 分辨率
    d_res = 2*dMax/(dDiv-1)
    ang_res = 2*np.pi/(angsDiv-1)

    # 亮度模板，起初为全黑，当经过某点，亮度 +1
    template = np.zeros((dDiv, angsDiv), dtype = np.uint8)
    
    # 亮度叠加计算
    for xx in range(len(x)):
        _, _d = getLine(x[xx], y[xx], angsDiv)
        _n = ((_d+dMax)/d_res).astype(np.int64)
        angle = np.arange(0, angsDiv)
        for p in zip(angle, _n):
            template[p[1], p[0]] = template[p[1], p[0]] + 1
    return template

if __name__ == "__main__":
    # 读取图片
    imgPath = "C:\\Users\\zhangwei156\\Desktop\\figure\\lena.bmp"
    grayImg = cv2.imread(imgPath, 0)
    
    # 提取边缘
    edgeImg = cv2.Canny(grayImg, 300, 500)
    
    # 设置离散的精度
    angsDiv = 500
    dDiv = 1000
    # 霍夫变换
    forceImg = Hough(edgeImg, angsDiv, dDiv)
    
    plt.imshow(edgeImg)
    plt.show()
    plt.imshow(forceImg)
    plt.show()

结果如下：

另！验证代码也贴上吧：

# 一些后面要用到的常数
y, x = np.where(edgeImg != 0)
dMax = np.sqrt(np.max(y)**2+np.max(x)**2)
d_res = 2*dMax/(dDiv-1)
ang_res = 2*np.pi/(angsDiv-1)

# 这是只考察了最大点，即亮度最大的点
ind = np.where(forceImg == np.max(forceImg))

# theta 和 d 的真实值
theta = (ind[1])* ang_res
d = -dMax + (ind[0]) * d_res

# 对应原空间的直线
xx = np.arange(512)
i = 0
yy = (d[i] - xx*np.cos(theta[i]))/np.sin(theta[i])

plt.plot(xx, yy, "r", linewidth = 0.3)
plt.imshow(edgeImg)
plt.show()

结果：

看起来，好像是那么回事！ 换个图片试试：

Nice!

说明：

为了编程方便，同时能形象展示亮度信息，我把角度 $\theta\in[0, 2\pi]$ 和距离 $d\in[-dMax,dMax]$ 范围各自等分了共 $angsDiv$ 份和 $dDiv$ 份。$forceImg$ 图像横纵轴坐标表示的是第几份(全是整数)，不是直接的 $\theta$ 和 d。如果求解真正的值，对于$d = -dMax + y * \frac{dMax - (-dMax)}{dDiv-1}$；$\theta$ 同理。
对于如上自编程的实现，分辨率 angsDiv, dDiv 设置真的很重要，设置过小的话。结果可能看起来好像不对，就是精度问题。况且里面还用到了许多四舍五入。
分辨率越大，计算效率越低。已经感受到了Hough变换耗时的地方。不过，一些库（如openCV）里应该会有精巧方法，对此有所优化。
一个可供改进的点：可以考虑将$\theta，d$ 相近的划归到一条直线，通过设定阈值实现，这样结果会更合理。
上面只使用了亮度最大的点，即只考虑原图中最直线的一条直线。
还有另一种编程思路，即借助另一种变换实现。即 Radon 变换。

3. Hough 圆形检测

3.1 传统圆形检测

跟直线检测同样的步骤。先来个简单的问题。

3.1.1 问题

如何基于上述 Hough 变换特点。检测四个点，即 $P_1(0,2), P_2(1,1), P_3(1,3), P_4(2,2)$ 在同一个圆上。

3.1.2 思路

前面1.2.3已经提及到，圆一般方程为

$$ (x-a)^2+(y-b)^2 = r^2 $$

选择 a,b,r 作为参数，变换到参数空间后，即对应一个圆锥，圆锥方程为

$$ (a-x_0)^2+(b-y_0)^2 = r^2 $$

即 过原空间中一点 $(x_0, y_0)$的所有圆对应于参数空间中一个圆锥 。

有意思之处：

空间中两点，对应参数空间两个圆锥，两圆锥相交是条曲线，即意味着：原空间里两点可以确定无数圆。符合实际。而且参数空间里这条相交的线有个特点，其在xoy平面投影是一条直线。对应到原空间的意思就是：经过这两点所有圆的圆心在一条线上。符合实际。

空间中三点，对应参数空间三个圆锥，三圆锥相交得到一点，该点的 $(a,b,r)$ 就对应着：原空间里有一个圆心为 $(a,b)$ ，半径为 r 的圆，经过这三点。如果空间中三点共线呢？对应到参数空间圆锥什么情况？自行考虑。

3.1.3 编程

虽然思路很清晰，即：编程表示出圆锥，然后求圆锥交线(点)，计算亮度。

但编程实践存在的问题是：

$r$ 虽然有点范围 $(r\ge 0)$，之所以说”有点“，在于 r 可以取 $+\infty$ 。另外 $a, b$ 可没边啊。
即便表示出了圆锥，那每次亮度计算都要遍历下三位空间像素点，计算量岂不是爆炸。

暂时撂下这些问题，针对现有问题，先实现之。针对上面两点：

本例问题中，只有四点，计算量尚可。
假如这四点是在一张 $W\times H$ 尺寸的图中，而且假设我们检测的要求是：整个圆要全部呈现在该图中。那么就可以知道 $a,b,r$ 的范围了。

首先，圆心确保在图中，就需要 $0 < a < W, 0 < b < H$

其次，半径有最大值，即斜对角距离，即 \(0

其实确定的这个范围，比预设的这个条件更宽泛些，比如该条件下，可以表示一部分非全部呈现在图中的圆。

思路确定后，编程实现如下：

import numpy as np
import cv2
import matplotlib.pyplot as plt

# 获得过原空间一点(x0,y0) 对应 参数空间中的圆锥
def get(x0, y0):
    mesh = np.meshgrid(a, b)
    a_x = mesh[0].flatten()
    b_y = mesh[1].flatten()
    r_z = np.sqrt((a_x - x0)**2 + (b_y - y0)**2)
    
    return r_z


p = np.array([[0,20],[10,10],[10,30],[20,20]])

# 假设有这个限制，即图高10，宽10，而圆“完整地”在图中。
W = 50
H = 50

# 精度
a_div = 100
b_div = 100
r_div = 100

# 离散取值
a = np.linspace(0, W, a_div)
b = np.linspace(0, H, b_div)
r = np.linspace(0, np.sqrt(W**2+H**2), r_div)
# 分辨率
a_res = W/(a_div-1)
b_res = H/(b_div-1)
r_res = np.sqrt(W**2+H**2)/(r_div-1)

# 三维亮度模板
template = np.zeros((a_div, b_div, r_div), dtype = np.uint8)

# 亮度叠加
for xx in p:
    _r = get(xx[0], xx[1])
    
    _x = np.arange(0, a_div)
    _y = np.arange(0, b_div)
    Mesh = np.meshgrid(_x, _y)
    x = Mesh[0].flatten()
    y = Mesh[1].flatten()
    z = (_r/r_res).astype(np.uint16)
    for _p in zip(x,y,z):
        template[_p[0], _p[1], _p[2]] = template[_p[0], _p[1], _p[2]] + 1


# 最终结果即三维亮度模板 template

3.1.4 结果

如果print(np.max(template)) ，理论上输出结果应该是 4。

因为直接观察也知道这四个点在一个圆上（圆心为 $(10,20)$，半径为 $10$），因此参数空间里四个圆锥应该相交到 $(10,20,10)$ 点上，也即模板template 在 (10, 20,10) 坐标处的亮度应该为 4。

但是，结果输出却是2。事实上，这是由于精度造成的。

尽管三维 template 不好展示，我们可以用 xoy 平面截取它。可以预想到的是，如果这么做的话，去如此截取四个圆锥，截面上会呈现四个圆。下面借助此思路，考察下精度（分辨率）的影响。

额外贴上检查代码

# 打印最大亮度 N
N = np.max(template)
print(N)

# 考察截取面情况(这是使用第10个层面)
plt.imshow(template[:,:,10])
plt.show()

# 最大亮度索引
ind = np.where(template == N)

# 真实值
a_ = ind[0]*a_res
b_ = ind[1]*b_res
r_ = ind[2]*r_res

# 打印最大亮度对应 a,b,r 的真实值
for p in (zip(a_,b_,r_)):
    print("a = %.2f,  b = %.2f,  r = %.2f"%(p[0],p[1],p[2]))

下面展示实验结果：

再贴一下，每张图对应的最大亮度：

图	分辨率	最大亮度 N	正确与否	(a,b,r)
(a)	50，50，100	2	×	一大堆，错的离谱
(b)	100，100，100	2	×	一大堆，错的离谱
(c)	1000，1000，100	3	×	一大堆，但近似地都在 (10,20,10) 周围
(d)	1000，1000，500	4	√	准确！（10，20，10）

这些已经很能说明问题了，当分辨率较小时，交点模模糊糊，而且”那么粗“的线，马赛克质量。当分辨率上去后，如最后一个图，交点就会很明确，且结果会准确。同时带来的问题是计算耗时变久。

BTW，上面图也展示了可以考察一部分不完全在图中的圆形。所以具体应用时，完全可以根据实际情况来手动规定 r 值，甚至圆心位置 (a,b) 范围。

3.1.5 说明

上面其实已经展示了如何中规中矩地应用前面提到的Hough变换原理，在一张图中进行圆检测。但是，仅仅检测 4 个点，程序计算时间就到了不可忍受的地步。

一个自己想到的解决思路时，适当调低分辨率，最终结果再使用聚类算法。如(c)中再使用一次聚类算法，就可以得到那一点。当然这种做法依旧起不到很大功效。毕竟计算复杂度在那摆着。

下面就是，用一些巧妙方法来提高计算效率。就称为进阶Hough圆形检测吧。

3.2 进阶圆形检测（梯度法）

该思路就要参考论文了。

论文：[1] Yuen H K , Princen J P , J Illingworth, et al. A Comparative study of Hough transform methods for circle finding[J]. Image and Vision Computing, 1990.

3.2.1 思路

先忘掉前面的 Hough 圆形检测，况且还涉及到坐标系、公式啥的，太复杂。

先考虑这么一种检测图像中圆形思路：先找圆心，再确定半径。

圆有这样一个性质，即：圆上一点与圆心的连线垂直于过该点的切线。 说人话就是，如下图所示

过切点 p ，且垂直切线的垂线4，经过圆心 o. 圆上任一点都是如此，如另一切点 q 等等

反过来说，如果一些点是在圆上，那么过该点垂直切线的垂线会相交于圆心。 如图中线3和线4 交于圆心o.

而且，最重要的一点是，不同于我们平时见到的点，一个常规的点是没有斜率信息的。但！图像上一点是可以求出该点“切线”的，即点的梯度信息（注意：这个是原图像的梯度信息，不是提取边缘后二值图的）。

这里的切线实际上是伪切线。因为它只反映了该点相对于周边点的变化趋势信息（梯度），或者术语点说就是只反映了局部信息。而且是特别小的局部。所以单看一点的话，偏差会非常大。非常大的机率，依据该梯度算得的垂线不经过实际圆心。

但是！虽然点梯度信息有时不是那么靠谱，可，如果图中有很多点是在同一圆上，数量多了，便总会冒出那么几个靠谱的点。

不靠谱点的垂线千奇百怪、全看心情，因此极小概率会发生三个或者更多数量的不靠谱垂线交于一个小小的点上。

但如果茫茫中存在三四个甚至更多点数量的靠谱点，它们相交于一点，那这圆心不就确定了嘛。

“相交于一点” 这个词如此熟悉，前面多次提到。编程的时候，自然而然想到用“亮度”表示。这里同样，亮度越大，该点是圆心的概率越高。

确定好圆心后，求解圆心到各个点的距离，距离一样的就是在一个圆上，这个距离就是半径。

即便如此，如果最终结果不那么理想，也情有可原，因为图像中点的梯度信息跟实际切线差距真的很大。（自我感觉，实际圆半径越大，差距越大）。它们是局部与整体的差距，是一滴水与西湖的差距，是一片落叶与秋天的差距。但也有一叶知秋，不是吗。

3.2.2 编程

下面就是激动人心的手撸代码阶段

关键点：（1）使用模板，获取亮度信息，亮度越高，是圆心几率越大。（2）梯度在灰度图上求，半径在二值图中确定。

求圆心

（1）获取图像中每点的梯度（主要是方向信息，即该点切线的垂线斜率）

该步可以使用 sobel 算子分别求取每个点的 x 方向和 y 方向上梯度$G_x,G_y$ ，求方向，合成下就好。该梯度方向就是垂线方向。至于 sobel 算子，比较简单基础，但用处很多，只简单说明

$3\times 3$ 大小的 sobel 核如下

$$ S_x = \left[ \begin{matrix} -1 &0 &+1\\ -2 &0 &+2\\ -1 &0 &+1 \end{matrix} \right]\ \ \ \ and\ \ \ \ \ S_y = \left[ \begin{matrix} -1 &-2 &-1\\ 0 &0 &0\\ +1 &+2 &+1 \end{matrix} \right] $$

从该核分布也能看出其思想：求解梯度时，是考虑到了周边 8 点对该点影响。而且进行了加权，求取水平方向时，$S_x$ 紧邻该点左右像素加权值最大（即影响最大，符合实际情况），该点竖直方向零加权（这也突出了求解水平方向梯度，跟竖直方向没有半毛钱关系）。$S_y$ 同理。

求解是通过卷积实现，即：

$$ G_x = S_x * A \ \ \ \ \ and \ \ \ \ \ G_y = S_x * A $$

其中，A 为图像。当然也可以考虑旋转下核，如旋转 45度，135度等，就可以获得斜向梯度。Sobel算子可以检测边缘。Canny边缘检测就是使用的Sobel算子。 Sobel算子很重要，它相当于数学中的求导，应用也十分广泛。挖坑：sobel 算子基本原理与实现。

回到问题，该步编程实现如下：（使用自带的sobel函数）

def getK(img):
    # 求取 x,y 方向梯度
    Gx=cv2.Sobel(img, cv2.CV_64F, 1, 0)
    Gy=cv2.Sobel(img, cv2.CV_64F, 0, 1)
    
    # 求取垂线斜率，为了分母不为0，用了一个很笨的操作
    k = Gy/(Gx+0.0000000001)
    return k

说明：这里 Sobel 算子计算后数值设置为了 cv2.CV_64F 这个很关键，因为 $G_x, G_y$ 可以为正或为负数。目前这么设置后面不用考虑因设置了如默认的-1 而发生的截断。

（2）求解垂线：即过边界点$(x_0, y_0)$，斜率为 k 的直线

该直线方程就是简单的：

$$ y=k(x-x_0)+y_o $$

（3）求解所有边界点对应的垂线的交点

跟前面用到多次的技巧一样，使用模板，亮度叠加，最亮的地方即相交次数最多的地方。

该阶段代码如下：

import numpy as np
import cv2
import matplotlib.pyplot as plt

def getK(img):
    # 求取 x,y 方向梯度
    Gx=cv2.Sobel(img, cv2.CV_64F, 1, 0)
    Gy=cv2.Sobel(img, cv2.CV_64F, 0, 1)
    
    # 求取垂线斜率，为了分母不为0，用了一个很笨的操作
    k = Gy/(Gx+0.0000000001)
    return k

# 读取图像
imgPath = "C:\\Users\\zhangwei156\\Desktop\\figure\\logo.png"
grayImg = cv2.imread(imgPath, 0)
grayImg=cv2.GaussianBlur(grayImg,(5,5),0)

# 提取边缘
edgeImg = cv2.Canny(grayImg, 300, 500)

H, W = grayImg.shape

# 确定灰度图中每一点“切线”垂线的斜率
k = getK(grayImg)

# 边界点坐标
y, x = np.where(edgeImg != 0)

# 经过边界点(x0,y0)，斜率为 k 的垂线方程 y = k(x-x0) + y0，后面就是求不同线交点，即求亮度叠加
# 已是常规操作，前面已经重复多遍
# 精度
x_div = W*5
y_div = H*5
# 离散取值
x_ = np.linspace(0, W, x_div)
# 分辨率
x_res = W/(x_div-1)
y_res = H/(y_div-1)
# 亮度模板
template = np.zeros((y_div, x_div), dtype = np.uint8)
# 亮度叠加
for p in zip(x, y):
    # 垂线方程
    y_ = k[p[1],p[0]] * (x_ - p[0]) + p[1]
    # 垂线经过的点
    x_n = (x_/x_res).astype(np.uint64)
    y_n = (y_/x_res).astype(np.uint64)
    # 删去y超出范围的点，然后亮度叠加
    ind = np.where(y_n >= y_div)
    y_n = np.delete(y_n,ind,None)
    x_n = np.delete(x_n,ind,None)
    template[y_n, x_n] = template[y_n, x_n] + 1

# 亮度图 template
forceImg = template
plt.imshow(forceImg)
plt.show()

# 圆心坐标为 O
O = np.where(forceImg == np.max(forceImg))
# 真实值
O_y = O[0]*y_res
O_x = O[1]*x_res

for p in (zip(O_x,O_y)):
    print("O_x = %d,  O_y = %d"%(p[0],p[1]))

结果

至此，得到圆心坐标O（即亮度最大的地方）。

求解半径

使用了直接遍历的方法，如下（这里可以直接略过，因为下面方法太不准确了）

# 规定半径最大、最小值和分辨率
rMin = 0
rMax = max(H,W)*2
rRes = 5

# 限定一个圆上，一段弧长需要至少多少点
circlePointsNum_min = 7
# 储存符合条件的半径
rList = []
for _o in (zip(O_x,O_y)):
    # 求边界各点到圆心o距离
    R = np.sqrt((x - _o[0])**2 + (y - _o[1])**2)
    # 假设 r 进行一个一个地尝试
    for r_ in np.arange(rMin, rMax, rRes):
        if r_ == 0:
            continue
        indTemp = np.where(np.abs(R - r_) < rRes/2)
        # 当半径为r时，圆上边界点数目符合要求，就储存该半径值
        if indTemp[0].shape[0]/r_ > circlePointsNum_min:
            rList.append(r_)

# 显示所有圆形
plt.imshow(grayImg)
for _o in (zip(O_x,O_y)):
    angle = np.linspace(0,2*np.pi,50)
    for _r in rList:
        xxx = _o[0] + _r*np.cos(angle)
        yyy = _o[1] + _r*np.sin(angle)
        plt.plot(xxx,yyy, color="r")
        plt.scatter(_o[0],_o[1])
plt.show()

结果

通过修改 rRes 和 circlePointNum_min 值，可以得到还行的结果。

验证一个猜想

前面有提到，实际圆半径越大，可能圆心准度越小。猜想（a）图总调不到理想结果的原因，除了半径限制有问题，另一个圆心不准确可能就与图尺寸较大相关（或者说圆形太大）。为了验证，把(a)图 grayImg = cv2.resize(grayImg, (int(H/5),int(W/5))) 尺寸缩小一下，看结果。

果然！确实如此。此种方法，此种编程下，实际圆尺寸越小，得到的结果越准确。所以有时缩小一下图像是好事，而且减少计算量。

3.2.3 结果说明

上图例（a）中难以得到精确结果的原因是：圆心不精确

上述代码存在问题一便是：直接且只取了圆心亮度最大值。理论上还存在多圆多圆心。
存在一些某名奇妙的圆的原因是：半径限制条件有问题

上述半径阶段代码不够智能。第二阶段问题解决地有很大瑕疵和改进空间。

折腾地头疼，不再研究。
圆形检测时，有时结果错的离谱并不奇怪，相反，应当是相当正常的现象。局部反应整体，有时出现错误是可以理解的。
整个思想是先确定圆心：滤波、梯度计算的改进（如参考更多点信息，扩大局部求解范围）可以提高该部分精度。后确定半径：考虑圆上点的特点或圆弧特点来智能地确定半径。

3.2.4 写在后面

前面折腾代码的目的在于理解其原理。难做实际应用。实际应用参考openCV。

前面提到忽略之前提到的 Hough变换，现已经明白了梯度法原理。再回头看，所以，梯度法体现的 Hough变换思想在哪？自行见解。

3.3 openCV对应代码解读

重在理解原理。需要再解读之。

3.4 Hough 变换椭圆检测

考虑旋转、平移，椭圆的一般方程为

$$ x^2 +Ay^2+ Bxy + Cx + Dy + E = 0 $$

也就是说，需要五个参数。换句话说，霍夫参数空间上升到了五个维度。基于前面圆检测的计算量，即便针对椭圆有改进算法，估计也快不到哪里去。实际应用很受限，理解原理也就差不多了。基本原理同直线和圆大致相似。

再换句话说，这是礼貌且仪式性的一节。

4. 广义霍夫变换

将 Hough 变换检测推广到了任意形状。基本思想是模板匹配。

同样，计算量很大，应用很受局限，但是思想很有趣，实现过程很有趣。有空闲时，会另起一篇补充。

0comment(s)