图像分割的水平集方法

图像分割的水平集方法–理论

在图像分割的传统方法中，基于偏微分方程的方法在众多方法中效果显著。但PDE方程的求解不易，因此有学者针对图像分割的 PDE 开发了水平集方法。本文主要从理论的角度理清水平集方法的来龙去脉，为实际方程中应用水平集方法提供理论支持。

曲线演化问题

对于封闭曲线序列 $C(p,t),t\ge 0$，随时间按如下演化
$$\frac{\partial C(p,t)}{\partial t}=\text{V}=\alpha (p,t)\text{T}+\beta (p,t)\text{N}(1)$$
$\alpha ,\beta$ 分别是切向速率和法向速率。在局部位置，曲线 $C$ 上的坐标可以使用函数表达，如下
$$y=\gamma (x,t)(2)$$
即上式亦可表示曲线 $C$，$\gamma$ 关于 $t$ 的变化也可表示曲线 $C$ 的变化。对于曲线 $C$ 的某一点坐标 $(x,\gamma (x))$,切矢量为 $C_x=(1,{\gamma }_x)$.

单位切矢量为 $\text{T}=\frac{(1,{\gamma }_x)}{\sqrt{1+{\gamma }_x^2}}$

单位法矢量为 $\text{N}=\frac{(-{\gamma }_x,1)}{\sqrt{1+{\gamma }_x^2}}$

则曲线 $C$ 演化时，横纵坐标 $x,y$的变化为
$$\begin{gathered} \frac{dy}{dt}=\alpha {\text{T}}^{(y)}+\beta {\text{N}}^{(y)} \\ =\alpha \frac{{\gamma }_x}{\sqrt{1+{\gamma }_x^2}}+\beta \frac{1}{\sqrt{1+{\gamma }_x^2}} \end{gathered}$$

$$\begin{gathered} \frac{dx}{dt}=\alpha {\text{T}}^{(x)}+\beta {\text{N}}^{(x)} \\ =\alpha \frac{1}{\sqrt{1+{\gamma }_x^2}}-\beta \frac{{\gamma }_x}{\sqrt{1+{\gamma }_x^2}} \end{gathered}$$
然而，另一方面，
$$\frac{dy}{dt}=\frac{dy}{dx}\frac{dx}{dt}+\frac{\partial y}{\partial t}={\gamma }_x\frac{dx}{dt}+{\gamma }_t$$

联立上面三式，有
$${\gamma }_t=\frac{dy}{dt}-{\gamma }_x\frac{dx}{dt}=\beta \sqrt{1+{\gamma }_x^2}$$
$\gamma$ 关于 $t$ 的变化也可表示曲线 $C$ 的变化，从上式可以看出，曲线的变化只与法向分量 $\beta$ 有关。因此研究曲线的演化只需考虑法向。可以简单写为
$$\frac{\partial C}{\partial t}=\beta \text{N}(3)$$
此即为曲线演化方程。但是 $C$ 没有具体表达式， 这个方程不好求解。于是有学者开发了水平集方法。

水平集方法

水平集方法是将空间中的曲线使用二维函数表示，将曲线的变化直接使用函数表示，使得空间的曲线具象化，让我们便于采用数学工具对曲线进行干预和计算。
对于空间中的曲线 $C$，可隐式表达为
$$C=\{(x,y),u(x,y)=c\}$$
曲线 $C$ 称为函数 $u(x,y)$ 的水平集，函数 $u(x,y)$为曲线 $C$ 的嵌入函数。则曲线 $C$ 的某种变化可以归结为函数 $u(x,y)$ 发生的变化。为了表现曲线的变化，加上时间线
$$C=\{(x,y,t),u(x,y,t)=c\}$$
函数 $u(x,y)$ 关于时间的变化为
$$\frac{du}{dt}=\frac{\partial u}{\partial t}+\nabla u\cdot \frac{\partial (x,y)}{\partial t}=0$$
注：趋于稳态时，函数$u$关于时间的变化率为0.
根据上面的曲线演化问题，有$\frac{\partial (x,y)}{\partial t}=\frac{\partial C}{\partial t}=\text{V}$.于是

$\frac{\partial u}{\partial t}=-\nabla u\cdot \text{V}=-|\nabla u|\cdot \frac{\nabla u}{|\nabla u|}\cdot \text{V}=|\nabla u|\cdot \text{N}\cdot \text{V}=\beta |\nabla u|(4)$.
此处 $\text{N}=-\frac{\nabla u}{|\nabla u|}$是按水平集内部取负，外部取正。
$(4)$ 式即为水平集方法的曲线演化方程式。上述的推导过程与$c$值无关，因此我们一般取零水平集，即 $c=0.$

1. 函数 $u(x,y)$如何选取

