关于水平集函数的重新初始化过程

转载自：水平集

 关于水平1.对于曲线的水平集演化的理解1.1 基于水平集方法的曲线演化原理^[1][2]

曲线演化理论（Curve Evolution）是指仅利用曲线的单位法向量和曲率等几何测度来研究曲线随时间形变的理论。其基本公式如下：

                      (1)

其中，C(s,t) = (x(s,t),y(s,t))，为动态曲线，s表示曲线的参变量，t表示时间变量；T 表示单位切向量，N表示向内单位法向量，α，β分别表示切向量和法向量的形变速度。如果曲线演化过程中指关心曲线的形状和几何特征，这(1)可以表示为：

                          (2)

若β= ak，其中k为曲率，则曲线演化过程中就会变的趋于光滑。

定义与C相关的水平集函数：

           (3)

则C的水平集表示为：

                            (4)

则对t求导，利用全微分有：

                           (5)

将(2)式代入，又则有：

                             (6)

推导到这可以得到，与曲线演化相对应的水平集演化函数。其中，V = β=a k(φ)。k 是水平集函数的泛函。

1.2  水平集方法与曲线演化的关系

  事实上，利用水平集方法进行曲线演化其关键在于公式(3)及(5)。公式(5)建立了水平集函数演化的条件，即曲线演化。图1.(a)显示了曲线的水平集演化过程，表示了曲线演化和水平集函数演化的关系。水平集函数一般初始化为与C对应的距离函数如公式(3)所示。

图 1 曲线的水平集演化过程

2. 重新初始化水平集函数的原因及其实现 2.1    需要进行重新初始化的原因

1.2所述过程表明演化过程的关键在于：

条件1：φⁱ =0能够精确地表示Cⁱ

然而，容易忽略的事实是，水平集函数经过演化方程(6) 以后就不再是距离函数（但是其零水平集当期曲线，只不过整体不是距离函数）。由于1）曲线演化过程中β=a k(φ)，是关于φ的函数（曲率k与φ之间的关系，以及φ不为距离函数对k的影响还要进一步看文献）。因此不正确的φ，就会导致不正确的曲线演化，进行导致不正确的水平集演化，形成恶性循环；2）直接影响公式(6)的演化，事实上会造成V在不同位置的扭曲（因为乘上了|∆φ|，而|∆φ|不等于1，造成在不同位置上V受到|∆φ|的影响）。

为了克服这一影响，需要在水平集函数演化结束后对其进行重新初始化，将其变化为距离函数，如图1.(b)所示。φ成为距离函数满足的条件是：（宏观理解，距离函数很容易导出条件2，但是条件2而是不是充分条件，不是很清楚，这一部分数学东西就深了）

条件2：

2.2    重新初始化水平集函数的实现

为了利用水平集方法进行精确的曲线演化，重新初始化水平集函数必须同时满足条件1和2，文献[3]提出如下重新初始化方法：

                         (7)

其中，是指经过公式(6)演化后的结果，S定义为：

                            (8)

方程(8)，保证了条件1的实现，即重新初始化后的水平集函数与演化后的水平集函数的零水平集保持一致，也就符号函数起的作用了。而保证了条件2的实现。

文献[3]中阐述了其数值实现，具体：

                       (9)

                     (10)

     (11)

                  (12)

其中， 可以解释为什么公式(11)中要减1，对照公式(12)和(7)可以得到结论。而对于(11)中取a⁺, b^-以及max的应用，事实上有克服干扰的作用：在j固定的情况下，φ是一个关于i的双曲线函数，如图2所示，取最大值方法可以有效的去除一些扰动（自己理解，没找到相关资料）。在i固定的情况下有类似的结果。

 

(a)                               (b)

图2 max的作用体现

2.3    相应 matlab 函数的解释

见creaseg 软件中的creaseg_chanvese.m 文件中的sussman函数。解释如下：

D = sussman(D,dt)

输入：D 为要进行重新初始化的水平集函数；

dt 为时间步长；时间步长的确定是根据库朗条件设定的，在cv模型的文档将具体涉及。

 输出：D 为处理后的水平集函数；

 

整体思路与公式(9)-(10)是完全一致的，可以很容易看懂。下面介绍下里面用到的函数：

求φ的偏导数用到的平移函数，例如shiftL 其作用如图3所示，其他平移函数都有类似作用可以得到a b c d

图 3 b的获取

  a_p 和 a_n 就是分别存储a 中大于零和小于零的值，不满足各自条件的设为0。就可以得到公式中的a⁺ a^-。

转载请标注

[1]     朱国普. 基于活动轮廓模型的图像分割. 哈尔滨工业大学博士论文,导师:曾庆双.2007

[2]     S.Osher, J.A.Sethian. Fronts Propagating withCurvature-Dependent Speed:Algorithms Based on Hamilton-Jacobi Formulations.Journal of Computational Physics. 1988,7(3):12~49.

[3]     M. Sussman, P. Smereka, S. Osher . a level ste approachfor computing solutions to incompressible two-phase flow . Journal ofComputational Physics, 1994, 114(1): 146-159.