数字图像处理(4版)——第 10 章——图像分割(上)(Rafael C.Gonzalez&Richard E. Woods)

上一章的内容开启了图像处理方法的转变，从输入和输出均为图像的方法，过渡到输入为图像但输出为从这些图像中提取的属性的方法。本章中的大多数分割算法都基于图像强度值的两个基本属性之一：不连续性(discontinuity)和相似性(similarity)。第一类方法是基于强度值的突变(例如边缘(edges))将图像分割成区域。第二类方法则基于一组预定义的标准将图像分割成相似的区域。阈值分割、区域增长以及区域分割和合并都是此类方法的示例。在本章中，我们展示了如何通过结合不同类别的方法(例如将边缘检测与阈值分割相结合)来提高分割性能。我们还讨论了使用聚类和超像素进行图像分割的方法，并介绍了图割，这是一种非常适合提取图像主要区域的方法。接下来讨论了基于形态学的图像分割方法，该方法融合了本章第一部分介绍的几种分割技术的特性。本章最后简要讨论了利用运动线索(motion cues)进行分割的方法。

10.1. 基础知识(Fundamentals)

令 R 表示由一幅图片占据的整个空间区域。我们可以将图像分割视为将 R 分割成满足如下条件的 n 个子区域 $R_{1} ,R_{2} ,{\hdots} , R_{n}$ 的过程：

(a) $\bigcup_{i=1}^{n}R_{i} = R$ 。

(b) $R_{i}$ 是一个连通区域( 对 i = 1, 2 ,… , n ) 。

(c) $R_{i}{\bigcup}R_{j} = \emptyset$ (对所有 i 和 j 且 i ≠ j )。

(d) $Q(R_{i}) = \text{TRUE}$ ( 对 i = 0 , 1, 2 ,… , n ) 。