一般选符号距离函数
$$u(x,y)=\{\begin{array}{r}d[(x,y),C],(x,y)在C外部\\ -d[(x,y),C],(x,y)在C内部\end{array}$$
原因是符号距离函数有如下性质：$|\nabla u|\equiv 1$（这个性质的好处目前未知）
性质证明：$\nabla u=(u_{\xi },u_{\eta })=(0,u_{\eta })$，只有法向变化，因此$|△u|=|△d|=|△\eta |$
$|\nabla u|=|u_{\eta }|=|li{m}_{△\eta \to 0}\frac{△u}{△\eta }|=1$

2. 演化问题

采用水平集演化曲线有一个问题，水平集演化方程与曲线演化方程中 $\beta$ 取值要求不一致。曲线演化方程中的 $\beta$ 的$(x,y)$ 要在曲线上，而水平集掩护方程的 $\beta$ 的 $(x,y)$ 点定义在整个平面 $\Omega$上， 即 ${\beta }_u$是${\beta }_C$的延拓。

直接将${\beta }_u$用作${\beta }_C$也可以，这种方法称为自然延拓，但是自然延拓无法保证演化过程中嵌入函数$u$始终保持为带符号的距离函数，使得演化可能向不期待的方向运行。

当嵌入函数偏离距离函数的性质后，某些局部可能由于$|\nabla u$ 远大于1而出现尖峰或深谷，可能由于$|\nabla u$ 远小于1而出现平坦，这会导致数值计算的迭代过程趋于不稳定。因此若干次迭代后需要重新初始化，把嵌入函数 $u$ 拉回到距离函数性质上来。

这种重新初始化可以通过求解
$$\frac{\partial u}{\partial t}=sgn(u^{(n)})(1-|\nabla u|)$$
完成，此方程达到稳态后，解 $u$ 显然满足距离函数的性质 $|\nabla u|=1$。再将经过此迭代后的函数放入水平集演化中，可继续进行演化。缺点是具体迭代多少次后进行重新初始化，需要通过实验得到。

变分水平集方法

水平集方法先利用变分法最小化曲线 $C$ 的能量泛函，得到关于 $C$ 的运动方程后，再引入嵌入函数，得到关于 $u$ 的PDE.

变分水平集方法是对水平集方法的改进，并不一定是比水平集方法好，适用地方不同。变分水平集方法使用的前提是问题来自于最小化曲线 $C$额能量泛函$E(C)$。不是这样来的就不适用变分水平集方法了。变分水平集方法通过引入嵌入函数 $u$，再利用 Heaviside 函数，将能量泛函 $E(C)$ 改为 $E(u)$ 通过变分法得到关于 $u$的 PDE 从而进行求解。
方法如下：
首先定义特殊函数 Heaviside 函数
$$H(z)=\{\begin{array}{r}1,z\ge 0\\ 0,z<0\end{array}$$
则
$$E(u)={\oint }_{C}g(C)ds={\iint }_{\Omega }g(x,y)|\nabla H(u)|dxdy$$
也就是要求在曲线上进行积分（$\nabla H$只有在曲线上有值，其他地方为0）
$\nabla H(u)=\delta (u)\nabla u,\delta (z)=\frac{dH(z)}{dz}$， $\delta (z)$是Dirac函数。
于是 $E(u)={\iint }_{\Omega }g(x,y)\delta (u)|\nabla u|dxdy$.
利用变分法，
$$\begin{gathered} E(u+ϵv)={\iint }_{\Omega }g\delta (u+ϵv)|\nabla (u+ϵv)|dxdy \\ {E}_{ϵ}^{{}^{\prime }}={\iint }_{\Omega }g\cdot \delta (u+ϵv)\frac{\nabla u+ϵ\nabla v}{|\nabla u+ϵ\nabla v|}\nabla vdxdy \\ =\iint g\cdot \delta (u)\cdot \frac{\nabla u\cdot \nabla v}{|\nabla u|}dxdy \\ =\iint div(g\cdot \delta (u)\cdot \frac{\nabla u}{|\nabla u|}\cdot vdxdy \end{gathered}$$
因此对应的梯度下降流为$\frac{\partial u}{\partial t}=\delta (u)div(g\frac{\nabla u}{|\nabla u|})$.但是在实际计算中，$\delta (u)$ 不好计算，一般用${\delta }_{ϵ}$作为近似。
$${\delta }_{ϵ}(z)=\frac{d}{dz}{H}_{ϵ}(z)$$
这里$H_{ϵ}$是正则化的 Heaviside 函数，满足 $H_{ϵ}(z)\to H(z),whenϵ\to 0$ 即可。
感觉一般还是使用水平集方法多一些。

改进的变分水平集方法

在介绍水平集方法时，由于演化过程中可能会偏离距离函数，因此需要重新初始化。改进的变分水平及方法就是把重新初始化的部分作为整个能量泛函的一部分，在演化过程中同时优化，从而避免了重新初始化。
具体地，定义
$$P(u):={\iint }_{\Omega }\frac{1}{2}(|\nabla u|-1{)}^{2}dxdy$$
利用变分法，
$$\begin{gathered} P(u+ϵv)={\iint }_{\Omega }\frac{1}{2}(|\nabla (u+ϵv)|-1{)}^{2}dxdy \\ {P}_{ϵ=0}^{{}^{\prime }}={\iint }_{\Omega }(|\nabla (u+ϵv)|-1)\frac{\nabla (u+ϵv)}{|\nabla (u+ϵv)|}\nabla vdxdy \\ ={\iint }_{\Omega }div(\nabla u-\frac{\nabla u}{|\nabla u|})vdxdy \end{gathered}$$
因此对应的梯度下降流为$\frac{\partial u}{\partial t}=\Delta u-div(\frac{\nabla u}{|\nabla u|})$.
将此项加入到演化方程中，得到改进的变分水平集方法
$$\frac{\partial u}{\partial t}=\mu [\Delta u-div(\frac{\nabla u}{|\nabla u|})]+{\delta }_{ϵ}(u)div(g\frac{\nabla u}{|\nabla u|})$$

得到水平集演化方程（以上所有方法的统称）后，再利用数值计算的方法（如有限差分法，迎风格式等等）对方程进行离散计算，从而迭代计算得到最终分割结果，即水平集。

参考资料：《图像处理的偏微分方程方法》，王大凯，候榆青等，科学出版社。