水平集算法笔记

水平集算法报告

水平集方法

 水平集方法的主要思想是将不断变形的维曲线（面）作为零水平集嵌入到更高一维的维函数中，由封闭超曲面的演化方程可以得到函数的演化方程，而嵌入的封闭曲线（面）则总是保持为函数的零水平截面上的点集，因此，只要获取演化函数在零水平截面上点集的位置，也就得到了变形曲线（面）的演化结果。

  如图1，一条曲线在运动过程中分裂为两条曲线，这样的拓扑结构变化不可能用一条连续的参数化曲线的运动来表示。然而，上述曲线的拓扑结构变化却可以表示成一个连续变化的曲面与一个固定的平面的交线的变化，曲面本身不发生拓扑结构的变化，从而使复杂的曲线运动过程变为更高一维的函数的演化过程。

图1

  该图同时也模拟了水平集演化的过程。给定一个初始的封闭曲线，该曲线可以沿着其法线方向朝外或朝内以一定的速度演化，引入时间变量，则可以形成随时间变化的曲线族。可以把该曲线族看作是更高维空间曲面函数的零水平集，此时（水平集函数）在切面上的虚线即为隐含的曲线，其对应着二维平面上的封闭曲线。由于是以一种隐含方式来表达闭合曲线的，因此在处理该曲线的演化问题时，水平集方法并不是去追踪演化后曲线的位置，而是遵循演化规律即水平集演化方程，来不断的更新水平集函数，从而使零水平集发生变化，从而达到演化隐含在水平集函数中闭合曲线的目的。这种方式的最大优势是：即使隐含在水平集函数中的闭合曲线发生了拓扑结构的变化，水平集函数仍然保持为一个有效的函数，并存在稳定的解。同样，这一策略对于更高维度的空间同样适用。

推导过程

设是单位法向矢量，则曲线演化方程可以写为：

                                 （1）

     给定水平集函数，则时刻演化的曲线为其零水平集，由复合求导法则：

                            （2）

其中为的梯度。

     设是的弧长参数，根据曲线只沿其法线方向演化，沿着切线方向的变化量，即：

                            （3）

由（2）式可知，垂直于切线，即和法线在同一方向上。如果规定在零水平集内为负值，在外为正值，则水平集曲线的内单位法向量为：

                                  （4）

将（4）式与（1）式带入（2）式得如下的水平集演化方程：

                                  （5）

由（5）式可知，水平集方法的实质就是求解一个随时间变化的微分方程。为了保证计算的精度，水平集函数需要具备以下两个性质：

   （1）水平集函数的光滑性。由于演化时需要计算水平集函数的梯度，而梯度对数值变化是非常敏感的。并且大部分方法中的速度函数都会包含有常数项和曲率项，曲率项在数值求解时需要计算水平集函数的二阶导数，常数项也要计算一阶导数。因此水平集函数的光滑性对计算精度有重大影响。

   （2）水平集函数应保持为符号距离函数。一方面，曲线的隐式表达式与它的符号距离函数是等价的。另一方面，从数值计算的角度，由于符号距离函数的梯度，由此可以保证离散网络的大小为，使得数值计算获得较高的精度。

综合以上两点，在水平集开始演化之前，需要把水平集函数初始化为符号距离函数。以二维空间为例，设为初始曲线，则零时刻的水平集的函数为：

               (6)

其中是一个取值为或的符号函数，如果点在的外（内）部，则，反之。表示点到初始曲线的最短距离。

对于规则的初始曲线如圆或矩形来说，距离函数的求解相对比较简单。如果对于任意形状的初始曲线来说，直接计算距离函数的复杂度还是相当大的，并且在演化一段时间过后，水平集函数会发生震荡从而失去光滑性和距离函数的特性。为了解决这个问题，我们每隔一段时间便需要重构水平集函数。

数值计算

由于水平集方法在演化过程中始终保持为一个有效函数，因此可以使用离散网格的形式来表示水平集函数。以二维函数为例，设离散网格的间隔为，时间步长为，则在时刻网格点处的水平集函数为（缩写为），则演化方程（5）可以离散为：

                                        （7）

其中表示时刻速度函数在网格点处的值。对于（7）式可以采用迎风有限差分法进行求解。

 

首先定义一阶中心差分、一阶前向差分和一阶后向差分六个算子如下：

                    （8）

则（7）式可以改写为：

                  （9）

其中和可分别表示如下：

  （10）

 将（9）式应用于图像分割时，速度函数可以写成如下形式：

                              （11）

其中代表常量演化速度，是曲率演化速度，是水平对流速度，。因此（9）式可进一步写为：

       （12）

通过（12）式可以不断更新水平集函数。因为使用有限步长法，所以需要时间步长的选择，即在给定网络间隔的情况下，要满足如下CFL（Courant-Friendrichs-Levy）条件：

                            （13）

为了保证演化的稳定性和收敛性，CFL条件对于时间步长和速度给出了一个上限。    

值得注意的是式（12）是在图像全局域上进行计算的，因此对于尺度较大的图像来说其计算量是非常大的。为了解决这个问题，我们利用窄带方法实现水平集演化的快速数值计算。

     

            （a）                             （b）

图2，窄带示意图

图2给出了窄带的示意图，图中的黑线即为零水平集，图2（a）是三维曲面上的窄带示意图，图2（b）对应二维平面。在数值演化过程中只更新零水平集周围即图中阴影区域内水平集函数的值，这可以使演化过程的计算量大大减小。

窄带法存在一个问题就是：在若干次迭代之后，零水平集可能会落在窄带范围之外，此时就需要更新窄带水平集函数，这包括两个步骤：

（1）更新窄带内外边界，使得零水平集重新落在窄带区域内，同时保持窄带宽度为2；

（2）重新初始化水平集函数的符号距离函数。

窄带更新的次数与窄带宽度有关，如果宽度过小，则水平集可能会经过很少的迭代次数之后就落在窄带之外，导致窄带需要重新更新。如果窄带过宽，则需要计算的数值个数会大量增加。因此，在演化之前需要选择合适的窄带宽度，比较合理的取值范围为。