第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小节)较大的边缘点看作角点。典型的角点检测器仍多基于像素灰度差(梯度)。
利用二阶导数算子检测角点的原理与利用一阶导数算子检测边缘有些类似。在2D图像中，由各个二阶导数组成的对称矩阵(也称为局部结构矩阵)可写成(下标指示偏导数的方向)
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)<G0，其他(3．1．8)
式(3．1．8)是根据核同值区原理获得的，即核同值区的面积越小，边缘的响应就越大。
当图像中没有噪声时，完全不需要几何阈值。但当图像中有噪声时，将阈值G设为3Smax/4可给出最优的噪声消除性能。考虑有一个垂直的阶跃状边缘，则S的值总会在其某一边小于Smax/2。如果边缘是弯曲的，则小于Smax/2的S值会出现在凹的一边。如果边缘不是理想的阶跃状边缘(有一定坡度)，则S的值会更小(R的值会更大)，这样边缘检测不到的可能性更小。
上面介绍的方法一般已可以给出相当好的结果，但一个更稳定的计算C(·;·)的公式是
C(x0,y0; x,y)=exp-f(x0,y0)-f(x,y)T2(3．1．9)
式(3．1．9)对应的曲线如图3．1．4中的曲线b(曲线a对应式(3．1．6))所示。可见，式(3．1．9)给出了式(3．1．6)的一个平滑版本。它允许像素的灰度有一定的变化而不会对C(·;·)造成太大的影响。
图3．1．4不同的C(·;·)函数曲线示例
例3．1．1角点检测实例
图3．1．5给出了两个用最小核同值区算子检测图像中角点得到的结果。
图3．1．5用最小核同值区算子检测到的角点
□
3． 边缘方向的检测
在对边缘的检测中，很多情况下不仅需要考虑边缘的强度，还需要考虑边缘的方向。首先，如果要除去非最大边缘值，则必须要借助边缘的方向。另外，如果需要确定边缘的位置到亚像素精度，常常需要利用边缘方向的信息。最后，在许多应用中常把边缘的位置、强度和方向结合使用。对大多数现有的边缘检测算子，边缘方向是借助边缘强度来确定的。根据核同值区的原理确定边缘方向需根据边缘点的种类采用不同的方法。
根据核同值区的特点，可将边缘点分成两类。参见图3．1．6，其中左图给出图像中一个典型的直线边缘。图中方格里的阴影用来近似地区分具有不同灰度的像素。对3个感兴趣点的局部核同值区分别显示在右边的3个3×3模板中。其中“○”代表核同值区的重心，“+”代表模板的核。区域A和B中都是同一类的边缘点，边缘都通过核同值区的重心，只是模板的核分别落在了边缘的两边。区域C中是另一类边缘点，这里模板的核与核同值区的重心重合。
图3．1．6主要的边缘种类示意
在区域A和区域B中，边缘落在像素之间，从核同值区的重心到模板核的矢量与边缘的局部方向(几乎)垂直。在区域C中，边缘通过像素的中心而不是像素之间的边界，且边缘两边有较高的反差。这种情况对应像素内部边缘，此时，所获得的核同值区是沿边缘方向的细条，如图3．1．6所示。通过寻找最长的对称轴就可以确定边缘方向。具体可通过如下求和公式所得到的结果来估计：
Fx(x0,y0)=∑(x,y)∈N(x0,y0)(x-x0)2C(x0,y0; x,y)(3．1．10)
Fy(x0,y0)=∑(x,y)∈N(x0,y0)(y-y0)2C(x0,y0; x,y)(3．1．11)
Fxy(x0,y0)=∑(x,y)∈N(x0,y0)(x-x0)(y-y0)C(x0,y0; x,y)(3．1．12)
使用Fy(x0，y0)和Fx(x0，y0)的比值就可以确定边缘的朝向，而利用Fxy(x0，y0)的符号可帮助区别对角朝向的边缘具有正的或负的梯度。
剩下的问题就是如何自动地确定哪个图像点属于哪种情况。首先，如果核同值区的面积比模板的直径小，那么应该对应像素内部边缘的情况。如果核同值区的面积比模板的直径大，那么可以确定出核同值区的重心，并可根据像素之间边缘的情况计算边缘的方向。当然，如果重心处在离核不到一个像素的位置，那么这更有可能属于像素内部边缘的情况。当用较大的模板使得处于中间灰度的条带比一个像素还宽时，这种情况就会发生。
4． 最小核同值区算子的特点
与其他边缘和角点检测算子相比，最小核同值区算子有一些独特的地方。
(1) 在用最小核同值区算子对边缘和角点进行检测增强时没有计算微分，这可帮助解释为什么在有噪声时最小核同值区算子的性能会较好。这个特点再加上最小核同值区算子的非线性响应特点都有利于减少噪声。为理解这一点，可考虑一个混有独立分布的高斯噪声的输入信号。只要噪声相对于核同值区面积比较小，就不会影响基于核同值区面积所做的判断。另外，在面积计算中，对各个像素值的求和操作进一步减少了噪声的影响。
(2) 从图3．1．1可以看出，当边缘变得模糊时，在边缘中心的核同值区的面积将减少。所以，最小核同值区算子对边缘的响应将随着边缘的平滑或模糊而增强。这个有趣的现象对一般的边缘检测算子是不常见的，但它在实际中很有用。
(3) 最小核同值区检测算子能提供不依赖于模板尺寸的边缘精度，这与大多数边缘检测算子会随图像或模板尺度的变化而改变所检测出的边缘位置不同。这是因为最小核同值区面积的计算是一个相对的概念，与模板的绝对尺寸无关，所以最小核同值区算子的性能不受模板尺寸的影响。这是一个很有用的期望特性。
(4) 最小核同值区算子选择控制参数很简单，且任意性比较小(只取决于模板的尺寸)，所以比较容易实现自动化选取。
最后，如果对边缘位置的精度要求比用全部像素作为计算单元所能获得的精度还要高，则可采用下面的方法来改进以获得亚像素的精度(可参见4．3节)。对每个边缘点，先确定在该点的边缘方向，然后在与该点边缘垂直的方向上细化边缘。对这样剩下来的边缘点用3个点的二阶曲线来拟合初始的边缘响应，在这条拟合线上的转向点(与细化后边缘点的中心距离应该小于半个像素)可取作边缘的准确位置。
3．1．3哈里斯兴趣点算子
哈里斯兴趣点算子也称哈里斯兴趣点检测器。检测器的表达矩阵可借助图像中局部模板中两个方向的梯度Gx和Gy来定义。一种常用的哈里斯矩阵可写成
H=∑G2x∑GxGy∑GxGy∑G2y(3．1．13)
1． 角点检测
哈里斯角点检测器通过计算像素邻域中灰度值平方差的和来检测角点［Sonka 2008］。在检测角点时，可用下式计算角点强度(注意，行列式det和秩trace都不受坐标轴旋转的影响)：
C=det(H)trace(H)(3．1．14)
理想情况下考虑圆形的局部模板。对模板中只有直线的情况，det(H)=0，所以C=0。如果模板中有一个锐角(两条边之间的夹角小于90°)的角点，如图3．1．7(a)所示，则哈里斯矩阵可写成
H=l2g2sin2θl2g2sinθcosθl2g2sinθcosθl2g2cos2θ+l1g2(3．1．15)
其中，l1和l2分别为角的两条边的长度，g表示边两侧的灰度对比度，在整个模板中为常数。
图3．1．7角点与模板的各种位置关系
根据式(3．1．15)可算得
det(H)=l1l2g4sin2θ(3．1．16)
trace(H)=(l1+l2)g2(3．1．17)
代入式(3．1．14)得到角点强度
C=l1l2l1+l2g2sin2θ(3．1．18)
其中包括3项： 依赖于模板中边长度的强度因子λ=l1l2/(l1+l2)，依赖于灰度差的对比度因子g2，依赖于锐角度数的形状因子sin2θ。
前面已指出，对比度因子是一个常数。形状因子依赖于夹角θ。在θ=π/2时，形状因子取得最大值1； 而在θ=0和θ=π时，形状因子都取得最小值（=0）。由式(3．1．18)可知，对直线，角点强度为零。
强度因子与模板中两条边的长度都有关。如果设l1与l2之和为一个常数L，则强度因子λ=(Ll2-l22)/L，并在 l1=l2=L/2时取得极大值。这表明为获得大的角点强度，需要把角点的两条边对称地放入模板区域，如图3．1．7(b)所示，即角点落在角的中分线(也是圆直径)上。为了获得最大的角点强度，要让角点的两条边在模板区域中都最长，这种情况如图3．1．7(c)所示，即将角点沿中分直径线移动，直到角点落在圆形模板的边界上。
由上面的分析可知，如果角点是直角的角点，则角点强度最大的位置是在圆形模板的边界上，如图3．1．7(d)所示。此时，角点两条边与圆形模板边界的两个交点间的直径是与角平分线垂直的。进一步地，对钝角角点也可进行类似的分析，而结论也是角点两条边和圆形模板边界的两个交点间的直径都是与角平分线垂直的，如图3．1．7(e)所示。
2． 交叉点检测
哈里斯兴趣点算子除了可帮助检测各种角点外，还可帮助检测其他兴趣点，如交叉点和T形交点。这里交叉点可以是两条互相垂直的直线的交点(如图3．1．8(a)所示)，也可以是两条互相不垂直的直线的交点(如图3．1．8(b)所示)。类似地，构成T形交点的两条直线可以互相垂直(如图3．1．8(c)所示)，也可以不互相垂直(如图3．1．8(d)所示)。在图3．1．8中，相同的数字表示所指示的区域具有相同的灰度，不同的数字表示所指示的区域具有不同的灰度。
图3．1．8交叉点和T形交点
在计算交叉点的强度时，仍可以使用式(3．1．18)，只是这里l1和l2的值分别是两个方向直线的总长度(交叉点两边线段之和)。另外需要注意，在交叉点处，沿两个方向的对比度符号都会反转。但这对交叉点强度的计算没有影响，因为在式(3．1．18)中使用了g的平方作为对比度因子。所以，如果使交叉点与模板中心点重合，那么l1和l2的值将分别是角点时的两倍，这也是交叉点强度最大的位置。顺便指出，这个位置是二阶导数的过零点。
T形交点可以看作比角点和交叉点更一般的兴趣点，因为它涉及3个具有不同灰度的区域。为考虑交点处有两种对比度的情况，需要将式(3．1．18)推广为
C=l1l2g21g22l1g21+l2g22sin2θ(3．1．19)
图3．1．9T形交点示例和检测
强度最大的点
T形交点可以有许多不同的构型，这里仅考虑一种一个(斜)弱边缘接触到一个(水平)强边缘但没有穿过该强边缘的情况，如图3．1．9所示。
在图3．1．9中，T形交点是不对称的。由式(3．1．19)可知，其最大值在l1|g1|=l2|g2|时取得。这表明检测强度最大的点处于弱边缘上但不处于强边缘上，如图3．1．9中的圆点所示。换句话说，检测强度最大的点并不处于T形交点的几何位置上，因为受到灰度的影响而产生了偏移。
3．2图 割 方 法
图割方法是一类基于图论的图像分割技术，本质上采用了基于边缘的串行分割思路，所以属于串行边界类。下面先给出用图割方法进行图像分割的主要步骤，再具体对每个步骤进行解释。
用图割方法进行图像分割的主要步骤如下。
(1) 将待分割图像I映射为一个对弧加权的有向图G，G在尺寸上和维数上都与I对应。
(2) 确定目标和背景的种子像素，并针对它们构建两个特殊的图结点，即源结点s和汇结点t； 然后将所有其他像素对应的图结点都与源结点或汇结点分别连接。
(3) 计算弧代价函数，并对图G中的各个弧赋予一定的弧代价。
(4) 使用最大流图优化算法来确定对图G的图割，从而区分对应目标和背景像素的结点。
下面对4个步骤的一些原理和方法分别进行简单介绍。
1． 构建有向图G
一个图结构可表示为G=［N，A］，其中，N是一个有限非空的结点集合，A是一个无序结点对的集合。集合A中的每个结点对(ni，nj)称为一段弧(ni∈N，nj∈N)。如果图中的弧是有向的，即从一个结点指向另一个结点，则称该弧为有向弧，称该图为有向图。对任一段弧(ni，nj)都可定义一个代价(或费用)，记为C(ni，nj)，它可看作对弧的加权。
对给定的待分割图像I，要将其转化表示为一个对弧加权的图G。其中，将图像I中的每个像素看成图G中的一个结点，即结点集合N由所有像素构成； 而将像素之间的邻接关系用图G中的弧来表示，即结点对集合A表示像素之间的(加权)联系。在图G中，需要确定目标种子结点和背景种子结点，它们应分别是最终分割结果中目标集合O和背景集合B的一部分。这个工作目前采用的方法常常借助人机交互来进行。根据所确定的种子，可以对应构建两个特殊的图结点： 源结点s和汇结点t(这里取源结点对应目标种子，汇结点对应背景种子)，它们分别对应源集合S和汇集合T。初始时，源结点和汇结点都与所有对应种子像素和背景像素的结点相连接。
2． 确定种子像素并连接图结点
如上构建的弧加权图可表示成： Gst=［N∪{s，t}，A］，其中结点集N对应图像I中的像素，s和t是两个特殊的终端结点。在Gst中的弧集合A中的元素可以分为两类： 连接一对相邻像素的弧以及将像素和终端结点连接起来的弧。可使用割(cut)将其穿过/跨越的弧切断。在Gst中的一个割集合可将图中的结点分成两组，它的代价是这个割集合所对应的弧集合的代价之和。代价最小的割集合被称为最小st割，它将结点分成两组不重叠的子集S(所有与源连接的结点，s∈S)和T(所有与汇连接的结点，t∈T)，且从s到t没有有向的通路。
图3．2．1是一个解释上面构建有向图和连接图结点的示意图。图3．2．1(a)给出了一幅待分割的图像I，且已确定了目标种子o(属于源集合S)和背景种子b(属于汇集合T)。图3．2．1(b)是所构建的图Gst(其中，用虚线表示连接相邻像素的弧，用不同粗细的实线表示连接像素和终端结点的弧)，其中，oN，bN，o∩b=。在分割开始前，所有结点都同时在源集合和汇集合之中。图3．2．1(c)在图Gst中给出一个st割(将相应的弧割断了)，它将结点分为两组，其中，所有目标结点都仅在源集合中，而所有背景结点都仅在汇集合中。图3．2．1(d)给出了根据这个st割将各个结点映射回图像而得到的分割结果，其中所有的目标像素都标记为白色，所有的背景像素都标记为黑色。
图3．2．1图割分割示意图
3． 计算弧代价
代价最小st割的代价是其所对应的所有弧的代价之和。下面考虑各种弧的代价函数的计算。
给各个像素ip一个二值的标号Lp∈{o，b}，其中，o和b分别代表目标和背景的标号。标号矢量L={L1，L2，…，Lp，…，L|I|}定义所得到的二值分割结果。为计算一段弧的代价，既要考虑该弧两个端结点所对应像素的灰度，也要考虑该弧两个端结点所对应像素之间的灰度差。这里，需要最小化以获得最优标号的代价函数，该函数可以定义为一个用λ加权的区域性质项R(L)和边界性质项F(L)的组合，即
C(L)=λR(L)+F(L)(3．2．1)
其中
R(L)=∑p∈IRp(Lp)(3．2．2)
F(L)=∑(p,q)∈AF(p,q)δ(Lp,Lq)(3．2．3)
δ(Lp,Lq)=1，Lp≠Lq0，其他(3．2．4)
先考虑连接一对相邻像素之间的弧。一方面，给定一个像素，根据其灰度将其标为目标或背景都会有一定的代价。这里，Rp(o)可看作将像素p标为目标的代价，Rp(b)可看作将像素p标为背景的代价。当亮的目标处于暗背景中时，Rp(o)的值在暗像素(低Ip值)处大而在亮像素处小。反之，当暗的目标处在亮背景上时，Rp(b)的值在亮像素(高Ip值)处大而在暗像素处小。另一方面，对两个相邻的像素p和q，根据其灰度将其赋予不同的标号也会有一定的代价。如果它们都属于目标或背景，则弧(p，q)的代价F(p，q)应比较大； 如果它们一个属于目标而另一个属于背景，即跨越目标和背景的边界，则弧(p，q)的代价F(p，q)应比较小。例如，可以取两个相邻像素p和q之间弧(p，q)的代价F(p，q)与它们之间的梯度幅度(边缘两边的绝对灰度差|f(p)-f(q)|)成反比。
再考虑将像素和终端结点连接起来的弧。两个终端结点之间借助通路上对应相邻像素之间的各段弧连接起来的总代价应是这些相邻像素之间弧的代价之和。实际中常常可再加1以使弧没有饱和(B和O分别代表背景像素集合和目标像素集合)：
K=1+maxp∈B∑q: (p,q)∈OF(p,q)(3．2．5)
将上述分析结果结合起来，赋予各种弧的代价函数如表3．2．1所示。
表3．2．1各种弧的代价函数
弧
(p，q)
(s，p)(p，t)
代价C(p，q)(p，q)∈NλRp(b)p∈I，p(O∪B)λRp(o)p∈I，p(O∪B)
Kp∈O0p∈O
0p∈BKp∈B
4． 计算最小st割
计算最小st割的问题可借助它的对偶——计算最大流问题来进行，即可通过搜索从s到t的最大流来解决。式(3．2．5)的弧代价也可解释成从s到p∈O(或从p∈B到t)的最大流容量。在最大流算法中，从s到t的“最大量的水流”通过图中的各段弧构成的通路来输送，而通过各段弧的流量由它的容量或弧代价所决定。从s到t的最大流能使图中的一组弧达到饱和，这些饱和的弧(对应最小割)会将结点分为不重合的两个集合S和T。最大流的值对应最小割的代价。
最小st割问题和它的对偶最大流问题都是传统的组合问题，可用多项式时间算法解决。有许多算法都可用来解决这个组合优化的问题，如增强通路算法、优化最大流的压入和重标记(pushrelabel)算法等［Sonka 2014］。
增强通路算法［Ahuja 1993］考虑推动从s到t的流直至达到最大流量以计算最短的s→t通路。算法开始时将流的状态初始化为0，在将流增强并使通路逐渐饱和的过程中，将流分布的当前状态保存在一个残留图Gr中。当Gr的拓扑与Gst相同时，弧的值在当前流的状态下保留了余下弧的流量。在各个迭代步骤，算法沿着残留图中未饱和的弧所构成的通路来确定最短的s→t通路。流过一条通路的流量借助推动最大可能的流而增加，直到这条通路上的至少一个弧能达到饱和。每个流量增加的步骤都能增加从s到t的总流量。持续增加每条通路的流量直到都不能再增加(不能定义新的s→t通路组成非饱和的弧)，则达到最大流，最优化过程结束。此时的饱和弧集合(对应最小st割)将分别属于S和T的图结点分了开来。在确定最短的s→t通路时可以采用宽度优先搜索的策略。开始先搜索长度最短的s→t通路，一旦所有长度为k的通路都饱和了，则接着搜索长度为k+1的通路，直到穷尽所有长度的通路。可以证明该算法具有收敛性，算法的复杂度是O(|N||A|2)，其中|N|是结点总数而|A|是弧总数。
图3．2．2给出了一个使用增强通路最大流算法来确定最小割的各个步骤的示例。其中，图3．2．2(a)所示为与图3．2．1(a)相对应的原始图像； 图3．2．2(b)所示为根据4连接计算得到的图像灰度的梯度幅度(|f(p)-f(q)|)； 图3．2．2(c)所示为构建图Gst而使用的结点标记，其中源结点记为s，汇结点记为t，其余结点按扫描顺序依次记为a到g。图3．2．2(d)中用粗线表示出饱和图的弧； 图3．2．2(e)加上了分开S和T结点的最小st割； 图3．2．2(f)是最终的图像分割结果(一旦获得了st割，则将结点映射回到图像中即可)。
图3．2．2(g)所示为根据图3．2．2(c)的标记而画出的图Gst，弧的代价标在弧旁； 图3．2．2(h)所示为将图3．2．2(b)的梯度幅度值画在图Gst中的结果； 图3．2．2(i)所示为根据表3．2．1构建的图Gst，其中，λ=0，F(p，q)=max|f(p)-f(q)| - |f(p)-f(q)|(由于有3条弧的流量为0，所以实际上该图中只有一条从s到t的通路)； 图3．2．2(j)所示为当唯一具有非饱和s→t连接的最短通路被确定且饱和后的残留图Gr，此时已没有非饱和的s→t通路了，分割结束。这里，图3．2．2(i)中沿最短通路的流量增加2后从结点a到结点b的弧饱和，则图3．2．2(j)中该通路上各段弧的残留流量都减少2。注意，图3．2．2(j)是与图3．2．2(d)相对应的。
图割方法的一个重要特性是提供了一种借助交互以有效地改进先前获得的分割结果的能力。假设用户确定了初始的种子，且有了代价函数，但图割产生的结果还未达到要求(所得到的割还不能将目标结点和背景结点完全分开)。用户可通过增加目标或背景的种子进行改进。此时不需要从头计算，可以使用先前已经得到的结果来初始化下一个优化过程并继续分割。
图3．2．2利用图割和最大流优化的图像分割过程
图3．2．3最大流计算示意图
例3．2．1增强通路算法示例
使用增强通路算法来计算最大流和最小割的过程可以借助图3．2．3来示例说明。图中给出一个图G={N，A}，即结点集合N和有向弧集合A构成的图。每条连接结点p和q的弧c上的(最大)流量为cpq(非负值)，与弧代价C(p，q)对应。结点集合N有两个特殊的结点，分别是源(起点)s和汇(终点)t。这是一个有向图，在一般情况下允许有些结点之间存在两个方向的弧。最大流计算问题可描述成在图中找到一条从源到汇的通路，这条通路要在各条弧所允许的流量约束下能通过最大的流。换句话说，可以考虑不断增加流量来获得需要的能通过最大流的通路。
这样找到的从源到汇的通路可能不止一条。但在每条这样的通路上，至少有一条弧的流量会达到饱和，否则就可以继续增加流量，因此该通路就不是最大流通路。根据这个分析，最大流计算问题也可从考虑弧的饱和来解决。设“割”是图中将源和汇结点分离开的弧集合，则如果把这个弧集合从图中除去，从源结点到汇结点就没有通路了。每个弧集合对应一个代价，即集合中所有弧的代价(这里对应流量cpq)之和。现在，最大流计算问题就成为找到代价最小的最小割的计算问题了。
图3．2．4给出了一个增强通路算法在计算过程中各个步骤的示例，所用的图G={N，A}与图3．2．3相同，对弧的标记形式为“当前流/饱和流”。它首先从图中任意选择一条从源到汇的通路，并不断增加其中的流量，直到最大流量被通路中允许流量最小的弧所限制(该弧达到饱和)。此时，将这个流量从这条通路中的其他弧的允许流量中减去，饱和弧的允许流量会成为0。重复这个过程，选择第二条从源到汇的通路，增加流量直到有一条弧饱和，更新各条弧的允许流量。这样重复进行，直到再也找不到不含饱和弧的、从源到汇的通路时为止。这些通路的总流量就是从源到汇的最大流量，而这些饱和的弧合起来就构成最小割。如果在选择通路时每次都选择能使剩下的流量达到最大的通路，则该算法可保证收敛。
图3．2．4增强通路算法计算过程
在图3．2．4中，图3．2．4(a)中画出了第一条从源到汇的通路(上部)，包括3条(粗线)弧，即从s到①到④到t。增加流量，直到第三条弧(④到t)饱和(红色)。图3．2．4(b)中画出了第二条从源到汇的通路(中部)，也包括3条(粗线)弧，即从s到②到⑤到t。同样增加流量，这次直到第二条弧(②到⑤)饱和(红色)。在图3．2．4(c)中，将刚才已饱和的弧用另一条弧(②到⑥)替换，得到一条新的通路，也包括3条弧，即从s到②到⑥到t。增加流量，使②到⑥的弧饱和，此时第一条弧(s到②)上的流量是原饱和弧(②到⑤)与现饱和弧(②到⑥)的总和。图3．2．4(d)中画出了第四条从源到汇的通路(下部)，包括3条弧，即从s到③到⑥到t。增加流量，使③到⑥的弧饱和，则第三条弧(⑥到t)上的流量是原经过②到⑥的流量与现经过③到⑥的流量的总和。图3.2.4(e)给出第五条通路，也是唯一的长度为4的通路(前面4条通路的长度均为3)。该通路包括的4条弧为： s到①、①到④、④到⑤、⑤到t，其中第三条弧(④到⑤)先饱和。图3．2．4(f)中共有5条饱和的弧，截断这些弧的割如图3.2.4(f)中的条带所示。该条带将所有从s通往t的通路全部割断，对应最小割。它将所有结点分成了两组： 一组包括源结点s不包括汇结点t，另一组包括汇结点t不包括源结点s。□
图割方法在医学图像分割中得到了广泛应用。例如，［杨2022］采用了一种结合级联DenseUNet和图割的方法分割CT图像中的肝脏肿瘤区域。首先运用级联的DenseUNet获取肿瘤初始分割结果(包含肿瘤组织的感兴趣区域)，然后利用图像像素级和区域级特征，分别构建可有效区分肿瘤与非肿瘤的灰度模型和概率模型，并将其融入图割能量函数，进一步精确分割感兴趣区域中的肿瘤组织。
3．3特色的阈值化和聚类技术
第2章介绍并行区域分割方法时主要讨论了基本的阈值化和聚类方法。本节介绍两种有一定代表性的借助特殊理论和方法确定阈值的技术,以及一种近年常用的确定空间聚类以分割图像的技术。
(1) 借助小波变换的多分辨率特性来帮助进行阈值选取。
(2) 借助过渡区选取阈值进行分割。
(3) 借助均移计算特征空间的聚类中心来进行分割。
3．3．1多分辨率阈值选取
利用图像的直方图帮助选取阈值是常用的方法，其中的关键是确定峰点和谷点。由于场景的复杂性、图像成像过程中各种干扰因素的存在等原因，峰点和谷点的有无检测和位置确定通常比较困难。峰点和谷点的检测与直方图的尺度有密切的联系。一般在较大尺度下常能较可靠地检测到真正的峰点和谷点，但在大尺度下对峰点和谷点的定位不易准确。相反，在较小尺度下对真正峰点和谷点的定位常比较准确，但在小尺度下误检或漏检的比例会增加。
图像在小波变换后可分解为一系列尺度不同的分量(见上册11．4节)。对小波变换后的图像直方图也可进行多分辨率分析(见上册16．3节)，具体就是先在较大尺度下检测出真正的峰点和谷点，再在较小尺度下对这些峰点和谷点进行较精确的定位。这里主要有如下两个步骤。
1． 确定分割区域的类数
首先利用在较大尺度(粗分辨率)下的直方图细节信息确定分割区域的类数。引入相当于低通滤波器的尺度函数(x)，则图像直方图H(x)的低通分量可表示为
S2i［H(x)］=H(x)2i(x)(3．3．1)
设原始图像直方图信号的分辨率为1，最低分辨率尺度为2I，则尺度在21~2I的各阶小波变换为{W2i［H(x)］，1≤i≤I}。可以证明，对信号在尺度为2I时被平滑掉的高频部分可以用尺度在21~2I的小波变换来恢复。这里集合{S2i［H(x)］，W2i［H(x)］，1≤i≤I}就是直方图的多分辨率小波分解表达。
先在分辨率为21时确定初始的分割区域类数，即判断直方图中独立的峰的个数。这里要求独立的峰应满足3个条件：
(1) 具有一定的灰度范围；
(2) 具有一定的峰下面积；
(3) 具有一定的峰谷差。
2． 确定最优阈值
确定类数后，可利用多分辨率的层次结构在直方图的相邻峰之间确定最优阈值。这个过程首先在最低分辨率一层进行，然后逐渐向高层推进，直到最高分辨率层。选择高斯函数作为平滑函数θ(x)，令小波函数为ψ(x)，有
ψ(x)=d2θ(x)dx2(3．3．2)
则ψ(x)对应的二进小波变换为
W2i［H(x)］=(2i)2d2dx2［H(x)θ2i(x)］(3．3．3)
由式（3.3.3）可知，小波变换的零交叉位置对应在分辨率2i时的低通信号H(x)θ2i(x)的剧烈变化点。当尺度2i减小时，信号的局部细节增多； 而当尺度2i增加时，信号中结构较大的轮廓比较明显。
对图像的直方图来说，式(3．3．1)中的S2i［H(x)］是一个最低分辨率下的近似信号，式(3．3．3)中的W2i［H(x)］代表不同分辨率下的细节信号，它们联合构成直方图的多分辨率小波分解表示。给定直方图，可考虑其多分辨率小波分解表示的零交叉点和极值点来确定直方图的峰点和谷点。具体可利用以下4个准则(参见图3．3．1)。
图3．3．1直方图的峰点和谷点的确定
(1) 用从负值变化到正值的零交叉点确定峰的起点(图3．3．1中各个s点)；
(2) 用从正值变化到负值的零交叉点确定峰的终点(图3．3．1中各个e点)；
(3) 用起点和终点间的最大值点确定峰的位置(图3．3．1中各个T点)；
(4) 用前一个峰的终点和后一个峰的起点间的最小值点确定这两个峰之间谷点的位置(图3．3．1中各个B点)。
随着分辨率的逐渐增加，阈值数目也会逐渐增加。这里可用最小距离判据来解决两个相邻尺度之间各阈值并非一一对应的问题。设在两相邻尺度2i+1和2i所对应的阈值分别为Ti+1j和Tik，考察下列条件(Ni是在尺度2i所具有的阈值数目)：
dis(Ti+1j,Tik)=min{dis(Ti+1j,Til)，l=0,1,…,Ni}(3．3．4)
当l=k时取得最小值，这表明在尺度2i+1的阈值Ti+1j对应在尺度2i的阈值Tik。所以，可先对在最低分辨率一层选取的所有阈值逐层跟踪，然后选取在最高分辨率一层的对应阈值作为最优阈值。
3．3．2借助过渡区选择阈值
一般在讨论基于区域和基于边界的算法时认为区域的并集覆盖了整个图像而边界本身是没有宽度的。然而实际数字图像中的边界是有宽度的，它本身也是图像中的一个区域。这是一个特殊的区域，一方面它将不同的区域分隔开来，具有边界的特点； 另一方面它的面积不为零，具有区域的特点。这类特殊区域可称为过渡区。下面介绍一种先计算图像中目标和背景之间的过渡区，再进一步选取阈值进行分割的方法［Zhang 1991c］。
1． 过渡区和有效平均梯度
过渡区可借助对图像有效平均梯度(EAG)的计算和对图像灰度的剪切操作来确定。设以f(x，y)代表2D空间的图像函数，再设g(x，y)代表f(x，y)的梯度图(可通过将梯度算子作用于f(x，y)得到)，则EAG可定义为(Z代表整数集合)
EAG=TGTP(3．3．5)
其中
TG=∑x,y∈Zg(x,y)(3．3．6)
为梯度图的总梯度值，而
TP=∑x,y∈Zp(x,y)(3．3．7)
为非零梯度像素的总数，因为这里p(x，y)定义为
p(x,y)=1，g(x,y)>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)<L(3．3．10)
若对这样剪切后的图像求梯度，则其梯度函数必然与剪切值L有关，由此得到的EAG也变成剪切值L的函数EAG(L)。注意，EAG(L)与剪切的方式也有关，对应低端和高端剪切的EAG(L)分别为EAGlow(L)和EAGhigh(L)。
2． 有效平均梯度的极值点和过渡区边界
图3．3．2对EAGlow(L)曲线是单峰曲线的解释
典型的EAGlow(L)和EAGhigh(L)曲线都是单峰曲线，即它们都各有一个极值，这可以借助对TG和TP的分析得到。这里仅考虑EAGlow(L)，参见图3．3．2［Zhang 1993e］，其中，EAGlow(L)是TGlow(L)与TPlow(L)的比。TGlow(L)和TPlow(L)都随L的增加而减少。TPlow(L)减少是因为大的L会剪切掉更多的像素，而TGlow(L)的减少有两个原因： 一是像素个数的减少，二是剩下的像素之间的对比度也会减少。当L从0开始增加时，TGlow(L)和TPlow(L)曲线都从它们各自的最大值下降。开始时，TGlow(L)曲线下降得比较慢，因为剪切掉的像素都属于背景(梯度较小)； 而TPlow(L)曲线下降得相对较快，因为此时剪切掉的像素个数较多。这两个因素的共同作用会使EAGlow(L)值逐步增加并达到一个极大值。然后，TGlow(L)曲线会比TPlow(L)曲线下降得快，因为更多的具有大梯度的像素会被剪切掉，结果EAGlow(L)值减少，并最后在Lmax处趋向0。
设EAGlow(L)和EAGhigh(L)曲线的极值点分别为Llow和Lhigh，则
Llow=arg{maxL［EAGlow(L)］}(3．3．11)
Lhigh=arg{maxL［EAGhigh(L)］}(3．3．12)
图3．3．3过渡区示例
上面计算出的EAGlow(L)和EAGhigh(L)曲线的极值点对应图像灰度值集合中的两个特殊值，由它们可确定过渡区。事实上，过渡区是一个由两个边界圈定的2D区域，其中像素的灰度值是由两个1D灰度空间的边界灰度值所限定的(见图3．3．3)。这两个边界的灰度值分别是Lhigh和Llow，它们在灰度值上限定了过渡区的灰度范围。
这两个极值点具有3个重要的性质(严格的证明可见［章2001c］或［Zhang 1991c］)。
(1) 对每个过渡区，Llow和Lhigh总是存在并且分别只存在一个。
(2) Llow和Lhigh所对应的灰度值都具有明显的像素特性区别能力。
(3) 对同一个过渡区，Lhigh不会比Llow小，在实际图像中，Lhigh总是大于Llow。
由于过渡区处于目标和背景之间，而目标和背景之间的边界又在过渡区之中，所以可借助过渡区来帮助选取阈值。首先，因为过渡区所包含像素的灰度值一般在目标和背景区域内部像素的灰度值之间，所以可根据这些像素确定一个阈值以进行分割。例如，可取过渡区内像素的平均灰度值或过渡区内像素的直方图的极值。其次，由于Lhigh和Llow限定了边界线灰度值的上下界，阈值还可直接借助它们来计算［Zhang 1993e］。
前面指出的两个极值点的3个重要性质在图像中有不止一个过渡区时也成立。如图3．3．4所示，其中图3．3．4(a)中的剖面有两个阶跃(对应两个过渡区)，反映在梯度曲线上有两个峰。图3．3．4(b)给出了由此得到的有两组极值点的EAGhigh(L)和EAGlow(L)曲线。可以看出，对同一个过渡区，极值点的上面3个重要性质仍成立［章1996c］。所以，上面基于过渡区的方法不仅可用于确定单个阈值对图像进行二值分割，还可以确定多个阈值对图像进行多值分割。另外，由图3．3．4(b)可见，如果将两个过渡区混淆了，则有可能出现Lhigh<Llow的反常情况。根据前面关于极值点性质的讨论，对同一个过渡区Lhigh不会比Llow小，据此可确定对应同一个过渡区的Lhigh和Llow。
彩图
图3．3．4多过渡区时的情况
3．3．3借助均移方法确定聚类
在空间聚类方法中，需要确定聚类的均值或聚类的中心。作为K均值法的一种替代方法，均移方法采用寻找点密度的最频值(峰)的技术，它的一个潜在好处是还可以自动地选择聚类的数目。均移代表偏移的均值矢量，是一种非参数技术，可用于分析复杂的多模特征空间并确定特征聚类。它假设聚类在其中心部分的分布要密，通过迭代计算密度核的均值(对应聚类重心，也是给定窗时的最频值)来达到目的。
下面借助图3．3．5介绍均移方法的原理和步骤，其中各图中的圆点表示2D特征空间(实际可更高维)中的特征点。首先随机选择一个初始的感兴趣区域(初始窗)并确定其重心(如图3．3．5(a)所示)。也可看作以该点为中心画个球(在2D情况下画个圆)。该球或圆的半径应能包含一定数量的数据点，但不能把所有数据点都包进来。接下来，搜索周围点密度更大的感兴趣区域并确定其重心(相当于移动球的中心到一个新的位置，该位置是在这个半径中所有点的平均位置)，然后将窗移动到该重心确定的位置，这里原重心和新重心之间的位移矢量对应的就是均移(如图3．3．5(b)所示)。重复上面的过程，不断地移动均值(结果就是球体会逐步向具有较大密度的区域靠近)，直到收敛(如图3．3．5(c)所示)。这里最后的重心位置确定了局部密度的极大值，即局部概率密度函数的最频值。
彩图
图3．3．5均移方法的原理示意图
在均移方法中，需要确定多变量密度核估计器。这里核函数的作用是使得特征点随着与均值的距离x不同，对均值偏移的贡献也不同。实际中使用的是放射对称的核K(x)，它满足
K(x)=ck(x2)(3．3．13)
其中，c是一个大于0的常数，用来使对K(x)的求和为1； k为剖面核函数。两个典型的核是正态核KN(x)和Epanechnikov核KE(x)。正态核定义为(它常被对称性地截断以获得有限支撑的核)
KN(x)=cexp-12x2(3．3．14)
它的核剖面是
kN(x)=exp-12x2，x≥0(3．3．15)
Epanechnikov(叶帕涅奇尼科夫)核定义为
KE(x)=c(1-x2),x≤10,其他(3．3．16)
它的核剖面在边界处不可微分：
kE(x)=1-x，0≤x≤10，x>1(3．3．17)
给定dD空间的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=1Kx-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． 分水岭
分水岭是一个地形学的概念。图像也可看成3D地形的表达，即2D的地基(对应图像坐标空间xy)加上第三维的高度(对应图像灰度f)。
参见图3．4．1，设想有一个简单的圆形目标，其灰度从边缘向中间逐步增加，如把图像看成3D地形，则该目标类似一座山峰(圆锥体)，如图3．4．1(a)所示。如果从上向下成像，则得到一幅有圆形目标的图像，如图3．4．1(b)所示。现在考虑有两个圆形的目标，对应两个圆锥体，它们相距很近且有部分重叠，如图3．4．1(c)所示。如果对这两个圆锥体从上向下成像，则得到一幅有两个重叠圆形目标的图像，如图3．4．1(d)所示。假设有水同时从两个山峰上流下来，则水交汇的地方为一条直线，该直线可称为分水岭。在分水岭位置将两个重叠圆形目标分开应可给出一种最优的目标分割结果。
图3．4．1将图像看成地形图用分水岭进行分割
上面借助降水法(水从高下降)引出了分水岭的概念。实际中建立不同目标之间分水岭的过程常借助涨水法(水从低上涨)来讨论。参见图3．4．2(可看作灰度图像的一个剖面)，这里为了简便，仅画出了各个目标的2D剖面。假设有水从各谷底孔涌出并且水位逐渐增高。如果从两个相邻谷底涌出的水的水位高过其间的山峰，这些水就会汇合。如要阻止这些水汇合，就需在该山峰上修筑水坝，且水坝的高度要随水位的上升而增高。这个过程随着全部山峰都被水浸没而结束。在这个过程中修筑的各个水坝将整片土地分割成很多区域，这些水坝就构成了这片土地的分水岭。
图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)<n，即
T［n］={(u,v)|g(u,v)<n}(3．4．1)
梯度阈值从图像梯度范围的最低值以整数步长增加。在梯度阈值为n时，算法统计处于平面g(x，y)=n以下的像素集合T［n］。对Mi所在的区域，其中满足条件的坐标集合Cn(Mi)可看作一幅二值图像：
Cn(Mi)=C(Mi)∩T［n］(3．4．2)
即在(x，y)∈C(Mi)且(x，y)∈T［n］的地方，有Cn(Mi)=1，其他地方Cn(Mi)=0。也可这样说，对处于平面g(x，y)=n以下的像素，用“与”操作可将与最低点Mi对应的那些像素提取出来。
图3．4．4用1D的形式(可看作梯度图像的一个剖面)给出了对上面讨论的各个概念的一个直观解释，Cn(Mi)可由C(Mi)和T［n］求交集得到，Cn(Mi+1)可由C(Mi+1)和T［n］求交集得到。
图3．4．4计算Cn(Mi)的示意图
若用C［n］代表在梯度阈值为n时图像中所有满足梯度值小于n的像素集合：
C［n］=∪Ri=1Cn(Mi)(3．4．3)
则C［max+1］将是所有区域的并集(max是图像灰度范围的最高值)：
C［max+1］=∪Ri=1Cmax+1(Mi)(3．4．4)
可以证明，集合Cn(Mi)和T［n］中的所有元素在分割算法进行中，即在将n逐渐增加期间始终保持在原集合中。随着n的增加，这两个集合中的元素个数是单增不减的。所以C［n-1］是C［n］的子集。根据式(3．4．2)和式(3．4．3)，C［n］是T［n］的子集，所以C［n-1］又是T［n］的子集。由此可知，每个C［n-1］中的连通组元都包含在一个T［n］的连通组元中。
分水岭算法先初始化C［min+1］=T［min+1］，然后算法逐步迭代进行。设在步骤n时已建立了C［n-1］，下面考虑从C［n-1］得到C［n］的计算过程。令S代表T［n］中的连通组元集合，对每个连通组元s∈S［n］，有如下3种可能情况(参见图3．4．5)：
(1) s∩C［n-1］是一个空集(见图3．4．5(a))；
(2) s∩C［n-1］中包含C［n-1］中的一个连通组元(见图3．4．5(b))；
(3) s∩C［n-1］中包含C［n-1］中一个以上的连通组元(见图3．4．5(c))。
图3．4．5计算Cn(Mi)的示意图
从C［n-1］得到C［n］的计算过程取决于出现上面3种情况中的哪一种。情况(1)在遇到一个新的极小值时出现，在这种情况下，C［n］可由把连通组元s加到C［n-1］中得到。情况(2)在s属于某些极小值的区域时出现，在这种情况下，C［n］也可由把连通组元s加到C［n-1］中得到。情况(3)在遇到部分或整个能区分两个或以上区域的分界线时出现。如果n继续增加，那么不同的区域将会连通，这时就需要在s中建分水岭。实际中，通过用全是1的3×3结构元素(参见14．3．1小节的讨论)对s∩C［n-1］进行膨胀(限定在s中)就可获得宽度为1的边界。
例3．4．2分水岭算法分割实例
图3．4．6给出了一组对实际沙粒图像用分水岭算法分割的结果。图3．4．6(a)是原始图像，图3．4．6(b)是阈值化后的结果(相邻沙粒混在一起)，图3．4．6(c)是应用分水岭算法分割的结果，图3．4．6(d)是将图3．4．6(c)得到的轮廓线叠加到原始图像上的效果。
图3．4．6分水岭算法分割实例□
3．4．2算法改进和扩展
上面介绍的基本分水岭算法已得到广泛应用，在实用中，还有一些改进和推广。
1． 利用标号控制分割
分水岭分割算法的主要优点是对图像灰度的变化很敏感(逐级检测)，既可以检测出感兴趣区域的轮廓，又可检测出均匀区域中的低对比度变化。在实际中，由于受到图像中噪声和其他局部不规则结构的影响(在梯度图中这些问题更明显)，直接使用前述的算法步骤常常会导致过分割，即将图像分割得过细。直观地说，由于梯度噪声、量化误差及目标内部细密纹理的影响，在平坦区域内部可能会产生许多局部的“谷底”和“山峰”，经分水岭变换后形成小的区域，很容易导致过分割，使希望得到的正确轮廓被大量不相关的轮廓所淹没。
为解决上述过分割的问题，一般可在分水岭分割算法前先加一个预处理步骤(如可借助先验知识)，限定允许的分割区域的数目。
有一种常用的控制过分割的方法是利用标号。一个标号本身是图像中的一个连通组元，可以区分内部标号和外部标号，前者对应目标，后者对应背景。相同内部标号的像素应有相近灰度并组成一个连通的局部极小值区域。运用分水岭算法，将得到的分水岭作为外部标号。这些外部标号将图像分成多个区域，每个区域仅包含一个内部标号和一部分背景。这样分割问题简化成将这些区域一分为二的问题： 一个目标和它的背景。
选择标号有许多方法，简单的仅考虑灰度和连通性，复杂的还可考虑尺寸、形状、纹理、相互位置和距离等。具体可根据应用的先验知识来决定。利用标号可将先验知识加入分割过程，从这个角度说，分水岭算法提供了一个借用先验知识帮助分割的框架。
下面介绍一种利用标号克服过分割的具体方法，称为标号控制的分割［Jhne 2000］。该技术的基本思路是除了对输入图像中的轮廓进行增强以获得分割函数图像(前面求梯度图就是一种特例)外，还用强制最小值技术先对输入图像进行滤波处理。具体就是借助特征检测确定一个标号函数，这个标号函数标示出目标和对应的背景(也可看作标号图像)。将这些标号强制加到分割函数中作为极小值，然后计算分水岭得到分割结果。整个流程可参见图3．4．7，其中细线框代表与外界交互的模块。
图3．4．7标号控制分割方法的流程框图
这里目标的标号可用特征检测的方法从图像中提取出来。对合适特征的选取与先验知识或对图像目标性质的假设有关，常用的包括图像极值、平坦区域等，必要时也可手工确定标号。对每个目标区域都要确定一个标号，因为标号和最终的分割区域是一对一的。标号区域的尺寸可大可小(最小一个像素)，如果处理噪声较大的图像，标号区域也要较大。目标标号确定后，与分割函数图像结合使用。最后通过计算滤波后的分割函数图像的分水岭而得到目标的边界。
前面提到的强制最小值技术是一种基于测绘算子的方法。测绘算子也是一类形态学算子(更多介绍见第14章和第15章)，其主要特点是需使用两幅输入图像。先将一个形态学算子(基本的膨胀和腐蚀算子)作用于第一幅图像，再要求结果保持大于或等于第二幅图像。
在使用强制最小值技术时，认为与图像特征对应的标号已知。对每个像素(x，y)，其标号图像fL(x，y)可定义为
fL(x,y)=0，(x,y)属于标号图像255，其他(3．4．5)
强制最小值的计算分两个步骤进行：
(1) 首先逐点计算输入图像f(x，y)和标号图像fL(x，y)的最小值： fmin(x，y)=min{f(x，y)，fL(x，y)}=f∧fL。这里∧表示取最小值的操作，并在对应标号的地方保留最小值，这样得到的图像将总是小于或等于标号图像。如果需要强制的两个最小值都已经属于f(x，y)在0级时的一个最小值，这时就需要考虑fmin(x，y)+1=min{f(x+1，y+1)，fL(x，y)}+1=(f+1)∧fL而不是f∧fL。
(2) 从标号图像fL(x，y)出发腐蚀fmin(x，y)，直到稳定(类似15．4节中图像聚类快速分割时用到的最终腐蚀)，可以表示为
R*［(f+1)∧fL］(fL)(3．4．6)
图3．4．8给出了对一个1D信号计算强制最小值的示例。图3．4．8(a)中的信号f(x)有7个极小值，而标号信号fm(x)有两个极小值，这两个极小值借助腐蚀加到f(x)上。图3．4．8(b)给出了逐点计算f(x)和fL(x)的最小值得到的结果。图3．4．8(c)给出了从标号图像fL(x)出发腐蚀(f+1)∧fL直到稳定的结果。
图3．4．8强制最小值技术示例
例3．4．3使用标号控制分割
图3．4．9给出了一个使用标号控制分割技术的示例。生物细胞在图像上呈现为块状区域，且常接触或互相重叠。为将这样的区域区分开，可使用标号控制的分割技术。图3．4．9(a)表示部分重叠的两个区域。图3．4．9(b)表示经过距离变换的结果(即1．4节中的距离图)，因为它的局部极小值和需要分开的区域有一对一的关系，所以可用作标号函数(这里取围绕局部极小值的小区域为标号)。从图3．4．9(b)中检测出的分水岭可将两个部分覆盖的区域分割开来，结果见图3．4．9(c)。
图3．4．9分割互相重叠的块状区域□
2． 分水岭算法的扩展
分水岭算法不仅可以在图像域(灰度图和梯度图)中使用，还可以在特征域中使用。下面介绍一种在3D颜色直方图中用分水岭算法找出颜色聚类从而进行彩色图像分割的方法。其主要步骤为［Dai 2003］：
(1) 选择合适的颜色空间(这里选了L*a*b*，主要考虑颜色空间中各颜色之间的相对欧氏距离应与人的视觉特性有尽可能接近正比的关系，参见上册14．2．2小节)，做出待分割图像的3D颜色直方图。实际中为减少计算量，可将颜色空间量化以减少3D颜色直方图中直方条的数量。
(2) 将3D颜色直方图进行反转变换，使其中极小值的位置对应颜色空间中像素个数最多的颜色。对自然图像来说，同一个区域内部的颜色通常比较接近，反映在直方图上会形成山峰形状。当将直方图反转过来后，区域内部的颜色将对应局部极小值。
(3) 在不同的颜色聚类之间建立分水岭，将颜色直方图分割开来，获得颜色聚类。分水岭算法所确定的分水岭位置是直方图中对应图像中不同颜色区域的边界位置，在这些位置能将各个对应不同颜色区域的颜色聚类区分开。
(4) 将聚类结果映射回图像域中，得到图像分割后的各个区域。如果在直方图中找到了颜色聚类的边界，那么在原始图像中也就找到了不同颜色区域的轮廓。
(5) 对上面的结果进行后处理，最终得到分割图像。对颜色丰富的自然图像，常会出现过分割的情况，有些颜色区域非常小，不对应有意义的目标。此时可先确定一个面积阈值将这些过小区域除去。对属于这些小区域的像素可利用诸如最高置信度优先算法等将它们分配到其他区域中。
上述第(3)步中建立分水岭的方法是一种迭代的方法。先将反转直方图中的极小值点依次计算标号，然后取出尚没有标号的最小的直方条，对它根据其邻域(这里需采用3D空间的26邻域)的情况计算标号。如果所考虑的最小直方条的26邻域中出现了多于一种的标号，则可判定它对应颜色聚类之间的分水岭； 否则将其邻域中已标记直方条的标号值赋给它。这个过程反复进行，最后就可找出所有分水岭的位置。
例3．4．4分水岭算法分割彩色图像
图3．4．10给出了几组用分水岭算法分割彩色图像得到的结果［Dai 2006］。第1列是一组原始图像，第2列是后处理前的结果，第3列是后处理后的结果。□
图3．4．10分水岭算法分割彩色图像的结果
3． 分水岭方法与深度学习技术
深度学习技术与分水岭方法也有许多交集。例如，［张2022a］在根据高分遥感影像对地震灾区倒损建筑进行的研究中，比较了基于分水岭分割和基于UNet分割的效果。基于分水岭分割中，先去除阴影,再采取分水岭算法对建筑物分割分类，对建筑物倒损程度进行评估。基于UNet分割中，先对数据集进行训练，得到建筑物分割结果,再对地震前后分割结果做变化检测，进行建筑物倒损程度评估。实验结果： 对建筑物提取精度的评价如表3．4．1所示； 对建筑物倒损识别精度的评价如表3．4．2所示。该对比表明UNet分割的效果更优(Kappa系数是一种衡量分类精度的指标，其值越大说明精度越高，可用于一致性检验)。
表3．4．1建筑物提取精度的评价表
精度评价方法
分水岭分割
UNet分割
总体精度0．6900．803
用户精度0．7320．782
Kappa系数0．7150．837
表3．4．2建筑物倒损识别精度的评价表
精度评价方法
分水岭分割
UNet分割
总体精度0．7190．823
用户精度0．6940．841
Kappa系数0．7340．845
彩图
深度学习技术与分水岭方法也可以结合。例如，［朱2022］提出了一种基于深度卷积网络和分水岭分割的方法用于分割耕地的地块。实际中，基于深度卷积模型直接识别耕地区域时会丢失内部边界，而基于边缘检测模型识别耕地边界时则会同时得到大量无关边界； 此外，基于阈值提取地块的策略所提取的地块不够规整，存在内陷的问题。结合深度网络和分水岭进行的改进如下。
(1) 将耕地边界视作一种地物类别，在深度卷积网络中进行类别概率检测，帮助实现对耕地边界的语义识别。
(2) 基于改进后的DLinknetXt网络进行检测，其网络架构适合对耕地边界这类线性目标的提取，同时更换原始DLinknet网络的残差单元，有助于网络的特征提取能力。
(3) 基于分水岭分割对耕地地块进行提取，利用了区域分割方法获取封闭的边界。这种以区域为单元进行分割及合并的方式，解决了原有方法在像元尺度上基于阈值提取所遇到的提取地块存在内陷的问题，使地块更加规整准确。
总结和复习
随堂测试