2.4 Allen-Cahn去噪模型的水平集求解

水平集方法的主要思想是把二维平面的闭合曲线视为三维空间连续函数曲面φ（水平集函数）的一个水平层，通常是{φ=0}的一个水平层。这样，二维曲线的演化就转换为三维水平集函数曲面的演化。要使得水平集函数曲面的演化结果和闭合曲线的演化相关，水平集函数曲面的演化要遵循如下Hamilton-Jacobi类型的偏微分方程：

式中，F是和曲线在其法线方向的速度相关的某种数量，通常情况下F只在曲线的位置上有意义，也就是定义在零水平层上，因此需要将零水平层的F扩展到整个函数曲面定义域，是扩展速度场Fext。一般来说，有两种速度扩展的方法。一种是全局扩展算法，另一种是窄带算法，速度比全局算法快。这两种算法的过程基本上是一样的，只是要处理的数据量不一样。窄带算法需要增加与限制数据量的相关操作。窄带法最初由Chop提出，Adalsteinsson给出了详细的实现方法。基本思想是只演化位于零水平集周围很窄的一个带状区域水平集函数的值。基本的过程描述如下：

①初始化平面上的所有点，包括计算初始距离和执行速度扩展；

②用速度函数作迭代计算；

③跟踪曲线经过一步迭代后的新位置；

④检查终止条件是否满足，如果不满足则进行新一次的迭代计算。

水平集方法的计算速度是决定它使用范围的一个重要因素，因为它将影响到整个图像处理过程的进度，而初始化步骤是水平集方法中最为耗时的一个过程。初始化过程包括对图像中所有点的初始距离的计算和每个点在初始曲线上的最相近点的确定。即使在使用窄带方案时，也要经过几步迭代后才作重新初始化计算。如果图像数据量很大，整个计算时间会相当可观，这样水平集方法就不能用于实时计算了。

2.4.1 快速行进算法

在窄带算法的计算过程中，第一步就是初始化步骤。一种最简单的初始化方法就是通过计算当前点与曲线上所有点的距离，寻找距离的最小值，以得到与当前点距离最近的曲线上的点，同时记录下最近点的坐标。这种方法称为直接法。直接法的计算速度为O（N×n），这里N是平面上点的数量，n是曲线上点的数量。Adalsteinsson提出了用快速行进法为水平集作初始化计算，这是一种用来求解Eikonal方程的数值算法。算法中结合了Hamilton-Jacobi方程的黏性解法中的上推算法、水平集方法中的窄带算法和一个最小数据堆的数据结构。

设C是一个初始曲面，F是沿曲面法线方向的速度。考虑当曲线经过点（x，y，z）的特殊情况，这时函数T（x，y，z）满足方程|∇F|T=1。根据Osher-Sethian求解水平集的“熵守桓”差分方法，可得方程|∇F|T=1的数值解表达式如下：

式中，D-和D+分别是前向差分算子和后向差分算子。由于上式从本质上说是每个点的二次方程，因此可以用迭代的方法求解方程，然后选择最大可能的解作为方程最后的解，这和水平集方程的黏性解法是一致的。

由式（2.12）可知，边界传播方向是从T比较小的点流向T比较大的点，根据这样的特点，Sethian提出了快速行进算法来迅速传播边界值T。基本思想是在传播边界外围构造一个激活（Active）窄带，窄带内的点的到达时间未定，利用迎风机制（Upwind Scheme）将当前边界向外传播，就像水波扩散一样，凡是扩散到的点，就冻结其波前到达时间，然后再根据当前的波前构造新的激活带，如此循环，就可以得到整个平面上每点的到达时间。因此，建立快速行进算法的关键是：用上推的方式扫描当前曲线周围一个窄带范围的点，得到曲线演化最先到达的点，从而把曲线不断地向前推进。根据计算网格每点距离闭合曲线的远近，将网格点分为三类：

①激活点（Alive），已经被波前传播过的点，到达时间已知；

②活跃点（Active），当前波前即将经过的点，到达时间未知；

