水平集(level set)算法原理介绍

本篇文章，解释的是水平集算法最基础的原理。

• 1

水平集方法的解释 
有一个表面S，它与一个平面P相交，得到一个曲线C，这个C就是我们通过水平集得到的轮廓。 
在图像分割中，表面S是随着由图像派生得到的势（force）来更新。 
本文的思路是： 
1提出问题 
2提出解决方法 
3方法的局限性

跟踪界面 
首先，我们来想象水从一个小山的山顶往下流的画面。我们的目标是，在水往下流的时候，跟踪水前（water front，每一个点的水下一秒流动的方向叫作水前，我们跟踪的是所有点的水前。并且我们的初始条件是，所有初始位置点构成的曲线）。

如下图所示，这是一个V型小山，黑色部分是山顶，越往白色的地方，海拔越低。 
 
红色箭头是水前。

问题来了，如何在给定时间T时刻，得到所有的水前？

显式轮廓方法 
一个方法是，跟踪一系列点，沿着它们的法线方向（normal direction）来演化水前，并且猜测水前的停止位置。

如下图所示，在t=0时刻点集的水前位置，当t=1的时候，点集运动到了虚线位置。 

显式演化一个点集的水前，听起来是个好主意，但是它有很大的缺陷，使得我们不得不放弃这种可以直观理解的方法。

那究竟是什么缺陷呢？ 
如下图所示，这个倒V型的曲线，在它的顶端是一个角（corner）。离角最近的那两个点有着各自的水前，它们方向不一致。随着水前的演化，该角最后的状态将无法得知。红色是0时刻的点状态，粉色是1时刻的点状态。 
 
另外，如果水前是向外扩展的，那么初始的点集可能不够用来得到演化后的水前。需要使用插值法或者在演化出现冲突的时候删掉一些点。并且，点之间的距离要保证足够小，才能得到平滑的水前。这种得到演化水前的机制在应用时是很麻烦的。

在应用上，拓扑变化需要关注的问题更多，比如在分裂和合并的情况下。在我们的小山流水的例子中，左图的中间的红色矩形方框是初始点集，圈住的是山顶最高位置，即将向上下的白色山谷流下去。 
在右图中，我们看到，初始点集演化到一定阶段，就会分裂，并且分裂后会与本来在山谷里的曲线合并~问题很明白了，这个方法要解决的挑战是，如何确定在分裂时该插入的点，以及在分裂时该删掉的点呢？ 
  
所以，追踪点集的水前，尽管是一个很直观的办法，但却有很难解决的缺陷。

隐式轮廓方法 
这个方法的本质是演化一个表面S，而不是一个水前曲线C，且我们用隐式轮廓法得到的水前定义为这个表面S在高度h=0的所有点构成的集合（没有删除和插值这么麻烦的事情了，开心~）。于是根据定义，水前曲线就是零水平集φ=0。如下图所示，是从一个演化的表面得到轮廓。 
 
当表面φ（x,y,t）在演化时，它将呈现杯状，过后杯口或将变窄，或将消失。下图中的零水平集显示了轮廓曲线合并和分裂的情况，z=0表示的是山谷。 
  
（a）t=50，合并开始 （b）t=52，合并结束 
  
（c）t=90，分裂开始 （d）t=120，分裂结束

对于拓扑变化，我们不用再过多关注其它问题了。 
那么问题在于，我们如何确定函数φ呢？

水平集方程 
首先，我们来看看这个方法背后的数学理论~ 
点x(x,y)属于一个随着时间演化的曲线，x(t)为它在t时刻的位置。在任意时刻t，对于每一个点x(t)都是表面φ在高度为0的曲线上的点，即： 
φ(x(t),t) = 0

问题还是，这个函数φ到底是什么呢？ 
只要它给得出我们需要的零水平集，那么它可以是任何定义。。。（这就糟了，怎么造呢）

举个栗子，上面有一幅图，显示的初始轮廓是一个矩形。该表面的高度等于（x,y）到轮廓上点的最近距离，于是有φ(x,y,t=0)=±d，+d表示在轮廓外，-d表示在轮廓内。只要φ的零水平集可以匹配初始轮廓，那它可以被定义成任何函数。 
在时刻t=0处，给定初始函数φ，根据运动方程∂φ/∂t我们可以得到任意t时刻的φ。对于此，利用链式法则，有： 
 
把∂φ/∂x记作▽φ，速度x_t 的方向由表面的法向量给定，因此x_t =F(x(t))n，其中。上面的运动方程可重写为： 
 
最后一个方程定义了φ的运动。给定t=0时的φ以及它随时间演化的运动方程，我们可以通过演化初始函数得知任何t时刻的。这样我们就回答了最初的问题，即我们知道了φ是什么。 
φ有个有趣的特征，就是我们可以用下式得到表面的曲率： 
 
我们可以用这个来控制水前曲线的平滑度。

应用 
在计算机的世界，图像是具有的像素，而函数需要被离散化。也就是说，在 像素(i,j)处的值由估计。这个梯度由有限差分法（finite difference scheme）来估计,比如下： 
 
其中，对于给定点x，表示左有限差分，表示右有限差分。如下图所示，梯度计算的方式因水前方向的不同而不同，有点差分法也会考虑到水前方向。 
 
于是，之前的运动方程可化作： 
 
从此，我们可以用下式来更新表面： 
 
当我们计算曲率的时候，只依赖于表面φ。于是可以使用中心差分法： 
 
曲率的数值计算可利用下式： 

为了使水前曲线平滑，高的曲率应该被惩罚。为了使零水平集扩展，我们不通过减小φ，而是利用高曲率使得水前运动反向。我们可以利用曲率来更新： 

结果 
给定任意的初始函数φ，比如一个初始轮廓的距离转换，的数值计算方法，我们来展示一些轮廓演化的例子。 
第一个小栗子，一滴在扩展的水，前面有一个小障碍物，如下图所示： 
 
水前运动被障碍物停止下来，并被扰动。 
下面有个更复杂的例子。 
 
初始轮廓仍旧是个圆~~随着表面的演化，轮廓开始分裂，从第三图可看到已经完全分裂。

下面这个例子，利用真实的图片。不使用连续的正势和负势，有无曲率皆可以，我们利用图像派生出势(force)。这个势值在物体内应该很高，在十分接近物体的边缘时，应该很低（我们是从物体内部去得到它的边缘的）。图像的梯度表现了它的边缘所在。 
  
势可以是上面梯度图像的逆，也可以是它的高斯部分，如下： 
 
λ和σ是控制惩罚因子的参数。下图是使用梯度图像来演化轮廓： 
  

到此为止，你应该了解，水平集方法是关于演化表面的函数，而不是演化曲线的函数~~（尽管我们的目的是得到轮廓） 
这也是这个方法的优雅魅力之处！