
2.2 水平集方法的基本理论
2.2.1 曲线演化理论
曲线演化是利用偏微分进行图像分割的基本理论,它的中心思想是利用几何特征来研究曲线随时间发生的形变情况,曲线的单位法向矢量和曲率是常用的几何参数。
设C(p, t)=(x(p, t), y(p, t))为一条闭合的动态曲线。其中,p为曲线弧长的参数;t为时间参数(t≥0)。曲线C(p, t)随时间演化的过程可用如下的偏微分方程来表示:

式中,T和N分别为单位切矢量和内法矢量;α和β分别为切向速度和法向速度;C0(p)为曲线的初始位置。
给定函数y=γ(x),其含义是以x为自变量的曲线局部y函数,则C(p, t)以x为参数可以表示为

由此可知,曲线C的切矢量为Cx=(1, γx),则曲线C的单位切矢量和单位法矢量可以示为

因此,当曲线C按式(2-12)进行演化时,曲线上任一点的坐标(x, y)按如下方程发生形变:

将y=γ(x)表示为x和t的函数,然后对t求导,可得到

则

上式说明,曲线的形变仅仅受到运动速度法向分量β的影响,据此可进一步去除切向速度α而获得曲线演化的一般方程:

针对曲线演化方程(2-19)的求解一直是学者们研究的难点问题,直到Sethian[139]和Osher等人[140]于20世纪末提出了水平集方法,才使该问题的求解获得了严格的数学理论支撑。
2.2.2 水平集方法
水平集方法[139]是为了解决描述火焰的外形方程而提出的,由于火苗的拓扑结构的不断随意变化,直接用方程是无法描述的,因此提出水平集方法来描述随时间的变化而不断变化的曲线。它的本质是用N+1维的水平集函数来隐含表达N维的空间多维函数(闭合曲线、三维的立体曲面或者多维的超曲面),当嵌入在水平集函数中的闭合曲线(曲面)的拓扑结构发生变化(分裂和合并)时,水平集函数仍能保持其函数的稳定性和有效性。基于水平集的图像分割的本质就是将二维图像空间中的形变曲线隐式地映射到三维的函数曲面ϕ(x, y, t)进行演变,其零水平集就是投影到二维空间的目标轮廓。因此轮廓曲线在演化过程中的重点是在一个偏微分方程的演化作用下,不断地更新水平集函数,计算出演化结束后的水平集函数在该时刻的零水平集,相应得到分割结果,而非时刻追踪标记演化曲线的位置。
在二维平面上,定义ϕ(x, y):Ω→R为水平集函数,零水平集C(s)用来表示活动轮廓线:

在水平集函数中引入一个时间变量t,定义一个随时间变化的水平集函数ϕ(x, y, t),它为闭合曲线C(s, t)在t时刻的隐式表达,即

在水平集函数ϕ(x, y, t)的演化过程中,要保持零水平集的活动轮廓线:

利用微分原理,将水平集函数ϕ对t求导,可以得到

式中,为水平集函数ϕ的梯度。
由水平集函数的定义可得,ϕ在曲线C的切线方向的变化量为0,则

由上式可以看出,ϕ的梯度方向与闭合曲线C的切线方向是垂直的,也就是说,∇ϕ的方向与C的法向量是同方向的。因此,曲线C的单位法向量N可以表示为

可得

式(2-26)被称为曲线演化的水平集方法的偏微分方程式。
运用水平集方法对图像进行分割时,要对式(2-26)进行数值计算。在实际的计算中,首先要对水平集函数进行初始化,一般选取符号距离函数(Function Distance Signed, SDF)作为水平集函数的初始值,主要原因是SDF同时具备简单性和精确性的特点。设SDF为ϕ(x, y, t)=±d(x, y),其中,d(x, y)为图像域中的点(x, y)到给定的闭合曲线C0之间的距离,其符号根据点(x, y)与C0的位置而定,表示如下[141]:

水平集函数初始化采用SDF,其优势在于SDF的基本属性∇ϕ≡1,这表明ϕ(x, y, t)沿着不同方向的变化率在任意点都可保持均匀状态,这种特性使得数值计算保持较高的稳定性。
2.2.3 变分水平集方法
采用曲线演化理论解决图像分割问题,曲线运动方程首先通过最小化关于闭合曲线C的能量泛函,得到C的演化方程,然后引入嵌入函数ϕ,通过求解关于ϕ的偏微分方程,来达到图像分割的目的。
Zhao等人[125]提出了一种新的水平集方法,称之为变分水平集方法[142],是通过引入嵌入函数ϕ和Heaviside函数来改造关于曲线C的能量泛函E(C),得到关于ϕ的能量泛函E(ϕ)。Heaviside函数定义如下:

