第3章〓典型分割算法 教学视频 图像分割技术的研究已有60多年的历史，而且一直都在不断发展，已提出的各种图像分割方法非常多，10年前的统计就表明有关图像分割文献的数量已达近10万［章2014c］。本章仅介绍几种具有比较特殊思路的典型算法作为例子。在2．1节中曾将图像分割技术分为4大类，本章各节介绍的内容也分别对应这4大类分割技术。 根据上述讨论，本章各节将安排如下。 3．1节先介绍利用二阶导数对兴趣点的检测(包括计算拉普拉斯值和海森值)，再讨论既可用于边缘检测也可用于角点检测的最小核同值区算子的具体方法和特点，最后介绍哈里斯兴趣点算子，它们都属于并行边界类技术。 3．2节介绍一种结合了图像整体信息和局部信息，并利用了图结构的特殊方法——图割方法，它属于串行边界类技术。 3．3节介绍3种不同的属于并行区域类技术的方法，分别是借助小波变换的多分辨率特性阈值化方法、借助过渡区的阈值化方法和借助均移确定空间聚类的方法。 3．4节介绍一种结合了图像整体信息和局部信息的特殊方法——基于分水岭的方法，除了基本原理，还讨论了对基本算法的改进和扩展，以及基于深度学习的方法。虽然分水岭本身对应边界，但这种方法更多的是利用边界所包围的区域性质来工作的，所以是一种串行区域类技术。 3．1兴趣点检测 兴趣点或感兴趣点泛指图像中或目标上具有特定几何性质或属性性质的点，例如，角点、拐点、梯度极值点等。对这些点的检测多可并行进行，检测出这些点对获得目标边界常常能起很大作用，所以可看作基于边缘的并行分割基础。对兴趣点的检测有很多方法，下面先讨论利用二阶导数来检测角点，再介绍两种比较有特色的检测算子，可检测多种类型的兴趣点。 3．1．1二阶导数检测角点 若一个像素在其小邻域中具有两个明显不同的边缘方向，则常可被看作一个角点(或者说小邻域中有两个边缘段，其朝向更偏于互相垂直)，也有人将局部曲率(离散曲率计算可参见11．5．2小节)较大的边缘点看作角点。典型的角点检测器仍多基于像素灰度差(梯度)。 利用二阶导数算子检测角点的原理与利用一阶导数算子检测边缘有些类似。在2D图像中，由各个二阶导数组成的对称矩阵(也称为局部结构矩阵)可写成(下标指示偏导数的方向) M(2)=MxxMxyMyxMyy,Mxy=Myx(3．1．1) 它给出在被检测点(可设为坐标原点)的局部曲率的信息。上述矩阵有两个特征值，可考虑它们的3种组合情况： (1) 两个特征值都很小，表示被检测点局部的邻域平坦，被检测点非边缘点或角点； (2) 两个特征值一个很小而另一个很大，表示被检测点局部的邻域呈脊线状，沿一个方向平坦而沿另一个方向变化迅速，被检测点为边缘点； (3) 两个特征值都很大，表示被检测点为角点。 这种检测方法对平移和旋转变换不敏感，也对光照和视角变化有一定的稳定性。 如果旋转坐标系统，可将M(2)变换成对角形式： M～(2)=Mx～x～00My～y～=K100K2(3．1．2) 此时的二阶导数矩阵给出了在原点的主曲率，即K1和K2。 再回到如M(2)的矩阵，其秩和行列式都是旋转不变的。进一步可分别得到拉普拉斯值和海森值： Laplacian=Mxx+Myy=K1+K2(3．1．3) Hessian=det(M(2))=MxxMyy-M2xy=K1K2(3．1．4) 基于一个函数的二阶偏导数构成的方阵，拉普拉斯算子和海森算子分别计算拉普拉斯值和海森值。前者对边缘和直线都给出较强的响应，所以不太适合用来检测角点。后者描述了函数的局部曲率，虽然对边缘和直线没有响应，但是在角点邻域中有较强的响应，所以比较适合用来检测角点。不过海森算子在角点自身处响应为零而且在角点周边的符号是不一样的，所以需要较复杂的分析过程以确定角点的存在性并准确地对角点定位。为避免这个复杂的分析过程，可先计算曲率(K)与局部灰度梯度(g)的乘积： C=Kg=KM2x+M2y=MxxM2y-2MxyMxMy+MyyM2xM2x+M2y(3．1．5) 再沿边缘法线方向利用非最大消除来确定角点的位置。 3．1．2最小核同值区算子 最小核同值区(SUSAN)算子是一种很有特色的检测算子［Smith 1997］，它既可以检测边缘点，也可以检测角点。 1． 核同值区 先借助图3．1．1来解释检测的原理。这里设有一幅图像，取其顶部一个矩形区域，其上部为亮区域，下部为暗区域，分别代表背景和目标。现在考虑有一个圆形的模板(其大小由模板边界所限定)，其中心称为“核”，用“+”表示。图3．1．1中给出了该模板放在图像中6个不同位置的示意情况，从左边数过去，第1个模板全部在亮区域，第2个模板大部在亮区域，第3个模板一半在亮区域(另一半在暗区域)，第4个模板大部在暗区域，第5个模板全部在暗区域，第6个模板的1/4在暗区域。 图3．1．1圆形模板在图像中的不同位置 如果将模板中各个像素的灰度都与模板中心核像素的灰度进行比较，就会发现总有一部分模板像素的灰度与核像素的灰度相同或相似。这部分区域可称为核同值区(USAN)，即与核像素灰度具有相同数值的区域。核同值区包含了很多与图像结构有关的信息。利用这种区域的尺寸、重心、二阶矩等各种特征可以检测图像中的边缘和角点。从图3．1．1可见，当核像素处在图像中的灰度一致区域时，核同值区的面积会达到最大，第1个模板和第5个模板就属于这种情况。在第2个模板和第4个模板的情况下，核同值区的面积超过一半。这个面积当核处在直边缘处约为最大值的一半，第3个模板就属于这种情况。这个面积当核处在角点位置时更小，约为最大值的1/4，第6个模板就属于这种情况。 利用核同值区的面积具有上述变化的性质可检测边缘或角点。具体来说，核同值区面积较大(超过一半)时表明核像素处在图像中的灰度一致区域，在模板核接近边缘时减少为一半，而在接近角点时减少得更多，即在角点处取得最小值(约1/4)。如果将核同值区面积的倒数作为检测的输出，则可以通过计算输出的极大值来方便地确定角点的位置。 由上面的讨论可见，使用核同值区面积作为特征可起到增强边缘和角点的效果。基于这种模板的检测方式与其他常用的检测方式有许多不同之处，最明显的就是不需要计算微分，因而对噪声不太敏感。 2． 最小核同值区算子与角点检测 图3．1．237个像素的模板 借助核同值区的定义可讨论最小核同值区算子。最小核同值区算子采用圆形模板来得到各向同性的响应。在数字图像中，可用一个含有37个像素的模板来近似圆，如图3．1．2所示。这37个像素排成7行，分别有3、5、7、7、7、5、3个像素。这相当于一个半径约为3．4个像素的圆。如考虑到计算量，也有用简单的3×3模板来粗略近似圆形模板的。 设模板用N(x，y)表示，将其依次放在图像中每个点的位置。在每个位置，将模板内每个像素的灰度值与核的灰度值进行比较，计算其响应： C(x0,y0; x,y)=1，|f(x0,y0)-f(x,y)|≤T0，|f(x0,y0)-f(x,y)|>T(3．1．6) 图3．1．3函数C(·;·)示例 其中，(x0，y0)是核在图像中的位置坐标； (x，y)对应模板N(x，y)中除核以外的其他位置； f(x0，y0)和f(x，y)分别是在(x0，y0)和(x，y)处像素的灰度； T是一个灰度差的阈值； 函数C(·;·)代表输出响应。一般当检测到角点时，该函数类似一个门函数，如图3．1．3所示，设这里的阈值T取为10。 上述计算需要对模板中的每个像素进行，由此可得到一个输出的游程和： S(x0,y0)=∑(x,y)∈N(x0,y0)C(x0,y0; x,y)(3．1．7) 这个总和其实就是核同值区中的像素个数，或者说它给出了核同值区的面积。由前面的讨论可知，这个面积在角点处会达到最小。由式(3．1．6)和式(3．1．7)可知，阈值T既可用来帮助检测核同值区面积的最小值，也可以确定可消除的噪声的最大值。 实际应用最小核同值区算子时，需要将游程和S与一个固定的几何阈值G进行比较以做出判断。该阈值可设为3Smax/4，其中Smax是S所能取得的最大值(对37个像素的模板，最大值为36，所以3Smax/4=27)。初始的边缘响应R(x0，y0)根据下式得到： R(x0,y0)=G-S(x0,y0)，S(x0,y0)00，g(x,y)=0(3．3．8) 由此定义可知，在计算EAG时只用到了具有非零值梯度的像素。因为除去了零梯度像素的影响，所以称为“有效”梯度。EAG是图中非零梯度像素的平均梯度，它代表了图像中一个有选择的统计量。 进一步地，为了减少各种干扰的影响，定义以下特殊的剪切变换。它与一般剪切操作的不同之处是： 它把被剪切了的部分设成剪切值，从而避免了一般剪切操作在剪切边缘处产生大的反差而导致的不良影响。根据剪切部分的灰度值与全图灰度值的关系，这类剪切可分为低端剪切与高低端剪切两种。设L为剪切值，则低端剪切与高端剪切后的图可分别表示为 flow(x,y)=f(x,y)，f(x,y)>LL，f(x,y)≤L(3．3．9) fhigh(x,y)=L，f(x,y)≥Lf(x,y)，f(x,y)1(3．3．17) 给定dD空间的n个数据点xi，在点x算得的多变量密度核估计器为 f～h,k(x)=1nhd∑ni=1Kx-xih(3．3．18) 其中，h表示核尺寸，也称为核带宽。 由图3．3．5可见，需要确定fh,k(x)的梯度为零的点，即确定x以使fh,k(x)=0。均移方法就是不用估计概率密度函数而确定这些点位置的一种有效方法。这里把要解决的问题由估计密度变成估计密度梯度： f～h,k(x)=1nhd∑ni=1Kx-xih(3．3．19) 使用剖面为k(x)的核形式，并设它的导数k′(x)=-g(x)对所有x∈［0，∞］存在(除去有限个点)，则 Kx-xih=ckkx-xih2(3．3．20) 其中，ck是归一化常数。若K(x)=KE(x)，则剖面gE(x)是均匀/各向同性的； 若K(x)=KN(x)，则剖面gN(x)由与KN(x)相同的指数表达所定义。使用g(x)作为剖面的核G(x)=cgg(x2)，式(3．3．19)变成 f～h,k(x)=2cknh(d+2)∑ni=1(x-xi)k′x-xih2=2cknh(d+2)∑ni=1(x-xi)gx-xih2 =2cknh(d+2)∑ni=1gi∑ni=1xigi∑ni=1gi-x(3．3．21) 其中，∑ni=1gi设计为正，gi=g((x-xi)/h2)。 式(3．3．21)中的第一项正比于用核G算得的密度估计： f～h,G(x)=cgnh(d)∑ni=1gx-xih2(3．3．22) 第二项代表均移矢量mh,G(x)： mh,G(x)=∑ni=1xig((x-xi)/h2)∑ni=1g((x-xi)/h2)-x(3．3．23) 对核G依次确定的位置{yj}j=1,2,…为 yj+1=∑ni=1xig((yj-xi)/h2)∑ni=1g((yj-xi)/h2)-x(3．3．24) 其中，y1是核G的初始位置。 用核K计算得到的密度估计序列为 f～h,K(j)=f～h,K(yj)(3．3．25) 如果核K具有凸的和单减的剖面，则序列{yj}j=1,2,…和{f～h，K(j)}j=1,2,…都收敛，其中{f～h，K(j)}j=1,2,…是单增的，收敛速度依赖于所用的核。在离散数据上使用Epanechnikov核，在有限步迭代后就可以收敛。当考虑对数据加权时(如使用正态核)，均移过程无穷收敛。此时可用步长变化的下限值来停止计算。 均移方法的优缺点都与它对数据的全局表达有关。最主要的优点就是它的通用性。由于对噪声鲁棒，所以可用于各种实际场合。它可以处理任何聚类的形状和特征空间。它仅有的一个选择参数(核尺寸h)具有物理上可理解的意义。不过对h的选择也是对它应用的一个限制，因为确定一个合适的h并不是一件简单的事情。过大会导致最频值被聚合起来，而过小又会引入不重要的最频值并人为使得聚类被分裂。 3．4分水岭分割算法 分水岭(也称分水线/水线)分割算法借助地形学概念进行图像分割，近年来得到了广泛使用。该算法的实现可借助一些数学形态学的方法(参见第14章和第15章)。该算法的计算过程是串行的，虽然直接得到的是目标的边界，但在计算过程中利用的是区域一致性。 3．4．1基本原理和步骤 下面先介绍分水岭的概念，再给出对分水岭的计算方法。 1． 分水岭 分水岭是一个地形学的概念。图像也可看成3D地形的表达，即2D的地基(对应图像坐标空间xy)加上第三维的高度(对应图像灰度f)。 参见图3．4．1，设想有一个简单的圆形目标，其灰度从边缘向中间逐步增加，如把图像看成3D地形，则该目标类似一座山峰(圆锥体)，如图3．4．1(a)所示。如果从上向下成像，则得到一幅有圆形目标的图像，如图3．4．1(b)所示。现在考虑有两个圆形的目标，对应两个圆锥体，它们相距很近且有部分重叠，如图3．4．1(c)所示。如果对这两个圆锥体从上向下成像，则得到一幅有两个重叠圆形目标的图像，如图3．4．1(d)所示。假设有水同时从两个山峰上流下来，则水交汇的地方为一条直线，该直线可称为分水岭。在分水岭位置将两个重叠圆形目标分开应可给出一种最优的目标分割结果。 图3．4．1将图像看成地形图用分水岭进行分割 上面借助降水法(水从高下降)引出了分水岭的概念。实际中建立不同目标之间分水岭的过程常借助涨水法(水从低上涨)来讨论。参见图3．4．2(可看作灰度图像的一个剖面)，这里为了简便，仅画出了各个目标的2D剖面。假设有水从各谷底孔涌出并且水位逐渐增高。如果从两个相邻谷底涌出的水的水位高过其间的山峰，这些水就会汇合。如要阻止这些水汇合，就需在该山峰上修筑水坝，且水坝的高度要随水位的上升而增高。这个过程随着全部山峰都被水浸没而结束。在这个过程中修筑的各个水坝将整片土地分割成很多区域，这些水坝就构成了这片土地的分水岭。 图3．4．2涨水法分水岭示意 由上可见，如果能确定出分水岭的位置，则可将图像用一组各自封闭的曲线分割成不同的区域。分水岭图像分割算法就是通过确定分水岭的位置而进行图像分割的。一般考虑到各区域内部像素的灰度比较接近，而相邻区域像素之间的灰度差距比较大，所以可先计算一幅图像的梯度图，再寻找梯度图的分水岭。在梯度图中，小梯度值像素对应区域的内部，大梯度值像素对应区域的边界。分水岭算法要寻找大梯度值像素的位置，也就是寻找分割边界的位置。 例3．4．1用分水岭算法分割接触目标 图3．4．3给出了使用分水岭算法对相接触的圆形目标进行分割的示例。图3．4．3(a)是原始图像，图3．4．3(b)是用分水岭算法分割的结果。需要指出的是，由于目标有可能具有不同的尺寸，所以如果仅使用形态学方法(见第14章)中的膨胀和腐蚀并不能很好地将所有目标分割开，而分水岭算法则对不同尺寸的目标均适用。 图3．4．3对不同尺寸的接触圆目标的分割 □ 2． 分水岭计算步骤 上面的讨论实际上也指出了分水岭计算的思路，即逐渐增加一个灰度阈值，每当它大于一个局部极大值时，就把当时的二值图像(只区分陆地和水域，即大于灰度阈值和小于灰度阈值两部分)与前一个时刻(即灰度阈值取上一个值)的二值图像进行逻辑异或(XOR)操作，从而确定出灰度局部极大值的位置。根据所有灰度局部极大值的位置集合就可以确定分水岭。 根据这个思路，可用不同的方法实现分水岭算法。下面介绍一种借助数学形态学(参见第14章和第15章)中的膨胀运算思想迭代计算分水岭的方法［Gonzalez 2008］。 设给定一幅待分割图像f(x，y)，其梯度图像为g(x，y)，分水岭的计算是在梯度空间进行的。用M1，M2，…，MR表示g(x，y)中各局部极小值的像素位置，C(Mi)为与Mi对应的(分水岭所围绕的)区域中的像素坐标集合。用n表示当前的梯度阈值，T［n］代表记为(u，v)的像素集合，g(u，v)