(e) $Q( R_{i}{\bigcup}R_{j} ) = \text{TRUE}$ ( 对任意邻接区域 $R_{i}$ 和 $R_{j}$ ）。

其中， $Q(R_{k})$ 是定义于集合 $R_{k}$ 中的点上的一个逻辑谓词(predicate)(译注：谓词，即下判断的词，得出真假结论)，而符号 ∅ 表示空集。符号 ∪ 和 ∩ 分别表示集合并运算和集合交运算。同 2.6 节中的定义。对于区域 $R_{i}$ 和 $R_{j}$ ，若其并集构成一个连通集，如 2.5 节中所定义，则称二者是邻接的。若由两个区域之并集构成的集合不是连通的，则称其为脱离的(disjoint)。

条件 (a) 表明，分割必须是完备的，即每一个像素必须位于一个分割区域。条件 (b) 要求一个区域中的点在某个预定义的意义上是连通的(例如，点必须是8 连通的)。条件 (c) 指的是区域必须是脱离的(不相交的)。条件 (d) 处理分割区域中的像素必须满足的属性—— 例如，若 $R_{i}$ 中的所有像素都具有相同的强度值，则 $Q(R_{i}) = \text{TRUE }$ 。最后，条件 (e) 表明，两个邻接区域 $R_{i}$ 和 $R_{j}$ 在谓词 Q 的意义上必须不同。( 注：通常，Q可以是一个复合表达式，例如，“若区域 $R_{i}$ 中的像素平均强度值小于 $m_{i}$ AND 若其强度的标准偏差大于 $\sigma_{i}$ ， $Q(R_{i}) = \text{TRUE}$ ，” 其中， $m_{i}$ 和 $\sigma_{i}$ 是具体的常量。 )

因此，我们看到图像分割的根本问题在于将图像分割成满足上述条件的区域。单色图像分割算法通常基于两种基本类别之一，这两种类别都与强度值的特性有关：不连续性和相似性。在第一类方法中，我们假设区域边界彼此之间以及与背景之间有足够的差异，从而可以基于强度值的局部不连续性进行边界检测。基于边缘(edge)的分割是此类方法的主要方法。第二类方法中的基于区域的分割方法则基于将图像分割成根据一组预定义标准相似的区域。

图 10.1 展示了前面的概念。图 10.1(a) 显示了一幅图像，其中包含一个强度恒定的区域，该区域叠加在一个强度也恒定的较暗背景上。这两个区域共同构成了整幅图像。图 10.1(b) 显示了基于强度不连续性计算内部区域边界的结果。边界内外的点均为黑色(0值)，因为这些区域不存在强度不连续性。为了分割图像，我们赋予边界上或边界内的像素一个级别(例如，白色)，赋予边界外的所有点另一个级别(例如，黑色)。图 10.1(c) 显示了这种分割过程的结果。我们可以看到，本节开头所述的条件 (a) 到 (c) 均满足此结果。条件 (d) 的谓词为：如果一个像素位于边界上或边界内，则将其标记为白色；否则，将其标记为黑色。我们看到，对于图 10.1(c) 中标记为黑色或白色的点，该谓词为真。类似地，两个分割区域(对象和背景)满足条件 (e)。

接下来的三幅图展示了基于区域的分割方法。图 10.1(d) 与图 10.1(a) 类似，但内部区域的强度构成了一种纹理图案。图 10.1(e) 显示了计算该图像中强度不连续性的结果。由于强度存在大量虚假(spurious)变化，因此难以识别原图像的唯一边界，因为许多非零强度变化都与边界相连，因此基于边缘的分割方法并不适用。然而，我们注意到外部区域是恒定的，因此解决此分割问题只需要一个能够区分纹理区域和恒定区域的谓词。像素值的标准差正是实现这一目标的度量，因为它在纹理区域内非零，而在其他地方为零。图 10.1(f) 展示了将原图像分割成 8 × 8 大小的子区域的结果。每一个子区域如果其像素的标准差为正(即谓词为 TRUE )，则标记为白色；否则标记为零。由于 8 × 8 的正方形组被标记为相同的强度(较小的正方形会产生更平滑的区域边界)，因此结果在区域边缘呈现出“块状”外观。最后，请注意，这些结果也满足本节开头所述的五个分割条件。

-----------------------------图 10.1：(a) 恒定强度区域的图像。(b)基于强度不连续性的边界。(c) 分割结果。(d) 纹理区域的图像。(e) 强度不连续性计算结果(注意大量细小边缘)。( f ) 基于区域属性的分割结果。--------------------------------------

10.2. 点，线，和边缘检测(Point , Line, and Edge Detection )

本节重点介绍基于检测强度局部突变的分割方法。我们感兴趣的三种图像特征是孤点、线和边缘。边缘像素是指图像强度发生突变的像素，而边缘(或边缘线段)则是一组相连的边缘像素(关于连通性，请参见第 2.5 节)。边缘检测子(edge detectors)是用于检测边缘像素的局部图像处理工具。一条线可以看作是一条(通常)较细的边缘线段，其中线两侧的背景强度远高于或远低于线像素的强度。事实上，正如我们稍后将讨论的，线会产生所谓的“屋脊边缘(roof edges)”。最后，孤点可以看作是被背景(前景)像素包围的前景(背景)像素。

10.2.1 背景(Background)

正如我们在3.5节中看到的，局部平均可以平滑图像。鉴于平均类似于积分，因此很直观地，可以使用导数来检测强度的突变。稍后我们将阐明，一阶和二阶导数尤其适用于此目的。

数字函数的导数是根据有限差分定义的(译注：“差分”其实就是差值，古汉语中指“指比例分配”)。计算这些差分的方法有很多种，但正如3.6节所述，我们要求用于一阶导数的任何近似值：(1) 在恒定强度区域必须为零；(2) 在强度阶跃或斜坡的起始点必须非零；(3) 在强度斜坡上的点必须非零。类似地，我们要求用于二阶导数的近似值：(1) 在恒定强度区域必须为零；(2) 在强度阶跃或斜坡的起始点和结束点必须非零；(3) 在强度斜坡上必须为零。由于我们处理的是数值有限的数字量，因此最大可能的强度变化也是有限的，并且变化发生的最短距离是相邻像素之间。

我们通过将一维函数 f ( x + Δx ) 展开成关于 x 的 Taylor 级数，来获得其在任意点 x 处的一阶导数的近似值:

(10-1)

$\begin{array}{lll} f ( x + {\Delta}x ) &=&\displaystyle f ( x ) + {\Delta}x\frac{{\partial}{f(x)}}{{\partial}x}+\frac{({\Delta}x)^{2}}{2!}\frac{{\partial^{2}}{f(x)}}{{\partial}x^{2}}+\frac{({\Delta}x)^{3}}{3!}\frac{{\partial^{3}}{f(x)}}{{\partial}x^{3}}+\hdots \\ \\ &=&\displaystyle \sum_{n=0}^{\infty}\frac{({\Delta}x)^{n}}{n!}\frac{{\partial^{n}}{f(x)}}{{\partial}x^{n}} \end{array}$

其中，Δx 是 f 的样本之间的间隔。对于我们的目的，这个间隔我们用像素单位度量。因此，在本书中遵循这个约定，对于 x 之前的样本，Δx = 1 ，对于 x 之后的样本，Δx = -1 。当 Δx = 1 时，公式 (10-1) 成为

(10-2)

$\begin{array}{lll} f ( x +1 ) &=&\displaystyle f ( x ) + {\Delta}x\frac{{\partial}{f(x)}}{{\partial}x}+\frac{({\Delta}x)^{2}}{2!}\frac{{\partial^{2}}{f(x)}}{{\partial}x^{2}}+\frac{({\Delta}x)^{3}}{3!}\frac{{\partial^{3}}{f(x)}}{{\partial}x^{3}}+\hdots \\ \\ &=&\displaystyle \sum_{n=0}^{\infty}\frac{1}{n!}\frac{{\partial^{n}}{f(x)}}{{\partial}x^{n}} \end{array}$

类似地，当 Δx = -1 时，有

(10-3)

$\begin{array}{lll} f ( x -1 ) &=&\displaystyle f ( x ) - {\Delta}x\frac{{\partial}{f(x)}}{{\partial}x}+\frac{({\Delta}x)^{2}}{2!}\frac{{\partial^{2}}{f(x)}}{{\partial}x^{2}}-\frac{({\Delta}x)^{3}}{3!}\frac{{\partial^{3}}{f(x)}}{{\partial}x^{3}}+\hdots \\ \\ &=&\displaystyle \sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\frac{{\partial^{n}}{f(x)}}{{\partial}x^{n}} \end{array}$

接下来，我们仅使用Taylor级数中的几项来计算强度差。对于一阶导数，我们只使用线性项，并且可以通过三种方式之一来表示差分。

正向差分可通过公式 (10-2)求得：

(10-4) $\displaystyle \frac{{\partial}f(x)}{{\partial}x} = f^{'}(x) = f ( x + 1) - f ( x )$

其中，正如你所见，我们仅保留了线性项。后向差分的求取也类似，仅保留公式 (10-3) 中的线性项：

(10-5) $\displaystyle \frac{{\partial}f(x)}{{\partial}x} = f^{'}(x) = f ( x ) - f ( x-1 )$

中心差分(central difference)是通过从式 (10-2) 中减去式 (10-3) 得到的：

(10-6) $\displaystyle \frac{{\partial}f(x)}{{\partial}x} = f^{'}(x) = \frac{f ( x +1) - f ( x-1 )}{2}$

级数中未使用的高阶项代表了精确导数展开式与近似导数展开式之间的误差。一般来说，Taylor 级数展开式中用于表示导数的项越多，近似结果就越精确。包含更多项意味着近似过程中使用了更多点，从而产生更小的误差。然而，对于相同数量的点，中心差分法的误差反而更小(参见习题 10.1)。因此，导数通常用中心差分法表示。

基于中心差分的 ${{\partial^{2}}{f(x)}}/{{\partial}x^{2}}$ 的二阶导数可通过将公式 (10-2) 和公式 (10-3) 相加而求得：

(10-7) $\displaystyle \frac{{\partial^{2}}f(x)}{{\partial}x^{2}} = f^{''}(x) = f ( x + 1) -2f(x) +f ( x-1 )$

为了获得三阶、中心导数，我们需要在 x 的两侧各再加一个点。即，即我们需要 f ( x + 2) 和 f ( x - 2) 的展开式，我们分别从公式 (10-2) ( 取 Δx = 2 ) 和公式 (10-3) ( 取 Δx = -2 ) 所求得。该策略是将两个Taylor展开式合并，以消除所有低于三阶的导数。忽略所有高阶项后的结果是[参见问题 10.2(a)]

(10-8)

$\displaystyle \frac{{\partial^{3}}f(x)}{{\partial}x^{3}} = f^{'''}(x) = \frac{f ( x + 2) -2f(x+1)+0f(x) +2f ( x-1)-f(x-2)}{2}$

类似地[见公式 10.2(b)]，在忽略所有更高阶项项之后，第四个有限差分(我们在本书中使用的最高阶)公式为

(10-9)

$\displaystyle \frac{{\partial^{4}}f(x)}{{\partial}x^{4}} = f^{''''}(x) =f ( x + 2) - 4 f ( x + 1 ) + 6 f ( x ) - 4 f ( x - 1) + f ( x - 2)$

表 10.1 总结了前面讨论的前四个中心导数。请注意系数关于中心点的对称性。这种对称性是中心差分在相同数据点数量下比其他两种差分方法具有更低近似误差的根本原因。对于两个变量，我们将表 10.1 中的结果分别应用于每一个变量。例如，

(10-10) $\displaystyle \frac{{\partial^{2}}f(x,y)}{{\partial}x^{2}} = f ( x + 1 ,y ) - 2f ( x , y ) + f ( x - 1, y)$

和

(10-11) $\displaystyle \frac{{\partial^{2}}f(x,y)}{{\partial}y^{2}} = f ( x ,y + 1) - 2f ( x , y ) + f ( x , y - 1)$

-------------------------------表 10.1：前四个中心数字导数(有限差分)，样本均匀采集，间隔 Δx = 1 个单位。----------------------------------------------------

不难验证，公式 (10-4) 至 (10-7) 中的一阶和二阶导数满足本节开头所述的关于一阶和二阶导数的条件。为了说明这一点，请参考图 10.2。图 10.2(a) 显示了包含多个物体、一条线和一个孤点的图像。图 10.2(b) 显示了穿过图像中心(包括孤点)的水平强度分布(扫描线)。沿着扫描线，固体物体和背景之间的强度过渡呈现出两种类型的边缘：斜坡边缘(左侧)和阶梯边缘(右侧)。正如我们稍后将讨论的，涉及细长物体(例如线)的强度过渡通常被称为屋脊边缘(roof edges)。

图 10.2(c) 展示了一个简化的轮廓，其中仅包含足够多的点，以便我们能够手动分析一阶和二阶导数在遇到点、线和物体边缘时的行为。在该图中，斜坡过渡跨越四个像素，噪声点为单个像素，线宽为三个像素，阶跃边缘的过渡发生在相邻像素之间。为简单起见，强度级别数量限制为八个。

----------------------图 10.2：(a) 图像。(b) 包含箭头所示孤点的水平强度分布图。(c) 欠采样分布图；为清晰起见添加了虚线。方框中的数字是分布图中所示点的强度值。导数分别使用公式(10-4)计算一阶导数，使用公式(10-7)计算二阶导数。-------------

当我们从左到右遍历轮廓线时，考虑一阶导数和二阶导数的性质。最初，一阶导数在强度斜坡的起始点和整个斜坡上均非零，而二阶导数仅在斜坡的起始点和结束点非零。由于数字图像的边缘类似于这种过渡，我们得出结论：一阶导数产生“厚”边缘，而二阶导数产生更薄的边缘。接下来，我们遇到孤立的噪声点。在这里，二阶导数在该点的响应幅度远大于一阶导数。这并不意外，因为二阶导数在增强急剧变化方面比一阶导数更加激进。因此，我们可以预期二阶导数比一阶导数更能增强精细细节(包括噪声)。本例中的线条相当细，因此它也属于精细细节，我们再次看到二阶导数的幅度更大。最后，请注意，在斜坡和阶梯边缘，二阶导数在进入和离开边缘时符号相反(由负到正或由正到负)。这种“双边缘”效应是定位边缘的重要特征，我们将在本节后面部分进行展示。当我们进入边缘时，二阶导数的符号还可以用来判断边缘是从亮到暗(二阶导数为负)还是从暗到亮(二阶导数为正)。

综上所述，我们得出以下结论：(1) 一阶导数通常会产生较粗的边缘。(2) 二阶导数对精细细节(例如细线、孤点和噪声)的响应更强。(3) 二阶导数在强度渐变和阶跃变化处会产生双边响应。(4) 二阶导数的符号可用于判断边缘过渡是由亮到暗还是由暗到亮。

计算图像中每一个像素位置的一阶和二阶导数的首选方法是使用空间卷积。对于图 10.3 中的 3 × 3 滤波器核，具体步骤是计算核系数与核所包围区域内强度值的乘积之和，正如我们在 3.4 节中解释的那样。也就是说，滤波器在核中心点的响应是

(10-12) $\displaystyle Z = w_{1}z_{1} + w_{2} z_{2} + \hdots+ w_{9} z_{9} = \sum_{k=1}^{9}w_{k} z_{k}$

其中， $z_{k}$ 是其空间位置对应第 k 个核系数位置的像数之强度。

----------------图 10.3：一般的 3 × 3 空间滤波器核。w 是核系数(权重)。-------

10.2.2 孤点检测(Detection of isolated points)

根据前一节得出的结论，我们知道点检测应该基于二阶导数，根据第 3.6 节的讨论，这意味着使用 Laplace 算子：

(10-13) $\displaystyle {\nabla}^{2} f ( x ,y ) = \frac{{\partial}^{2}f}{{\partial}x^{2}} + \frac{{\partial}^{2} f}{{\partial}y^{2}}$

其中，偏导数是利用公式 (10-10) 和 (10-11) 中的二阶有限差法计算的。Laplace算子则为

(10-14)

$\displaystyle {\Delta}^{2}f ( x ,y ) = f ( x + 1 , y ) + f ( x - 1, y) + f ( x , y + 1) + f ( x , y - 1) - 4 f ( x ,y )$

如第 3.6 节所述，该表达式可以使用图 10.4(a) 中示例 10.1 中的Laplace核来实现。我们称在核中心位置 (x, y) 检测到一个点，如果滤波器在该点的响应绝对值超过指定的阈值。在输出图像中，这样的点标记为 1，所有其他点标记为 0，从而生成二值图像。换句话说，我们使用以下表达式：

(10-15) $g(x,y)=\left \{ \begin{array}{l} 1 \hspace{0.2cm}(\text{ if \; } |Z(x,y)|>T) \\\\ 0 \hspace{0.2cm}(\text{ otherwise\; }) \end{array} \right.$

其中，g( x ,y ) 是输出图像，T 是一个非负阈值，而 Z 由公式 (10-12) 给出。该公式简单地度量一个像素与其8个相邻像素之间的加权差。直观地说，其思想是孤点的强度与其周围环境差异很大，因此很容易被这种类型的核函数检测到。我们感兴趣的强度差是那些足够大(由T决定)以至于可以视为孤点的差。请注意，与导数核函数一样，系数之和为零，这表明滤波器的响应在强度恒定的区域为零。

例 10.1：图像中孤点的检测。

图 10.4(b) 是喷气发动机涡轮叶片的 X 射线图像。叶片存在孔隙，图像右上象限中有一个黑色像素。图 10.4(c) 是使用Laplace核对图像进行滤波的结果，图 10.4(d) 显示了公式 (10-15) 的结果，其中 T 等于图 10.4(c) 中图像最高绝对像素值的 90%。该像素清晰可见，位于箭头尖端（为了增强其可见性，像素已被放大）。这种检测方法的特殊之处在于，它基于探测器核区域内均匀背景中单个像素位置的突变强度变化。当此条件不满足时，本章讨论的其他方法更适合检测强度变化。

-----------------图 10.4：(a) 用于点检测的Laplace核。(b) 涡轮叶片的X射线图像，其中孔隙以单个黑色像素表示。(c) 核与图像卷积的结果。(d) 使用公式(10-15)的结果为单个点(箭头尖端处放大显示)。(原始图像由X-TEK Systems, Ltd.提供)----------------

10.2.3 线检测(Line Detection)

下一个复杂度级别的检测是线检测。根据本节前面的讨论，我们知道，对于线检测，二阶导数会产生更强的滤波器响应，并且生成的线比一阶导数更细。因此，我们也可以使用图 10.4(a) 中的Laplace核进行直线检测，但需要注意的是，必须正确处理二阶导数的双线效应。以下例子说明了该过程。

例 10.2：用 Laplace 核实现线检测。

图 10.5(a) 显示了电子电路引线键合掩模的 486 × 486(二值图像)部分，图 10.5(b) 显示了其Laplace变换图像。由于Laplace变换图像包含负值(参见例 3.18 后的讨论)，因此需要进行缩放才能显示。如图所示，中间灰色代表零值，较深的灰色代表负值，较浅的灰色代表正值。在放大区域中，双线效应清晰可见。

乍一看，似乎可以通过对Laplace算子图像取绝对值来处理负值。然而，如图 10.5(c) 所示，这种方法会使线条粗细增加一倍。更合适的方法是仅使用Laplace算子的正值(在噪声较大的情况下，我们使用超过正阈值的值来消除噪声引起的零附近的随机波动)。如图 10.5(d) 所示，这种方法可以得到更细的线条，通常更实用。注意图 10.5(b) 到 (d)，当线条相对于Laplace算子核的大小较宽时，线条之间会出现一个零“谷”。这并不意外。例如，当 3 × 3 的算子核位于一条 5 像素宽的恒定强度线上时，响应将为零，从而产生上述现象。当我们讨论线条检测时，通常假设线条相对于检测子的大小较细。不符合此假设的线条最好视为区域，并由下一节讨论的边缘检测方法进行处理。

-------------------------图 10.5：(a) 原始图像。(b) Laplace算子图像；放大区域显示了Laplace算子特有的正负双线效应。(c) Laplace算子的绝对值。(d) Laplace算子的正值。----------------------------------------------------------------------

图 10.4(a) 中的Laplace 检测子核是各向同性的，因此其响应与方向无关(相对于 3 × 3 核的四个方向：垂直、水平和两条对角线)。通常，人们感兴趣的是检测特定方向上的线条。考虑图 10.6 中的核。假设一幅具有恒定背景且包含各种线条(方向分别为 0°,±45° 和 90°) 的图像，使用第一个核进行滤波。最大响应将出现在图像中水平线穿过核中间行的位置。这很容易验证：绘制一个简单的全 1 阵列，并在阵列中水平穿过一条强度不同的线(例如，强度为 5 的线)。类似的实验表明，图 10.6 中的第二个核对方向为 +45° 的线条响应最佳；第三个核对垂直线条响应最佳；第四个核对方向为 −45° 的线条响应最佳。每个核函数的优先方向的权重系数(即 2)大于其他可能的方向。每个核函数的系数之和为零，表明在强度恒定的区域响应为零。

令 $Z_{1}$ , $Z_{2}$ , $Z_{3}$ ，和 $Z_{4}$ 分别表示与图 10.6 中从左到右的核相对应的响应，其中， $Z_{s}$ 由公式 (10-12) 给出。假设一幅图像用这四个核进行滤波，一次处理一个。若在图像中的一个已知点有 $| Z_{k}| > | Z_{j}|$ (对于所有 j ≠ k ) ，则称这一点为在核 k 的方向与线更可能相关的点。例如，若图像中的一个已知点满足 $|Z_{1}| > | Z_{j}|$ (对于j = 2 ，3 ，4 ) ，则称这个点为更可能与一条水平线相关的点。如果我们想要检测图像中所有沿给定核函数方向的线条，只需将该核函数应用于图像，并对结果的绝对值进行阈值处理，如公式 (10-15) 所示。阈值处理后剩余的非零点即为最强的响应，对于单像素粗的线条而言，这些响应最接近核函数定义的方向。以下示例说明了此过程。

------------------------------图 10.6：线检测核。检测角度是相对于图 2.19 中的坐标系而言的，正角度是相对于(垂直)x 轴逆时针测量的。---------------------------

例 10.3：按指定方向检测线条。

图 10.7(a) 显示了前一个例子中使用的图像。假设我们想要找到所有宽度为 1 个像素且方向为 +45° 的线条。为此，我们使用图 10.6(b) 中的卷积核。图 10.7(b) 是使用该卷积核对图像进行滤波的结果。与之前一样，图 10.7(b) 中比灰色背景更深的阴影对应于负值。图像中有两个主要线段方向为 +45°，分别位于左上角和右下角。图 10.7(c) 和 (d) 显示了图 10.7(b) 中对应于这两个区域的放大图。图 10.7(d) 中的线段比图 10.7(c) 中的线段更亮，因为图 10.7(a) 右下角的线段宽度为 1 个像素，而左上角的线段宽度则不是 1 个像素。该核函数经过“调谐”，能够检测 +45°方向上的单像素粗细的线条，因此我们预期当检测到此类线条时，其响应会更强。图10.7(e)显示了图10.7(b)中的正值。由于我们关注的是最大响应，因此我们令T等于254(图10.7(e)中的最大值减1)。图10.7( f )以白色显示了满足条件g > T的点，其中g是图10.7(e)中的图像。图中的孤立点也是对该核函数具有类似强响应的点。在原图像中，这些点及其相邻点的方向使得核函数在这些位置产生最大响应。这些孤立点可以使用图10.4(a)中的核函数进行检测并删除，或者可以使用形态学算子进行删除，如上一章所述。

------------------------------图 10.7：(a) 引线键合模板图像。(b) 使用图 10.6 中的 +45°线检测子内核处理后的结果。(c) (b) 左上区域的放大图。(d) (b) 右下区域的放大图。(e) 将 (b) 中的所有负值设为零后的图像。(f) 所有满足条件 g > T 的点(白色)，其中 g 为 (e) 中的图像，T = 254(图像中的最大像素值减 1)。( ( f ) 中的点已放大以便于观察。)----------------------------

10.2.4 边缘模型(Edge Models)

边缘检测是一种常用的图像分割方法，它基于图像中强度的突变(局部)变化。我们首先介绍几种边缘建模方法，然后讨论一些边缘检测方法。

边缘模型根据其强度分布进行分类。阶跃边缘的特征是两个强度级别之间的过渡，理想情况下，该过渡发生在一个像素的距离内。图 10.8(a) 显示了垂直阶跃边缘的截面以及穿过该边缘的水平强度分布。例如，在计算机生成的用于实体建模和动画等领域的图像中，就会出现阶跃边缘。只要不使用任何额外的处理(例如平滑处理)来使其看起来“真实”，这些清晰、理想的边缘就可以在一个像素的距离内出现。数字阶跃边缘经常用作算法开发中的边缘模型。例如，本节稍后讨论的 Canny 边缘检测算法最初就是基于阶跃边缘模型推导出来的。

实际上，数字图像的边缘往往模糊且存在噪声，模糊程度主要取决于聚焦机制(例如光学图像中的镜头)的限制，而噪声水平则主要取决于成像系统的电子元件。在这种情况下，边缘更接近于具有强度斜坡轮廓的模型，例如图 10.8(b) 中的边缘。斜坡的斜率与边缘的模糊程度成反比。在这个模型中，我们不再沿着轮廓线定义一个单一的“边缘点”。相反，边缘点现在是斜坡上的任何一点，而边缘线段则是一组相互连接的此类点。

第三种边缘类型是所谓的屋脊边缘，其特征如图 10.8(c) 所示。屋脊边缘是穿过某个区域的线条模型，其底边(宽度)由线条的粗细和清晰度决定。极限情况下，当其底边宽度为一个像素时，屋脊边缘就相当于图像中穿过某个区域的一条像素粗的线条。例如，在深度成像中，当细长物体(例如管道)比背景(例如墙壁)更靠近传感器时，就会出现屋脊边缘。管道看起来更亮，因此会形成类似于图 10.8(c) 中模型的图像。屋脊边缘也经常出现在其他领域，例如线图数字化和卫星图像中，在这些领域中，细长特征(例如道路)可以用这种边缘类型来建模。

-----------------------图 10.8：从左到右依次为阶跃、斜坡和屋顶边缘的模型(理想表示)及其相应的强度分布图。-------------------------------------------------

同时包含三种边缘类型的图像并不罕见。尽管模糊和噪声会导致边缘形状偏离理想形状，但清晰度较高且噪声适中的图像中的边缘确实与图 10.8 中的边缘模型特征相似，如图 10.9 中的轮廓所示。图 10.8 中的模型使我们能够为图像处理算法开发中的边缘编写数学表达式。这些算法的性能取决于实际边缘与算法开发中所用模型之间的差异。

-----------------------图 10.9：一张 1508 × 1970 像素的图像，显示了(放大后的)实际斜波(左下)、阶跃(右上)和屋顶边缘的轮廓。轮廓线由深至浅排列，位于小圆圈内。斜波和阶跃的轮廓线分别跨越 9 像素和 2 像素。屋顶边缘的底部为 3 像素。(原图由范德比尔特大学(Vanderbilt University)的 David R. Pickens 博士提供。)----------------

图 10.10(a) 显示了提取图 10.8(b) 中线段的图像。图 10.10(b) 显示了一条水平强度分布图。该图还显示了强度分布图的一阶导数和二阶导数。沿着强度分布图从左向右移动，我们注意到，在斜坡的起始点和斜坡上的某些点处，一阶导数为正；在强度恒定的区域，一阶导数为零。二阶导数在斜坡的起始点处为正；在斜坡的结束点处为负；在斜坡上的某些点处以及强度恒定的区域，二阶导数为零。对于由亮到暗过渡的边缘，上述导数的符号会相反。零强度轴与连接二阶导数极值的直线的交点称为二阶导数的零交叉点。

------------------------图 10.10：(a) 两个强度恒定的区域由理想的斜坡边缘隔开。(b) 边缘附近的细节，显示了水平强度分布及其一阶和二阶导数。-------------------

从这些观察结果中我们可以得出结论：一阶导数的大小可用于检测图像中某一点是否存在边缘。类似地，二阶导数的符号可用于确定边缘像素位于边缘的暗侧还是亮侧。二阶导数在边缘附近的另外两个特性是：(1) 对于图像中的每一条边缘，它都会产生两个值；(2) 它的零交叉点可用于定位较厚边缘的中心，我们将在本节后面部分进行展示。一些边缘模型利用平滑的过渡进入和退出斜坡(参见问题 10.9)。然而，使用这些模型得出的结论与使用理想斜坡得出的结论相同，而且使用后者可以简化理论公式。最后，尽管到目前为止，我们只关注了一维水平轮廓，但类似的论证也适用于图像中任意方向的边缘。我们只需在任意所需点定义一个垂直于边缘方向的轮廓，并以与刚才讨论的垂直边缘相同的方式解释结果即可。

例 10.4：噪声边缘区域的一阶导数和二阶导数的行为。

图 10.8 中的边缘模型不含噪声。图 10.11 第一列的图像片段显示了四个斜坡边缘的特写，这些边缘从左侧的黑色区域过渡到右侧的白色区域(请注意，从黑色到白色的整个过渡过程是一个单独的边缘)。左上角的图像片段不含噪声。第一列中的其他三个图像分别添加了均值为零、标准差别为 0.1、1.0 和 10.0 强度等级的加性Gauss噪声。每一个图像下方的图表是穿过图像中心的水平强度分布图。所有图像的强度分辨率均为 8 位，其中 0 和 255 分别代表黑色和白色。

考虑中间列顶部的图像。如结合图 10.10(b) 所述，左侧扫描线的导数在恒定区域为零。这些恒定区域在导数图像中显示为两条黑色条带。斜坡上各点的导数为常数，且等于斜坡的斜率。导数图像中的这些常数值以灰色显示。随着我们向下移动中间列，导数与无噪声情况的差异越来越大。事实上，很难将中间列中的最后一个轮廓与斜坡边缘的一阶导数关联起来。这些结果的有趣之处在于，左侧列图像中的噪声几乎无法用肉眼察觉。这些例子很好地说明了导数对噪声的敏感性。

正如预期的那样，二阶导数对噪声更加敏感。无噪声图像的二阶导数显示在右侧列的顶部。如图 10.10 所示，细白线和细黑线分别代表二阶导数的正分量和负分量。这些图像中的灰色表示零(如前所述，缩放会导致零显示为灰色)。唯一一张几乎与无噪声情况相似的含噪声二阶导数图像对应于标准差为 0.1 的噪声。其余的二阶导数图像和轮廓清晰地表明，检测它们的正分量和负分量确实非常困难，而这些分量才是二阶导数在边缘检测方面真正有用的特征。

视觉噪声如此之小，却能对用于边缘检测的两个关键导数产生如此显著的影响，这一点值得我们重视。尤其是在应用中可能存在类似我们刚才讨论的噪声水平时，图像平滑处理应作为使用这些导数之前需要认真考虑的问题。

----------------------------------图 10.11：第一列：取值范围为 [0, 255] 的 8 位图像，以及分别受均值为零、标准差为 0.0，0.1，1.0 和 10.0 强度级的Gauss噪声干扰的斜坡边缘强度分布图。第二列：一阶导数图像和强度分布图。第三列：二阶导数图像和强度分布图。-----------------------------------------------------------------

总之，边缘检测通常包括以下三个步骤：

(1) 籍图像平滑降噪。图 10.11 的第二列和第三列的结果说明了这一步骤的必要性。

(2) 检测边缘点。如前所述，这是一个局部操作，从图像中提取所有潜在的边缘点之候选点。

(3) 边缘局部化。此步骤的目标是从候选点中选择只属于构成边之点集的点。

本节余下部分将探讨实现这些目标的技术。

10.2.5 基本边缘检测(Basic Edge Detection)

如前所述，检测强度变化以寻找边缘可以使用一阶或二阶导数。我们首先讨论一阶导数，并在下一小节中讨论二阶导数。

10.2.5.1 图像梯度及其属性(The Image Gradient and It’s Properties)

求得一幅图像 f 在一个任意位置 ( x ,y ) 处的的边缘长度和方向的选择工具是梯度，表示为 ▽f ，并定义为向量

(10-16) $\displaystyle {\nabla}f( x,y ) \equiv \text{grad}[ f ( x ,y ) ] \equiv \begin{bmatrix} g_{x}(x,y)\\\\ g_{y}(x,y) \end{bmatrix} \equiv \begin{bmatrix} \frac{\displaystyle {\partial}f(x,y)}{\displaystyle {\partial}x}\\ \\ \frac{\displaystyle { \displaystyle \partial}f(x,y)}{ \displaystyle {\partial}y} \end{bmatrix}$

向量有一个著名的属性，即其指向f 在位置 ( x ,y ) 处最大变化率的方向(见问题 10.10)。公式 (10-16) 在一个任意点(但是单点)( x ,y ) 处是有效的。当对所有适用的 x 和 y 值进行计算时， $\nabla f ( x ,y )$ 就成了一幅图像，其每一个元素都是由公式 (10-16)给出的向量。这个梯度向量在点 ( x ,y ) 处的大小 M (x ,y ) (Magnitude)通过其 Euclid 向量范数进行计算：

(10-14) $\displaystyle M(x ,y) = {\Vert}{\nabla} f ( x ,y ){\Vert} = \sqrt{g_{x}^{2} (x,y)+g_{y}^{2}(x,y) }$

这个值是梯度向量在点 ( x ,y ) 方向上的改变率之值。注意， M (x ,y) ， ${\nabla}f ( x ,y )$ ， $g_{x}(x ,y)$ 和 $g_{y}(x ,y)$ 是与 f 大小相同的数组，当 x 和 y 允许在 f 上的所有像素位置处变化时可创立这些量。当这个含义清晰时，在常规实践中称 M(x ,y) 和 $\nabla f ( x ,y )$ 为梯度图像，或直接简称梯度。求和、平方、和平方根运算是逐元素进行的，正如 2.6节中所定义。

在点 (x ,y) 处的梯度向量的方向的计算公式为

(10-18) $\displaystyle \alpha(x,y) = {\tan}^{-1} \left [ \frac{g_{y}(x ,y )}{g_{x}(x , y)} \right ]$

角度是针对 x 轴按逆时针方向度量的(见图 2.19)。这也是一幅与 f 大小相同的图像，通过在所有 x 和 y 的适用值上 $g_{x}$ 和 $g_{y}$ 按逐元素相除而创建。下述例子说明，边缘在 (x ,y) 处的方向正交于梯度向量在这点的 α (x ,y ) 。

例 10.5：计算梯度。

图 10.12(a) 显示了一幅图像的放大区域，其中包含一段直线边缘。每一个方框对应一个像素，我们感兴趣的是获取图中方框高亮点处边缘的强度和方向。图中阴影像素的值为 0，白色像素的值为 1。在本例之后，我们将讨论一种利用以该点为中心的 3 × 3 邻域计算 x 和 y 方向导数的方法。该方法包括从邻域底部一行的像素值中减去顶部一行的像素值，从而获得 x 方向的偏导数。类似地，我们从邻域右侧一行的像素值中减去左侧一行的像素值，从而获得 y 方向的偏导数。然后，利用这些差值作为偏导数的估计值，可以得出在目标点处 ∂f/∂x = -2 和 ∂ f/∂y = 2 。则

$\displaystyle {\nabla}f {\equiv} \begin{bmatrix}g_{x}\\g_{y} \end{bmatrix} {\equiv}\begin{bmatrix}\frac{{\partial}f}{{\partial}x}\\ \\ \frac{{\partial}f}{{\partial}x} \end{bmatrix} =\begin{bmatrix}-2\\2 \end{bmatrix}$

根据这个计算值，在这个点处求得 ${\Vert}{\nabla}f{\Vert} = 2\sqrt{2}$ 。类似地，根据公式 (10-18)求得在同一点的梯度方向为： $\alpha = {\tan}^{-1}({g_{y}}/g_{x}) = -45^{\circ}$ ，这相当于在我们的图像坐标系中，沿 x 轴正方向(逆时针方向)测量的 135°(见图 2.19)。图 10.12(b) 显示了梯度矢量及其方向角。

如前所述，某点处边缘的方向与该点的梯度向量正交。因此，本例中边缘的方向角为 -90°= 135°- 90°= 45°，如图 10.12(c) 所示。图 10.12(a) 中的所有边缘点都具有相同的梯度，因此整个边缘线段的方向相同。梯度向量有时称为边缘法向量(edge normal)。当该向量除以其模长进行归一化后，所得向量称为边缘单位法向量(edge unit normal)。

----------------------图 10.12：利用梯度确定某一点的边缘强度和方向。注意，边缘方向垂直于计算梯度的点处的梯度向量方向。每个方格代表一个像素。(回顾图 2.19，坐标系原点位于左上角。)---------------------------------------------------------

10.2.5.2 梯度算子(Gradient Operators)

要获得图像的梯度，需要计算图像中每一个像素位置的偏导数 ∂f/∂x 和 ∂f/∂y 。对于梯度计算，我们通常使用前向或中心差法(参见表 10.1)。使用前向差我们求得

(10-19) $\displaystyle g_{x}( x ,y ) = \frac{{\partial}f(x ,y )}{{\partial}x} = f ( x + 1 ,y ) - f ( x ,y )$

和

(10-20) $\displaystyle g_{y}( x ,y ) = \frac{{\partial}f(x,y)}{{\partial}y}= f ( x ,y + 1 ) - f ( x ,y )$

这两个公式可以通过使用图 10.13 中的一维核对 f ( x ,y ) 进行滤波来实现，适用于 x 和 y 的所有值。

-----------图 10.13：用于实现公式 (10-19) 和 (10-20) 的一维核函数。-----------

当需要关注对角边缘方向时，我们需要二维核函数。Roberts 交叉梯度算子( Roberts [1965])是最早尝试使用具有对角偏好的二维核函数的方法之一。考虑图 10.14(a) 中的 3 × 3 区域。Roberts 算子是基于实现对角差

(10-21) $\displaystyle g_{x} = \frac{{\partial}f}{{\partial}x} = ( z_{9} - z_{5} )$

和

(10-22) $\displaystyle g_{y} = \frac{{\partial}f}{{\partial}y} = ( z_{8} - z_{6} )$

这些导数可以通过使用图 10.14(b) 和 (c) 中所示的核对图像进行滤波来实现。

-------------------------图 10.14：图像的 3 × 3 区域( z 为强度值)，以及用于计算标记为 $z_{5}$ 的点的梯度的各种核。----------------------------------------------

2×2 大小的核在概念上很简单，但它们在计算边缘方向方面不如关于其中心对称的核有用，其中最小的核的大小为 3×3。这些核考虑了中心点两侧数据的性质，因此携带了更多关于边缘方向的信息。用大小为 3×3 的核计算的最简单的偏导数数字近似值公式为

(10-23) $\displaystyle g_{x} = \frac{{\partial}f}{{\partial}x} = ( z_{7} + z_{8} + z_{9} ) - ( z_{1} + z_{2} + z_{3} )$

和

$\displaystyle g_{y} = \frac{{\partial}f}{{\partial}y} = ( z_{3} + z_{6} + z_{9} ) - ( z_{1} + z_{4} + z_{7} )$

在此公式中，3×3 区域中第三行与第一行之差近似表示 x 方向的导数，第三列与第一列之差近似表示 y 方向的导数。直观上，我们预期这些近似值比使用 Roberts 算子得到的近似值更精确。公式 (10-22) 和 (10-23) 可以通过使用图 10.14(d) 和 (e) 中的两个核函数对整幅图像进行滤波来实现。这些核函数称为 Prewitt 算子(Prewitt [1970])。

前两个方程稍作修改，中心系数的权重设为 2：

(10-24) $\displaystyle g_{x} = \frac{{\partial}f}{{\partial}x} = ( z_{7} + 2z_{8} + z_{9} ) - ( z_{1} + 2z_{2} + z_{3} )$

和

(10-25) $\displaystyle g_{y} = \frac{{\partial}f}{{\partial}y} = ( z_{3} + 2z_{6} + z_{9} ) - ( z_{1} + 2z_{4} + z_{7} )$

可以证明(参见习题 10.12)，在中心位置使用 2 可以实现图像平滑。图 10.14( f ) 和 ( g ) 显示了用于实现公式 (10-24) 和 (10-25) 的核函数。这些核函数被称为 Sobel 算子( Sobel [1970] )。

Prewitt 核函数比 Sobel 核函数更容易实现，但它们之间微小的计算差异通常不会造成问题。Sobel 核函数具有更好的噪声抑制(平滑)特性，因此更受欢迎，因为正如前面在图 10.11 的讨论中提到的，噪声抑制在处理导数时是一个重要问题。请注意，图 10.14 中所有核函数的系数之和为零，因此在强度恒定的区域响应为零，这与导数算子的预期一致。

图 10.14 中的任意一对卷积核与图像进行卷积，以获得每一个像素位置的梯度分量

$g_{x}$ 和 $g_{y}$ 。这两个偏导数数组随后用于估计边缘强度和方向。获得梯度幅值需要进行公式 (10-17) 中的计算。这种实现方式并非总是理想的，因为平方和平方根运算会带来较大的计算负担。因此，一种常用的方法是使用绝对值来近似梯度幅值:

(10-26) $M( x ,y ) \approx | g_{x} | + | g_{y} |$

该方程在计算上更具优势，并且仍然保留了强度级别的相对变化。这种优势的代价是，得到的滤波器通常不是各向同性的(isotropic)(旋转不变的)。然而，当使用诸如 Prewitt 和 Sobel 核之类的核来计算 $g_{x}$ 和 $g_{y}$ 时，这不是问题，因为这些核仅对垂直和水平边缘给出各向同性的结果。这意味着无论使用哪个方程，结果都只会在这两个方向上的边缘是各向同性的。也就是说，当使用 Sobel 核或 Prewitt 核时，公式 (10-17) 和 (10-26) 对垂直和水平边缘给出相同的结果(参见问题 10.11)。

图 10.14 中的 3 × 3 卷积核主要对垂直和水平边缘表现出最强的响应。图 10.15 中的 Kirsch 罗盘卷积核( Kirsch [1971])旨在检测所有八个罗盘方向上的边缘大小和方向(角度)。Kirsch 的方法并非使用公式 (10-17) 计算大小，使用公式 (10-18) 计算角度，而是通过将图像与所有八个卷积核进行卷积来确定边缘大小，并将图像中某一点的边缘大小赋给在该点产生最大卷积值的卷积核的响应。该点的边缘角度则对应于与该卷积核相关的方向。例如，如果图像中某一点的最大卷积值是由北 ( N ) 卷积核产生的，则该点的边缘大小将被赋给该卷积核的响应，方向则为 0°( 因为罗盘卷积核对之间存在 180° 的旋转；选择最大响应值总是会得到一个正数)。虽然在使用 Sobel 核函数时，我们通常将南北边缘视为垂直方向，但 N 和 S 罗盘核函数能够区分南北边缘，区别在于定义边缘的强度过渡方向。例如，假设强度值在 [0, 1] 范围内，图 10.8(a) 中的二值边缘由左侧的黑色 (0) 和右侧的白色 (1) 定义。当所有 Kirsch 核函数都应用于该边缘时，N 核函数将产生最高值，从而表明该边缘指向北方(在计算点处)。

-----------图 10.15：Kirsch 罗盘核。每一个核的最强响应方向标注在其下方。-------

例 10.6：二维梯度幅值和角度示意图。

图 10.16 展示了梯度的两个分量 $g_{x}$ 和 $g_{y}$ 的 Sobel 绝对值响应，以及由这两个分量之和构成的梯度图像。梯度水平分量和垂直分量的方向性在图 10.16(b) 和 (c) 中清晰可见。需要注意的是，回顾问题 3.32 中的一个重要结论，即使用系数之和为零的核函数会生成像素之和也为零的滤波图像。这通常意味着某些像素将为负值。类似地，如果核函数系数之和为 1，则原始图像和滤波图像中的像素之和将相同(参见问题 3.31)。例如，图 10.16(b) 中屋顶瓦片、水平砖缝和窗户的水平部分与其他边缘相比显得多么突出。相比之下，图 10.16(c) 则突出了立面和窗户等垂直分量的特征。当图像的主要特征是边缘时，例如梯度幅值图像，通常使用“边缘图”一词。图 10.16(a) 中图像的强度值被缩放到 [0, 1] 范围内。我们使用此范围内的值来简化本节讨论的各种边缘检测方法中的参数选择。

-------------------图 10.16：(a)尺寸为 834 × 1114 像素的图像，强度值缩放到 [0,1] 范围内。(b) $g_{x}$ ，即 x 方向的梯度分量，使用图 10.14( f ) 中的 Sobel 核对图像进行滤波得到。(c) $g_{y}$ ，使用图 10.14(g) 中的核得到。( d ) 梯度图像， $g = | g_{x} | + | g_{y} |$ 。 -------------------------------------------------

图 10.17 显示了使用公式 (10-18) 计算得到的梯度角度图像。一般来说，角度图像在边缘检测方面不如梯度幅值(magnitude)图像有效，但它们可以补充从图像中提取的梯度幅值信息。例如，图 10.16(a) 中的等强度区域，如倾斜屋顶的前缘和前墙顶部的水平条带，在图 10.17 中也保持不变，这表明这些区域中所有像素位置的梯度矢量方向相同。正如我们将在本节后面展示的那样，角度信息在 Canny 边缘检测算法(一种广泛使用的边缘检测方案)的实现中起着关键的辅助作用。

-------------------图 10.17：利用公式(10-18)计算梯度角图像。该图像中强度恒定的区域表明，在这些区域内，梯度矢量的方向在所有像素位置上都相同。----------------

图 10.16(a) 中的原始图像分辨率相当高，在图像采集距离下，墙体砖块对图像细节的贡献显著。这种精细的细节在边缘检测中通常是不利的，因为它容易产生噪声，而导数计算会增强噪声，从而增加主要边缘检测的难度。一种减少精细细节的方法是，在计算边缘之前对图像进行平滑处理。图 10.18 显示了与图 10.16 相同的图像序列，但首先使用 5 × 5 平均滤波器对原图像进行了平滑处理(关于平滑滤波器，请参见 3.5 节)。现在，每一个卷积核的响应几乎不再受砖块的影响，结果主要由图像中的主要边缘决定。

图 10.16 和 10.18 显示，水平和垂直 Sobel 核无法区分 ±45° 方向的边缘。如果需要强调特定对角方向的边缘，则应使用图 10.15 中的 Kirsch 核之一。图 10.19(a) 和 (b) 分别显示了 45° (西北) 和 -45° (西南) Kirsch 核的响应。这些核对对角线方向的选择性更强，这一点在这些图中显而易见。两种核对水平和垂直边缘的响应相似，但垂直方向的响应较弱。

-------------------------图 10.18：与图 10.16 中的序列相同，但在边缘检测之前，使用 5 × 5 平均核对原图像进行平滑处理。------------------------------------

-------------------------图 10.19：对角边缘检测。(a) 使用图 10.15(c) 中的 Kirsch 核的结果。(b) 使用图 10.15(d) 中的核的结果。两种情况下的输入图像均为图 10.18(a)。-------------------------------------------------------------------

10.2.5.3 梯度与阈值结合法(Combining the Gradient with Thresholding)

图 10.18 的结果表明，在计算梯度之前对图像进行平滑处理可以提高边缘检测的选择性。另一种实现相同目标的方法是将梯度图像进行阈值处理。例如，图 10.20(a) 显示了图 10.16(d) 中的梯度图像，该图像经过阈值处理，使得像素值大于或等于梯度图像最大值 33% 的像素显示为白色，而低于阈值的像素显示为黑色。将此图像与图 10.16(d) 进行比较，我们可以看到阈值处理后的图像中边缘数量更少，并且边缘更加清晰锐利(例如，屋顶瓦片的边缘)。另一方面，阈值处理后的图像中许多边缘，例如定义屋顶远边缘的斜线(见箭头所示)，出现了断裂。

当既要突出主要边缘又要尽可能保持连通性时，通常会同时使用平滑和阈值分割。图 10.20(b) 显示了对图 10.18(d) 进行阈值分割的结果，即平滑图像的梯度。该结果表明断裂边缘的数量有所减少；例如，比较图 10.20(a) 和 (b) 中箭头所指的对应边缘。

------------------------图 10.20：(a)图 10.16(d) 阈值处理后的结果，原图像的梯度。(b) 图 10.18(d) 阈值处理后的结果，平滑后图像的梯度。---------------------

10.2.6 边缘检测高级技术(More Advanced Techniques for Edge Detection)

前几节讨论的边缘检测方法基于使用一个或多个卷积核对图像进行滤波，而没有考虑边缘特征和噪声含量。本节将讨论更高级的技术，这些技术通过考虑图像噪声和边缘本身的特性等因素，力求改进简单的边缘检测方法。

10.2.6.1 Marr-Hildreth边缘检测子(The Marr-Hildreth Edge Detector)

最早将更复杂的分析方法成功融入边缘检测过程的尝试之一，归功于 Marr 和 Hildreth [1980]。当时使用的边缘检测方法基于小型算子，例如前面讨论过的 Sobel 核。Marr 和 Hildreth 指出：(1) 强度变化与图像尺度并非无关，这意味着检测强度变化需要使用不同大小的算子；(2) 强度突变会导致一阶导数出现峰值或谷值，或者等价地，导致二阶导数出现零交叉点(如图 10.10 所示)。

这些想法表明，用于边缘检测的算子应具备两个显著特征。首先，它应该是一个微分算子，能够计算图像中每一个点的一阶或二阶导数的数字近似值。其次，它应该能够进行“调整”，以适应任何所需的尺度，从而可以使用较大的算子来检测模糊的边缘，而使用较小的算子来检测清晰的细节。

Marr 和 Hildreth 表明，满足这些条件的最令人满意的算子是滤波器 $\nabla^{2}G$ ，其中，如 3.6 节中所定义， $\nabla^{2}$ 是Laplace 算子，而 G 是二维 Gauss 函数

(10-27) $G( x , y ) = e^{\displaystyle -\frac{x^{2}+y^{2}}{2{\sigma}^{2}}}$

其中，σ 是标准偏差(在这个背景下，有时候也称 σ 为空间常量)。我们将 Laplace 算子应用于公式 (10-27) 来求得 $\nabla^{2}G$ 的一种表达式：

(10-28) $\begin{array}{lll} \nabla^{2}G(x,y)&=&\displaystyle \frac{{\partial}^{2}G(x,y)}{{\partial}x^{2}}+\frac{{\partial}^{2}G(x,y)}{{\partial}y^{2}} \\ \\ &=&\displaystyle \frac{\partial}{{\partial}x}\left ( \frac{-x}{\sigma^{2}} e^{\displaystyle -\frac{x^{2}+y^{2}}{2{\sigma}^{2}}} \right )+\frac{\partial}{{\partial}y} \left ( \frac{-y}{\sigma^{2}} e^{\displaystyle -\frac{x^{2}+y^{2}}{2{\sigma}^{2}}} \right ) \\ \\ &=&\displaystyle \left (\frac{x^{2}}{\sigma^{4}}-\frac{1}{\sigma^{2}} \right ) e^{\displaystyle -\frac{x^{2}+y^{2}}{2{\sigma}^{2}}} + \left (\frac{y^{2}}{\sigma^{4}}-\frac{1}{\sigma^{2}} \right )e^{\displaystyle -\frac{x^{2}+y^{2}}{2{\sigma}^{2}}} \end{array}$

合并项，我们求得

(10-29) $\begin{array}{lll} \nabla^{2}G(x,y)=\left ( \displaystyle \frac{x^{2}+y^{2}-2\sigma^{2}}{\sigma^{4}} \right )e^{ \displaystyle -\frac{x^{2}+y^{2}}{2\sigma^{2}}} \end{array}$

这人表达式称为 Gauss 函数的 Laplace 算子(Laplacian of a Gaussian)(LoG)。

图 10.21(a) 至 (c) 分别展示了 LoG 函数负值的三维图、图像和截面图(注意 LoG 函数的零交叉点位于 $x^{2}+y^{2}=2\sigma^{2}$ ，这定义了一个以 Gauss 函数峰值为中心、半径为 $\sqrt{2}\sigma$ 的圆)。由于图 10.21(a) 所示的形状，LoG 函数有时也称为Mexico帽算子(Mexico hat operator)。图 10.21(d) 展示了一个 5 × 5 的核函数，它近似于图 10.21(a) 中的形状(通常情况下，我们会使用该核函数的负值）。这种近似并非唯一。其目的是捕捉 LoG 函数的基本形状；就图 10.21(a) 而言，这意味着一个正的中心项被一个相邻的负值区域包围，该负值区域的值随着与原点距离的增加而减小，以及一个零值外层区域。系数之和必须为零，这样核在强度恒定的区域的响应才为零。

-------------------------------图 10.21：(a) LoG 负像(Negative)的三维图。(b) LoG 负像的图像显示。(c) (a) 的横截面，显示零交叉点。(d) (a) 中形状的 5 × 5 核近似。实际应用中将使用该核的负像。---------------------------------------------

可以通过对公式 (10-29) 进行采样，并将系数缩放至总和为零，来生成任意大小(但 σ 固定)的滤波器核。生成 LoG 核的更有效方法是对公式 (10-27) 进行采样至所需大小，然后将所得数组与Laplace核(例如图 10.4(a) 中的核)进行卷积。由于将图像与系数总和为零的核进行卷积会得到元素总和也为零的图像(参见习题 3.32 和 10.16)，因此该方法自动满足 LoG 核系数之和为零的要求。我们将在本节后面讨论 LoG 滤波器的大小选择。

选择算子 $\nabla^{2}G$ 背后有两个基本思路。首先，该算子的 Gauss 部分会模糊图像，从而降低远小于σ 尺度的结构(包括噪声)的强度。与图10.18中使用的平均滤波器不同，Gauss函数在空间域和频域中都是平滑的(参见4.8节)，因此不太可能引入原图像中不存在的伪影(例如振铃)。另一个思路与 Laplace 算子 $\nabla^{2}$ 的二阶导数特性有关。虽然一阶导数可以用于检测强度的突变，但它们是方向性算子。在另一方面，Laplace 算子具有各向同性(旋转不变)的重要优势，这不仅符合人类视觉系统的特征( Marr [1982] )，而且对任何核方向上的强度变化都做出相同的响应，从而避免了使用多个核来计算图像中任何点的最强响应。

Marr-Hildreth 算法由 LoG 核与一个输入图像的卷积构成，即

(10-30) $g(x ,y) = [\nabla^{2}G( x , y )] {\bigstar} f( x , y )$

然后求得 g(x, y) 的过零点(或零交叉点)(zero crossings)，从而确定 f ( x , y ) 中边缘的位置。因为 Laplace 算子和卷积是线性过程，我们可以将 (10-30) 写成

(10-31) $g(x ,y) = \nabla^{2}[G( x , y ) {\bigstar} f ( x , y )]$

这表明我们可以先用Gauss滤波器平滑图像，然后计算结果的Laplace算子。这两个公式给出的结果完全相同。

Marr-Hildreth边缘检测算法可概括如下：

(1) 通过使用采样公式(10-27)而获得的 n × n Gauss低通核对输入图像进行滤波。

(2) 使用例如图 10.4(a) 中的 3 × 3 核，计算步骤 (1) 所得图像的Laplace算子。[步骤 (1) 和 (1) 实现了公式 (10-31)。]

(3) 求得步骤 2 中图像的零交叉点。

为了确定Gauss核的大小，回顾图 3.35 的讨论，Gauss函数在距离均值大于 3σ 处的值足够小，可以忽略不计。如 3.5 节所述，这意味着使用大小为 [6σ] × [6σ] 的 Gauss 核，其中 [6σ] 表示 6σ 的上限，即不小于 6σ 的最小整数。由于我们处理的是奇数维核，因此应使用满足此条件的最小奇数。使用小于此值的核会“截断”LoG 函数，截断程度与核的大小成反比。使用更大的核对结果几乎没有影响。

求取滤波图像 g(x ,y) 中任意像素 p 的零交叉点的一种方法是以 p 为中心，构建一个 3 × 3 的邻域。p 处的零交叉点意味着其至少两个相对相邻像素的符号不同。需要测试四种情况：左右、上下以及两条对角线。如果将 g(x ,y) 的值与阈值进行比较(一种常用方法），那么不仅相对相邻像素的符号必须不同，而且它们的数值差的绝对值也必须超过阈值，才能将 p 判定为零交叉点像素。我们在例 10.7 中演示了这种方法。

计算零交叉点是 Marr-Hildreth 边缘检测方法的关键特征。前一段讨论的方法因其实现简单且通常能取得良好结果而备受青睐。如果该方法在特定应用中找到的零交叉点位置精度不足，则可以采用 Huertas 和 Medioni [1986] 提出的亚像素(subpixel)精度零交叉点查找技术。

例 10.7：Marr-Hildreth边缘检测方法示意图。

图 10.22(a) 显示了之前使用的建筑图像，图 10.22(b) 是 Marr-Hildreth 算法步骤 1 和步骤 2 的结果，其中 σ = 4(约为图像短边的 0.5% ) 和 n = 25 以满足上述尺寸条件。与图 10.5 类似，该图像中的灰度是由于缩放造成的。图 10.22(c) 显示了使用前面讨论的 3 × 3 邻域方法，阈值为零时获得的零交叉点。请注意，所有边缘都形成了闭合回路。这种所谓的“意大利面条效应”是该方法在使用零阈值时的一个严重缺陷(参见问题 10.17)。我们通过使用正阈值来避免闭合回路边缘。

图 10.22(d) 展示了使用阈值约为 LoG 图像最大值 4% 的结果。大部分主要边缘都能被轻松检测到，而砖块和瓦片屋顶等“无关”特征则被滤除。这种性能几乎无法通过前面讨论的基于梯度的边缘检测技术获得。使用零交叉点进行边缘检测的另一个重要优势是，得到的边缘宽度仅为 1 像素。这一特性简化了后续的处理步骤，例如边缘连接。

--------------------图10.22：(a) 尺寸为 834 × 1114 像素的图像，强度值缩放至 [0, 1] 范围。(b) 使用 σ = 4 和 n = 25 的 Marr-Hildreth 算法步骤 (1) 和步骤 (2) 的结果。(c) 使用阈值 0 对 (b) 进行过零点检测(注意闭环边缘)。(d) 使用等于 (b) 中图像最大值 4% 的阈值找到的过零点。注意细边缘。-------------------------------

可以用Gauss差( DoG——difference of Gaussians)来近似公式(10-29)中的LoG函数：

(10-32)

$\displaystyle D_{G}( x , y ) = \frac{1}{2{\pi}{\sigma}_{1}^2} e^{\displaystyle -\frac{x^{2}+y^{2}}{2\sigma_{1}^{2}}} - \frac{1}{2{\pi}{ \sigma}_{2}^2} e^{\displaystyle -\frac{x^{2}+y^{2}}{2\sigma_{2}^{2}}} \hspace{0.2cm} ( \sigma_{1} > \sigma_{2} )$

实验结果表明，人类视觉系统中的某些“通道(channels)”对方向和频率具有选择性，并且可以使用标准差比为 1.75:1 的公式 (10-32) 进行建模。使用 1.6:1 的比例既保留了这些观察结果的基本特征，又提供了更接近 LoG 函数的“工程”近似( Marr 和 Hildreth [1980])。为了使 LoG 和 DoG 具有相同的过零点，必须根据以下公式选择 LoG 的 σ 值(参见问题 10.19)：

(10-33)

$\displaystyle \sigma^{2}=\frac{\sigma_{1}^{2}\sigma_{2}^{2}}{\sigma_{1}^{2}-\sigma_{2}^{2}}\ln \left [\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}\right ]$

虽然当 σ 取此值时，LoG 和 DoG 的零交叉点相同，但它们的振幅尺度不同。我们可以通过缩放这两个函数，使它们在原点处具有相同的值，从而使它们兼容。

图 10.23(a) 和 (b) 中的轮廓线分别以 1:1.75 和 1:1.6 的标准偏差比率生成(按照惯例，所示曲线已反转，如图 10.21 所示)。LoG 轮廓线为实线，DoG 轮廓线为虚线。所示曲线是穿过 LoG 和 DoG 阵列中心的强度轮廓线，分别通过对公式 (10-29) 和 (10-32) 进行采样生成。所有曲线在原点的振幅均已归一化为 1。如图 10.23(b) 所示，1:1.6 的比率对 LoG 和 DoG 函数的近似效果略好(例如，比较两幅图的底部波瓣)。

----------------------图 10.23：(a) 使用 σ 比为 1.75:1 的 LoG (实线) 和 DoG(虚线)轮廓的负片。(b) 使用 σ 比为 1.6:1 的轮廓。---------------------------------

Gauss核是可分离的(参见第 3.4 节)。因此，LoG 和 DoG 滤波操作都可以用一维卷积来实现，而无需直接使用二维卷积(参见习题 10.19)。对于大小为 M × N 的图像和大小为 n × n 的卷积核，这样做可以将每次卷积所需的乘法和加法次数从二维卷积的 $n^{2}MN$ 成正比减少到一维卷积的 nMN 成正比。这种实现方式的差异非常显著。例如，如果 n = 25，则一维实现所需的乘法和加法运算次数将比使用二维卷积少约 12 倍。

10.2.6.2 Canny边缘检测子(The Canny Edge Detector)

尽管算法更为复杂，但本节讨论的 Canny 边缘检测子( Canny [1986] )的性能总体上优于目前讨论过的其他边缘检测子。Canny 的方法基于三个基本目标：

(1) 低误关率。应该找到所有边，并且不应该有虚假响应。

(2) 边缘点应当具有良好局部性。定位到的边缘必须尽可能接近真实边缘。即，检测子标记为边缘的点与真实边缘中心之间的距离应尽可能小。

(3) 单边缘点响应。对于每一个真实的边缘点，检测子应该只返回一个点。也就是说，真实边缘周围的局部极大值数量应该尽可能少。这意味着检测器不应该在只有一个边缘点的情况下识别出多个边缘像素。

Canny 的工作精髓在于用数学方式表达上述三个准则，然后尝试求得这些公式的最优解。通常，求得满足所有上述目标的闭式(closed-form)解是困难的(甚至不可能的)。然而，使用数值优化方法处理受加性Gauss白噪声干扰的一维阶跃边缘(回想一下，白噪声是指在特定频带内频谱连续且均匀的噪声。Gauss白噪声是指振幅值分布呈高斯分布的白噪声。Gauss白噪声能够很好地近似许多现实世界的情况，并能生成易于数学处理的模型。它还有一个很有用的特性，那就是它的振幅值在统计上是独立的)，得出的结论是：Gauss 函数的一阶导数是最优阶跃边缘检测子的良好近似:

(10-34) $\displaystyle \frac{d}{dx}e^{\displaystyle -\frac{x^{2}}{2\sigma^{2}}}=\frac{-x}{\sigma^{2}}e^{\displaystyle -\frac{x^{2}}{2\sigma^{2}}}$

近似值仅比使用优化的数值解差约 20% (这种程度的差异在大多数应用中通常是视觉上无法察觉的)。

将上述结果推广到二维需要认识到，一维方法仍然适用于边缘法线方向(见图 10.12)。由于法线方向事先未知，因此需要将一维边缘检测器应用于所有可能的方向。这项任务可以通过以下方式近似完成：首先使用二维圆形 Gauss 函数平滑图像，计算结果的梯度，然后使用梯度的大小和方向来估计每个点的边缘强度和方向。

令 f (x , y ) 表示输入图像，G (x , y ) 表示 Gauss 函数：

(10-35) $\displaystyle G (x ,y) = e^{\displaystyle -\frac{x^{2}+y^{2}}{2\sigma^{2}}}$

将 f 与 G 进行卷积运算，我们就得到了一幅平滑图像 $f_{s}(x,y)$ ：

(10-36) $f_{s}(x ,y) = G (x ,y) {\bigstar} f (x ,y)$

接下来，需要计算梯度的大小和方向(角度)，如前所述：

(10-37) $M_{s}(x ,y) = {\Vert}{\nabla}f_{s} (x ,y){\Vert} = \sqrt{g_{x}^{2} (x,y)+g_{y}^{2} (x,y)}$

和

(10-38) $\displaystyle \alpha(x,y)=\tan^{-1}\left [\frac{g_{y}(x,y)}{g_{x}(x,y)} \right ]$

其中， $g_{x}(x ,y) = {\partial}f_{s}(x ,y)/{\partial}x$ ， $g_{y}(x ,y) = {\partial}f_{s}(x ,y)/{\partial}y$ 。图 10.14 中任意导数滤波器核对(pairs)可以用于求得 $g_{x}(x ,y)$ 和 $g_{y}(x ,y)$ 。公式 (10-36) 用一个其大小如下描述的 n × n Gauss 核来实现。谨记， ${\Vert}{\nabla}f_{s} (x ,y){\Vert}$ 和 $\alpha(x,y)$ 数组的大小与计算它们所依据的图像的大小相同。

梯度图像 ${\Vert}{\nabla}f_{s} (x ,y){\Vert}$ 典型的是在局部极大值周围存在宽阔的脊状结构。下一步是使这些脊状结构变薄。一种方法是使用非极大值抑制。该方法的本质是指定边缘法线(梯度向量)的若干离散方向。例如，在一个 3 × 3 的区域中，我们可以为穿过该区域中心点的边缘定义四个方向(注：每条边都有两种可能的方向。例如，法线方向为 0° 的边和法线方向为 180° 的边是同一条水平边 )：水平、垂直、+45° 和 −45°。图 10.24(a) 显示了水平边缘的两种可能方向的情况。由于我们需要将所有可能的边缘方向量化为四个范围，因此我们必须定义一个方向范围，在该范围内，我们将边缘视为水平边缘。我们根据边缘法线的方向来确定边缘方向，而边缘法线可以直接使用公式 (10-38) 从图像数据中获得。如图 10.24(b) 所示，如果边缘法线的方向在 -22° 到 22° 或 -157° 到 157° 的范围内，我们就称该边缘为水平边缘。图 10.24(c) 显示了所考虑的四个方向对应的角度范围。

------------------图 10.24：(a) 3×3 邻域内水平边的两种可能方向（阴影部分）。(b) 水平边的法线方向角 a 的取值范围(阴影部分)。(c) 3×3 邻域内四种边方向的法线角度范围。每一种边方向都有两个范围，分别用相应的阴影部分表示。---------------------

令 $d_{1}$ ， $d_{2}$ , $d_{3}$ 和 $d_{4}$ 表示刚才讨论的一个 3×3 区域的四个基本边缘方向：分别是水平方向 -45°，垂直方向 +45°。我们可以针对以 α 中任意点 ( x ,y ) 为中心的 3 × 3 区域制定如下非极大值抑制方案：

(1) 求得最接近 α ( x ,y ) 的方向 $d_{k}$ 。

(2) 令 K 表示 ${\Vert}{\nabla} f_{s}{\Vert}$ 在 ( x ,y ) 处的值。若在点( x ,y ) 沿 $d_{k}$ 方向的一个或全部两个领域处 K 的值小于 ${\Vert}{\nabla} f_{s}{\Vert}$ 的值，则令 $g_{N}(x,y)=0$ (抑制)；否则，令 $g_{N}(x,y)=K$ 。

我们重复 x 和 y 的所有值，执行这个过程就获得了一幅与 $f_{s}(x,y)$ 大小相同的非极大化抑制图像 $g_{N}(x,y)$ 。例如，针对图 10.24(a) ，令 ( x ,y ) 位于 $p_{5}$ 处，并假设水平边缘穿过 $p_{5}$ ，在第 (2) 中的感兴趣像素应当是 $p_{2}$ 和 $p_{8}$ 。图像 $g_{N}(x,y)$ 仅包含薄边缘；其等于其非极大边缘点受抑制的图像 ${\Vert}{\nabla} f_{s}(x,y){\Vert}$ 。

最后一步操作是对 $g_{N}(x,y)$ 进行阈值处理，以减少假边缘点。在 Marr-Hildreth 算法中，我们使用单个阈值来实现这一点，其中所有低于阈值的值都设置为 0。如果阈值设置得太低，仍然会存在一些假边缘(称为假正值(false positives))。如果阈值设置得太高，则会消除有效的边缘点(假负值(false negatives))。Canny 算法试图通过使用滞后阈值来改进这种情况，正如我们将在 10.3 节中讨论的那样，它使用两个阈值：低阈值 $T_{L}$ 和高阈值 $T_{H}$ 。实验证据( Canny [1986])表明，高阈值与低阈值的比值应在 2:1 到 3:1 的范围内。

我们可以将阈值分割操作可视化为创建两个额外的图像：

(10-39) $g_{NH}( x ,y ) = g_{N}( x ,y ) \ge T_{H }$

和

(10-40) $g_{NL}( x ,y ) = g_{N}( x ,y ) \ge T_{L}$

初始时， $g_{NH}( x ,y )$ 和 $g_{NL}( x ,y )$ 设为 0 。在阈值处理之后， $g_{NH}( x ,y )$ 通常将比 $g_{NL}( x ,y )$ 具有更小的非零像素，但 $g_{NH}( x ,y )$ 中的所有非零像素将包含于 $g_{NL}( x ,y )$ 中，因后面的这幅图像由一个更低的阈值构成。我们通过令

(10-41) $g_{NL}( x ,y ) = g_{NL}( x ,y ) - g_{NH}( x ,y )$

从而从 $g_{NL}( x ,y )$ 中消除来自 $g_{NH}( x ,y )$ 的所有非零元素。 $g_{NH}( x ,y )$ 和 $g_{NL}( x ,y )$ 中的非零像素可以分别视为“强”边缘像素和“弱”边缘像素。在经过阈值运算之后， $g_{NH}( x ,y )$ 中的所有强像素都假设为有效的边缘像素，因此立即标识。根据 $T_{H}$ 的值， $g_{NH}( x ,y )$ 中的边缘通常有间隙。较长的边缘使用下述过程形成：

(a) 在 $g_{NH}( x ,y )$ 中定位下一个未遍历的边缘像素 p 。

(b) 将 $g_{NL}( x ,y )$ 中使用(比如)8连通连接到 p 的所有弱像素标为有效边缘像素。

(c) 若 $g_{NH}( x ,y )$ 中所有非零像素都已遍历，则转向 (d) ，否则回到 (a) 。

(d) 将 $g_{NL}( x ,y )$ 中未标为有效边缘像素的所有像素置为 0 。

在此过程的最后，Canny 算法的最后图像输出是通过将来自 $g_{NL}( x ,y )$ 中的所有像素追加到 $g_{NH}( x ,y )$ 而形成的。

我们用两个附加图像 $g_{NH}( x ,y )$ 和 $g_{NL}( x ,y )$ 来简化这个讨论。在实践中，滞后阈值可以直接在非极大值抑制期间实现，阈值可以直接在 $g_{N}( x ,y )$ 上实现，方法是形成强像素列表和与其相连的弱像素列表。

总之，Canny 边缘检测算法由下列步骤构成：

(1) 用一个 Gauss 滤波器平滑输入图像。

(2) 计算梯度幅值和角度大小(方向)。

(3) 对梯度幅值图像应用非极大值抑制。

(4) 使用双阈值法和连通性分析来检测和连接边缘。

虽然非极大值抑制后的边缘比原梯度边缘更细，但前者仍然可能超过一个像素厚。为了获得一个像素厚的边缘，通常在步骤 4 之后再执行一次边缘细化算法(参见第 9.5 节)。

如前所述，平滑是通过将输入图像与大小为 n × n 的Gauss核进行卷积来实现的，该Gauss核的大小必须选定。一旦确定了 σ 值，我们就可以使用与 Marr-Hildreth 算法相关的讨论方法来确定一个奇数 n 值，该奇数 n 值能够使Gauss滤波器在指定的 σ 值下发挥“完全”的平滑能力。

关于实现的最后几点说明：如前文在讨论 Marr-Hildreth 边缘检测器时所述，公式 (10-35) 中的二维Gauss函数可以分解为两个一维Gauss函数的乘积。因此，Canny 算法的第一步可以表述为对图像的行(列)逐一进行一维卷积，然后再对结果的列(行)进行处理。此外，如果我们使用公式 (10-19) 和 (10-20) 中的近似值，也可以将第二步所需的梯度计算实现为一维卷积(参见习题 10.22)。

例 10.8：Canny边缘检测方法的例释和比较。

图 10.25(a) 显示了常见的建筑物图像。为了进行比较，图 10.25(b) 和 (c) 分别显示了使用阈值梯度法得到的图 10.20(b) 的结果，以及使用 Marr-Hildreth 检测子得到的图 10.22(d) 的结果。需要注意的是，生成这两幅图像时使用的参数旨在检测主要边缘，同时尽量减少“无关”特征，例如砖块和屋顶瓦片的边缘。

图 10.25(d) 展示了使用 Canny 算法得到的结果，参数为 $T_{L} = 0.04$ , $T_{H} = 0.10$ (低阈值的 2.5 倍) , σ = 4 以及大小为 25 × 25 的卷积核，该卷积核对应于不小于 6σ 的最小奇数。这些参数是通过实验选择的，旨在实现前一段中针对梯度图像和 Marr-Hildreth 图像所述的目标。将 Canny 图像与其他两幅图像进行比较，我们可以看到 Canny 算法在主要边缘的细节方面有显著改进，同时对无关特征的剔除也更加有效。例如，请注意，Canny 算法检测到了图像上部砖块周围混凝土带的两条边缘，而阈值梯度图像则丢失了这两条边缘，Marr-Hildreth 方法仅检测到了上方的边缘。在滤除无关细节方面，Canny 图像中没有出现任何由屋顶瓦片引起的边缘；但在另外两幅图像中并非如此。Canny 检测图像中的线条在连续性、纤细度和直线度方面也表现更佳。正是由于这些优势，Canny 算法才成为边缘检测的首选工具。

---------------------------图 10.25：(a) 原图像，尺寸为 834 × 1114 像素，强度值缩放至 [0, 1] 范围。(b) 平滑图像的阈值梯度。(c) 使用 Marr-Hildreth 算法获得的图像。(d) 使用 Canny 算法获得的图像。请注意，与其他两种算法相比，Canny 算法获得的图像质量有显著提升。--------------------------------------------------

例 10.9：本节讨论的三种主要边缘检测方法的另一个示例。

为了进一步比较本节讨论的三种主要边缘检测方法，请看图 10.26(a)，该图显示了一张 512 × 512 的头部 CT 图像。我们的目标是提取大脑外轮廓(图像中的灰色区域)、脊柱区域轮廓(位于鼻后方，朝向大脑前方)以及头部外轮廓的边缘。我们希望生成尽可能细且连续的轮廓，同时消除与眼睛和大脑区域灰质相关的边缘细节。

图 10.26(b) 显示了经过阈值处理的梯度图像，该图像首先使用 5 × 5 平均核进行平滑处理。达到所示结果所需的阈值为梯度图像最大值的 15%。图 10.26(c) 显示了使用 Marr-Hildreth 边缘检测算法获得的结果，阈值为 0.002，σ = 3，核大小为 19 × 19。图 10.26(d) 是使用 Canny 算法获得的， $T_L = 0.05$ , $T_{H} = 0.15$ (低阈值的 3 倍)，σ = 2，核大小为 13 × 13。

就边缘质量和去除无关细节的能力而言，图 10.26 的结果与前一个例子的结果和结论非常吻合。还需注意的是，Canny 算法是唯一能够为大脑后边界和脊髓最近边界生成完全连续边缘的算法。它也是唯一能够在去除原图像中所有与灰质相关的边缘的同时，找到最清晰轮廓的算法。

----------------------------图 10.36：(a) 尺寸为 512 × 512 像素的头部 CT 图像，强度值缩放至 [0, 1] 范围。(b) 平滑图像的阈值梯度。(c)使用 Marr-Hildreth 算法获得的图像。(d)使用 Canny 算法获得的图像。(原图像由Vanderbilt大学的 David R. Pickens 博士提供。) ------------------------------------------------------

Canny算法性能提升的代价是其实现方式比前面讨论的两种方法复杂得多。在某些应用中，例如实时工业图像处理，成本和速度要求通常决定了使用更简单的技术(主要是阈值梯度法)。当边缘质量是主要考虑因素时，Marr-Hildreth算法和Canny算法，尤其是后者，提供了更优的选择。

10.2.7 连接边缘点(Linking Edge Points)

在理想情况下，边缘检测应该只生成位于边缘上的像素集合。然而，在实际应用中，由于噪声、非均匀光照导致的边缘断裂以及其他会引入强度值不连续性的因素，这些像素很少能完整地表征边缘。因此，边缘检测之后通常会进行连接算法，旨在将边缘像素组装成有意义的边缘和/或区域边界。本节将讨论两种基本的边缘连接方法，这两种方法代表了实际应用中常用的技术。第一种方法需要了解局部区域( 例如，3×3 邻域 )中的边缘点信息，而第二种方法是一种全局方法，它处理整个边缘图。事实上，连接区域边界上的点也是下一章讨论的一些分割方法以及从分割图像中提取特征(我们将在第 11 章中讨论)的重要组成部分。因此，在接下来的两章中，你将会遇到更多边缘点连接方法。

10.2.7.1 局部处理法(Local Processing)

一种连接边缘点的简单方法是分析每一个被前几节讨论的技术判定为边缘点的点 ( x , y ) 周围小邻域内的像素特征。所有根据预定义标准相似的点都会被连接起来，形成一条由符合指定标准的共同属性像素组成的边缘。

在这种局部分析中，用于建立边缘像素相似性的两个主要属性是：(1) 梯度强度(幅值)和 (2) 梯度方向。第一个属性基于公式 (10-17)。令 $S_{xy}$ 表示图像中以点 ( x , y ) 为中心的邻域的坐标集合。对于 $S_{xy}$ 中一个坐标为 ( s , t ) 的边缘像素，若

(10-42) $|M( s , t ) - M( x , y )| {\le} E$

其中， E 是正阈值，则称此像素与点 (x , y) 处的像素在幅值上相似。

梯度向量方向的角度根据公式 (10-18) 计算而得。对于 $S_{xy}$ 中一个坐标为 ( s , t ) 的边缘像素，若

(10-43) $|\alpha( s , t ) - \alpha( x , y )| \le A$

其中，A 是一个正角度阈值，则称此像素与点 (x , y) 处的像素在角度上相似。如前所述，(x , y) 处的边缘方向与该点的梯度向量方向垂直。

对于 $S_{xy}$ 中一个坐标为 ( s , t ) 的像素，若同时满足幅度和方向的标准，则视为连与点(x , y) 相连。对每一个边缘像素重复此过程。随着邻域中心从一个像素移动到另一个像素，就会记录所有连接的点。一种简单的记录方法是为每组连接的边缘像素分配不同的强度值。

上述公式计算量很大，因为需要检查每一个点的所有邻域。一个特别适合实时应用的简化方法包括以下步骤：

(1) 计算输入图像 f ( x , y ) 梯度幅值 M( x , y ) 和角度数组 α ( x , y ) 。

(2) 生成一幅二值图像 g ( x , y )，其任意点 ( x , y ) 处的值计算公式为

$g(x,y)=\left \{ \begin{array}{l} 1 \hspace{0.2cm}(\text{if\;}M( x ,y )> T_M \text{\;AND\;} \alpha( x ,y )=A{\pm}T_{A}) \\ \\ 0 \hspace{0.2cm}(\text{otherwise}) \end{array} \right.$

其中， $T_{M}$ 是一个阈值，A 是一个具体的角度方向，且 ${\pm}T_{A}$ 定义了一个关于 A 的可接受方向“带(band)”。

(3) 扫描 g 的每一行，并将每一行中长度不超过指定长度 L 的所有间隙( 0 的集合)填充 (设置为 1 )。注意，根据定义，间隙的两端由一个或多个 1 限定。每一行单独处理，彼此之间没有“记忆”。

(4) 要检测其他方向的间隙 θ ，将 g 旋转该角度，并在步骤 (3) 中应用水平扫描程序。将结果旋转 -θ 。

当关注水平和垂直边缘连接时，步骤 (4) 就变成了一个简单的过程：将图像 g 旋转 90 度，扫描各行，然后将结果旋转回去。这是实践中最常见的应用，正如以下示例所示，这种方法可以产生良好的结果。通常，图像旋转是一个计算量很大的过程，因此，当需要沿多个角度方向进行连接时，将步骤 (3) 和 (4) 合并为一个径向扫描过程更为实用。

例 10.10：使用局部处理实现的边缘链接。

图 10.27(a) 显示了一张 534 × 566 像素的车辆后部图像。本例旨在说明如何使用前述算法寻找尺寸合适的矩形，使其适合作为车牌的候选区域。这些矩形的形成可以通过检测明显的水平和垂直边缘来实现。图 10.27(b) 显示了梯度幅值图像 M( x , y )。图 10.27(c) 和 (d) 显示了算法步骤 (3) 和步骤 (4) 的结果，这些结果通过以下设置获得： $T_{A}$ 等于最大梯度值的 30%，A = 90°， $T_{A}=45^{\circ}$ ，并填充所有小于等于 25 像素的间隙(约占图像宽度的 5% ) 。为了检测车牌框的圆角以及车辆的后窗，需要较大的角度范围。图 10.27(e) 是前两幅图像进行逻辑或运算的结果，图 10.27( f ) 则是通过对图 10.27(e) 进行细化处理(细化方法详见 9.5 节)得到的。如图 10.27( f ) 所示，图像中清晰地检测到了车牌对应的矩形区域。利用车牌宽高比具有特定比例( 例如，美国车牌的宽高比为 2:1 ) 这一特性，可以轻松地将车牌从图像中的所有矩形区域中分离出来。

------------------------图 10.27：(a) 车辆后部图像。(b) 梯度幅度图像。(c) 水平连接的边缘像素。(d) 垂直连接的边缘像素。(e) (c)和(d)的逻辑或运算结果。( f )使用形态学细化后的最终结果。(原始图像由Perceptics公司提供。)----------------------

10.2.7.2 籍Hough变换的全局处理法(Global Processing Using the Hough Transform)

上一节讨论的方法适用于已知像素所属对象信息的情况。然而，我们常常需要在非结构化环境中工作，此时我们只有边缘图，而对目标对象的位置一无所知。在这种情况下，所有像素都可能成为连接对象，因此必须根据预定义的全局属性来决定是否接受或排除。本节将介绍一种基于像素集合是否位于特定形状曲线上的方法。一旦检测到这些曲线，它们就构成了目标对象的边缘或区域边界。

给定图像中的 n 个点，假设我们要求得这些点的子集，这些点位于直线上。一种可能的解决方案是求得由每对点确定的所有直线，然后求得所有接近特定直线的点的子集。这种方法需要求得 $n(n - 1) {\;\sim\;} n^{2}$ 条直线，然后对每个点与所有直线进行 $(n)( n(n - 1)) {\;\sim \;} n^{3}$ 次比较。在大多数应用中，这是一项计算量巨大的任务。

Hough [1962] 提出了一种折衷的方案，通常称为 Hough 变换。令 $( x_{i}, y_{i} )$ 为 xy 平面上的一点，并考虑一条直线的斜截式公式： $y_{i} = ax_{i} + b$ 。有无限多条直线经过点 $( x_{i}, y_{i} )$ ，但对于各种取值的 a 和 b ，所有这些直线都满足公式： $y_{i} = ax_{i} + b$ 。然而，将这个公式写成 $b = - ax_{i} + y_{i }$ 并考虑 ab 平面(又称参数空间)便得到一个定点 $( x_{i}, y_{i} )$ 的单一直线公式。此外，第二点 $( x_{j}, y_{j} )$ 也有一条与其参数空间关联的单一直线，其与和 $( x_{i}, y_{i} )$ 关联的直线相交于参数空间中的某一点 $( a^{'} , b^{'} )$ ，其中， $a^{'}$ 和 $b^{'}$ 分别是 xy 平面中同时包含 $( x_{i}, y_{i} )$ 和 $( x_{j}, y_{j} )$ 直线的斜率和截距( 当然，我们假设直线不是平等的 )。实事上，这条直线上的所有点在参数空间中都有直线与其相交于点 $( a^{'} , b^{'} )$ 。图 10.28 图释了这些概念。

--------------------图 10.28：(a) xy 平面。(b)参数空间。---------------

在原则上，可以绘制出对应于 xy 平面上所有点 $( x_{k} , y_{k} )$ 的参数空间直线，并通过识别参数空间中大量参数空间直线相交的点来求得该平面上的主要直线。然而，这种方法的难点在于，当直线趋近于垂直方向时，斜率a趋近于无穷大。解决此难点的一种方法是使用直线的法线表示：

(10-44) $x\cos{\theta}+y\sin{\theta} = \rho$

图 10.29 (a) 说明了参数 ρ 和 θ 的几何解释。一条水平直线有 θ = 0°，且 ρ 等于正的 y 截距。类似地，一条垂直直线有 θ = 90°，且 ρ 等于正的 y 截距，或 θ = -90°，ρ 等于负的 y 截距 (我们将角度限定于 -90°≤ θ ≤ 90° )。图 10.29(b) 中的每一条正弦曲线都代表穿过 xy平面上的一个特定的点 $( x_{k} , y_{k} )$ 的直线族。图 10.29(b) 中的交点 $(\rho^{'} ,\theta^{'})$ 对应经过图 10.29(a) 中点 $( x_{i}, y_{i} )$ 和 $( x_{j}, y_{j} )$ 的直线。

Hough变换的计算优势(attractiveness)源于将 ρθ 参数空间细分为所谓的累加器单元( accumulator unit )，如图 10.29(c) 所示，其中， $( \rho_\text{min} , \rho_\text{max} )$ 和 $( \theta_\text{min} , \theta_\text{max} )$ 是参数值的期望值范围：-90°≤ θ ≤ 90° 且 -D ≤ ρ ≤ D ，其中，D 是一幅图像中对角之间的最大距离。坐标 ( i , j ) 处具有累加器值 A( i , j ) 的单位格对应与参数空间坐标 $( \rho_{i} , {\theta}_{j} )$ 相关的正方形。开始时，这些单元格置为零。则对于xy平面上的每一个非背景点 $( x_{k} , y_{k} )$ ，我们令 θ 等于每一个 θ 轴上允许子划分的值，并使用方程 $\rho = x_{k}\cos\theta + y_{k}\sin\theta$ 来求解对应的 ρ 。所求得的 ρ 值四舍五入为沿 ρ 轴的最接近的允许单元格。若所选择的 $\theta_{q}$ 导致了解 $\rho_{p}$ , ，则我们令 A( p , q ) = A( p , q ) + 1 。在这个过程的最后，一个单元格 A( i , j ) 中的一个 K 值意味着 xy平面上有 K 个点位于直线 $x\cos{\theta_{j}}+y\sin{\theta_{j}}=\rho_{j}$ 上。ρθ 平面上的细分数量决定了这些点的共线精度。可以证明(参见问题 10.27)，刚才讨论的方法中的计算次数与 n ( xy 平面上非背景点的数量)呈线性关系。

----------------------图 10.29：(a) xy平面上一条直线的 ( ρ , θ ) 参数化。(b) ρθ 平面上的正弦曲线；与经过 xy 平面上点 $( x_{i} , y_{i} )$ 和 $( x_{j} , y_{j} )$ 的直线对应的交点 $(\rho^{'},\theta^{'})$ 。将 ρθ 平面划分成累加器单位。--------------------------------

例 10.11：Hough变换的一些基本属性。

图 10.30 展示了基于公式 (10-44) 的Hough变换。图 10.30(a) 显示了一幅大小为 M × M ( M = 101 ) 的图像，其中包含五个标记的白点；图 10.30(b) 显示了将这五个点分别映射到 ρθ 平面上的结果，其中 ρ 轴和 θ 轴的划分单位均为 1 。θ 轴的取值范围为 ±90°，ρ 轴的取值范围为 $\pm\sqrt{2}M$ 。如图 10.30(b) 所示，每一条曲线都呈现出不同的正弦形状。由点 1 的映射得到的水平线是一条零振幅的正弦曲线。

图 10.30(b) 中标记为 A(不要与累加器值混淆)和 B 的点说明了Hough变换的共线性( collinearity )检测特性。例如，点 B 标记了对应于 xy 图像平面中点 2 , 3 和 4 的曲线的交点。点 A 的位置表明这三个点位于一条穿过原点 ( ρ = 0 ) 且方向为 −45° 的直线上 [参见图 10.29(a)]。类似地，参数空间中相交于点 B 的曲线表明点 2, 3 和 4 位于一条方向为 45° 的直线上，该直线到原点的距离为 ρ = 71 (图像原点到对角线距离的一半，四舍五入到最接近的整数)。最后，图 10.30(b) 中标记为 Q ,R 和 S 的点表明，Hough变换在参数空间的左右边界处呈现出反射邻接关系。这一特性是由于 ρ 和 θ 在 ±90° 边界处改变符号的方式所致。

--------------------图 10.30：(a) 尺寸为 101 × 101 像素的图像，包含五个白点(四个位于角上，一个位于中心)。(b)对应的参数空间。-------------------------------

尽管目前为止的研究重点都集中在直线上，但Hough变换适用于任何形如 g(v, c) = 0 的函数，其中 v 是坐标向量，c 是系数向量。例如，位于圆上的点

(10-45) $(x-c_{1} )^{2} + (x-c_{2} )^{2} = c_{3}^{2}$

可以使用前述讨论的方法进行检测。区别在于引入了三个参数 $c_{1}$ , $c_{2}$ , 和 $c_{3}$ ，从而形成一个具有立方体状单元的三维参数空间，以及形如 A ( i, j, k ) 的累加器。具体步骤是：递增 $c_{1}$ 和 $c_{2}$ ，求解满足公式 (10-45) 的 $c_{3}$ 值，并更新与三元组 $(c_{1},c_{2},c_{3})$ 关联的累加器单元。显然，Hough变换的复杂度取决于给定函数表示中的坐标和系数的数量。如前所述，Hough 变换可以推广到检测没有简单解析表示的曲线，也可以应用于灰度图像。

回到边缘连接问题，基于Hough变换的方法如下：

(1) 使用本节前面讨论的任意方法获取二值边缘图。

(2) 指定 ρθ 平面上的细分。

(3) 检查累加器单元中像素密度较高的区域。

(4) 检查所选单元中像素之间的关系(主要检查连续性)。

在这种情况下，连续性通常基于计算与给定累加器单元对应的不连续像素之间的距离。如果间隙长度小于指定的阈值，则与给定单元关联的线段中的间隙将被连接起来。能够根据方向对线段进行分组是一个全局概念，适用于整个图像，只需要检查与特定累加器单元关联的像素即可。以下示例说明了这些概念。

例 10.12：使用Hough变换实现边缘连接。

图 10.31(a) 显示了一张机场的航拍图像。本例的目标是使用Hough变换提取定义主跑道的两条边缘。此类问题的解决方案可能对自主空中导航等应用具有重要意义。

第一步是获取边缘图。图 10.31(b) 显示了使用 Canny 算法获得的边缘图，其参数和步骤与例 10.9 相同。为了计算Hough变换，可以使用前面讨论过的任何其他边缘检测技术获得类似的结果。图 10.31(c) 显示了使用 θ 方向 1°增量和 ρ 方向 1 像素增量获得的Hough参数空间。

目标跑道与正北方向的夹角约为 1°，因此我们选择对应于 ±90° 且计数最高的单元格，因为跑道是沿这些方向的最长线段。图 10.31(c) 边缘的小方框突出显示了这些单元格。如前所述，结合图 10.30(b)，Hough变换在边缘处表现出邻接性。另一种解释方式是，方向为 +90° 的直线和方向为 −90° 的直线是等价的(即，它们都是垂直的)。图 10.31(d) 显示了与刚才讨论的两个累加器单元格对应的直线，图 10.31(e) 显示了叠加在原图像上的直线。这些直线是通过连接所有不超过图像高度 20%(约 100 像素)的间隙而得到的。这些直线清晰地对应于我们感兴趣的跑道的边缘。

请注意，解决此问题所需的唯一信息是跑道的方向以及观察者相对于跑道的位置。换句话说，自主导航的车辆会知道，如果目标跑道朝北，且车辆的行驶方向也朝北，则跑道在图像中应呈垂直方向。其他相对方向的处理方式类似。世界各地跑道的方向信息可在飞行图表中找到，而行驶方向则可利用 GPS(全球定位系统)信息轻松获取。这些信息还可用于计算车辆与跑道之间的距离，从而估算诸如相对于图像大小的预期线段长度等参数，正如我们在本例中所做的那样。

-------------图 10.31：(a) 一张 502 × 564 像素的机场航拍图像。(b) 使用 Canny 算法获得的边缘图。(c) 霍夫变换参数空间(方框突出显示与长垂直线相关的点)。(d) 图像平面上与方框突出显示的点对应的直线。(e) 叠加在原始图像上的直线。------------