③远点（Faraway），远离波前的点，到达时间未知。

快速算法的主要技术是一个最小堆的数据结构。这个结构实际上是一个完全二叉树，二叉树中每个父节点的值要不大于其子节点的值。在实际应用中，一般会把二叉树存储在一个数组中，而子节点的位置正好是父节点的两倍。因此，存储最小值的根节点是数组的第一个元素。寻找父节点或子节点的位置只要花费O（1）的时间。在存储时间T的同时，需要存储点在网格中的位置坐标。快速行进算法通过不断寻找窄带中的最小值（FindSmallest）而达到求解的目的。FindSmallest操作包括去掉根结点和一次对二叉树的下推（DownHeap）扫描，DownHeap操作是为了保证剩余的点满足二叉树的性质。算法步骤为：标记当前点的邻点中当前状态是FarAway的状态为Trial，这是通过插入（Insert）操作完成的；同时还要重新计算当前点的邻点的到达时间。Insert增加节点的数量，并且用上推（UpHeap）操作把插入点放在合适的位置以保证二叉树的性质。快速行进算法中的很重要一点在于需要在计算过程中始终保持二叉树组的性质。这可以使我们在查找窄带中最小值点的过程中值只需要O（1）的时间就可以完成。这样，遍历整个图像点的过程需要花费O（Nlog（N））的时间。

2.4.2 Hermes算法

FMM算法是一种有效降低水平集计算开销的方法，它通过计算Eikonal方程得到网格点的行进速度，通过这种方式可以使计算复杂度从O（N2）降至O（Nk），其中N为水平（或垂直）方向的网格点数目，k为窄带的宽度，但该算法要求演化种子点的行进方向保持不变。Paragios提出的Hermes算法对经典FMM算法进行改进，在保留了FMM算法快速优点的同时，解除了行进速度方向保持不变的限制。

本章将速度函数应用到Hermes算法中，设计出式（2.10）的水平集方程数值算法。

式（2.10）在Hermes算法中的离散差分为

式中，h和τ分别为空间步长和时间步长，V为演化曲线的速度函数。对流项∇g·∇u的迎风型有限差分逼近和Δu项的拉普拉斯差分逼近为

式中，，D+和D-分别表示前向差分和后向差分算子。

行进速度函数V的有限差分表示为

2.4.3 快速水平集演化算法

根据Hermes算法和演化曲线速度函数V的描述，设计APNACM算法如下。

STEP1.输入含噪图像，设定移动界面宽度参数ξ和积分阈值参数γ。

STEP2.含噪输入图像中的所有网格点标记为远离点，并将速度值v置为∞。

STEP3.计算初始图像网格点（ih，jh）的u值，确定初始活动点集。如果网格点（ih，jh）的u≤0.5，且与其相邻的4个网格点中的任何一个u＞0.5，则认为该点是初始图像曲线所在网格点，并将该点标记为活动点。

STEP4.提取STEP3中生成的一个活动点v，根据式（2.16）依次对v相邻的4个网格点wi（i=1，2，3，4）计算wi的速度值vi：如果wi是窄带点，则更新窄带点wi的速度值；如果wi不是窄带点，则标记该点为窄带点，将其插入窄带点集。

STEP5.如果STEP3中生成的所有活动点都处理完毕，则执行STEP6；否则返回STEP3。

STEP6.从窄带点集中取出速度最快的窄带点w作为演化种子点，并标记w为活动点。

STEP7.如果水平集曲线积分值小于阈值参数γ，则转至STEP9，否则依次对w点相邻的4个网格点wi（i=1，2，3，4）进行标记和行进速度的更新：如果①wi是活动点，则返回STEP7；②wi是窄带点，则根据式（2.16）计算wi的新速度值vi，同时更新窄带点集中wi的相应速度值vi；③wi不是窄带点，则根据式（2.16）计算wi的新速度值vi，同时标记wi为窄带点，将其插入窄带点集。

STEP8.如果窄带点集不空，则返回STEP6，继续下一个窄带点演化；否则执行STEP9。

STEP9.输出去噪图像，结束。