在具体的计算中,常采用Hε来逼近Heaviside函数;在求解能量泛函过程中,使用Hε的导数δε(z)来进行计算:

变分水平集法和水平集法的主要区别在于:前者通过引入嵌入函数ϕ和Heaviside函数来改造关于曲线C的能量泛函E(C),得到关于ϕ的能量泛函E(ϕ),然后通过变分法得到关于ϕ的偏微分方程;后者直接通过用变分法最小化E(C)来获取关于C的演化方程,接着引入嵌入函数ϕ,得到ϕ偏微分方程。两者比较,前者具有稳定性高、对嵌入函数无须重新初始化等优点,而后者的适用面更广。
基于水平集的方法在演化一段时间后,会受到数字图像离散化的影响和噪声干扰,在一段时间后会发生震荡,使ϕ逐渐偏离SDF,并使轮廓C出现尖角或平坦区域,最终将导致计算结果偏离真实情况和图像分割失败,作为距离函数的水平集函数必须不断地重新初始化,直接代价就是牺牲算法的时间复杂度。为了解决此类问题不再进行初始化,Li等人[98]提出了一种改进的变分水平集方法,在原模型的基础上增加一个惩罚项,可以在水平集演化的过程中完全避免重新初始化。

梯度下降流为

由式(2-31)可以看出,最小化P(ϕ)意味着要使∇ϕ≡1,也就能够使水平集函数在演化过程中尽可能保持距离函数的稳定性,避免再次对其进行初始化,减少了算法的运行时间。在初始化时,可以利用简单的正负常数来代替符号距离函数,即

2.2.4 水平集方法的数值求解
图像分割模型中建立的偏微分方程大多是非线性的,其解析解的求取非常困难,只能借助数值计算方法得到对应偏微分方程的数值逼近或近似解。水平集方法的数值实现的关键就是要把定义于连续空间的偏微分方程进行离散化,使得计算机可执行计算。常见的偏微分方程数值解方法有有限差分法、有限元法、谱法等。在图像处理中,待分割的图像是在二维空间中以离散点的形式表达的,这主要是因为数字图像的空间定义域是二维等间隔离散化的点集,构成了有限差分法所需的等分网格,使用(i, j)作为网格的坐标。因此,有限差分法是图像处理中最为常用的一种数值求解方法。
有限差分法的基本思想是使用离散化的差分近似地代替微分。有限差分法是在离散情况下对于偏微分方程的近似求解,其精度必然低于偏微分方程的解析解。
对于偏导数,当步长Δx足够小时,将其用泰勒公式展开可得到如下式子:

则

式中,O(Δx)为Δx的同价无穷小量。
同理,可得到另外两个展开式:

通过以上分析,可知有限差分有三种格式,即前向差分(forward difference)、后向差分(backward difference)和中心差分(central difference),分别如下:

可以证明,前向差分和后向差分具有一阶精度,中心差分具有二阶精度。对于偏微分方程中的二阶偏导数,可采用中心差分先求两个半点()处的一阶偏导数,得

将上面两个一阶差分作中心差分,得到二阶偏导数的近似表达式为

同理可以得到的表达式,同样使用中心差分计算二阶偏导数为

其中,各“半点”处的函数值近似如下:

结合式(2-43)和式(2-44),可得二阶偏导数的近似表达式为

利用有限差分法的思想近似地表达偏微分方程中的导数,可将偏微分方程转化为计算机可计算的离散化的代数方程。选取一个合适的数值方案,就可以实现方程的演化。根据建立代数方程时选用的方法不同,数值实现方案可分为显式、隐式和半隐式三类。使用一维Burgers方程对几种常用的有限差分法的数值实现方案进行介绍。
1)显式方案
对方程左边进行关于时间变量t的向前差分,对方程右边可以采用关于空间变量x的向前、向后或中心差分。全都使用n时刻的数据进行表达,即

这样就可以求出第n+1时刻的未知数据,这种方案为显式方案。
2)隐式方案
根据上述方程建立差分方程时,方程右侧都使用n+1时刻的数据构造,即

在已知的情况下,得到具有N个未知数的N个方程构成的联立代数方程组,这样的方案被称为隐式方案。
3)半隐式方案
同样在上述方程中,与函数u有关的数据,一部分使用n+1时刻的数据,其余部分使用n时刻的数据,这种方案就是半隐式方案。

半隐式方案得到的代数方程组是线性的,但也能够近似为线性的差分方程,并具有计算简便和高稳定性的优点,因此被广泛应用于偏微分方程的数值计算中。