[图像处理]水平集(Level set)算法实现思路(简化)_水平集算法

[图像处理]水平集(Level set)算法实现思路(简化)

创建时间：2020年6月22日
修改时间：2021年6月12日

一、水平集介绍

       图像分割的首要任务是确定图像的轮廓。图像内物体的轮廓，其实际上是一条二维平面内的闭合曲线，而基于对这一数学事实的思考，研究者们开始探究从更高维度的思考角度审视图像分割的算法实现，水平集方法(Level Set Method，LSM)应运而生。
       LSM不直接对图像轮廓进行检测操作，而是将图像轮廓嵌入到更高维度的函数，此时，当前的图像的二维轮廓即为这个高纬度曲面函数的零水平集。这个高维函数称为水平集函数(Level Set Function, LSF)，图像轮廓被认为是LSF的零水平集，即当前轮廓(闭合曲线L)满足
L = {   ( x , y )   ∣   ϕ ( x , y ) = 0   } (1) L=\{~(x,y)~|~\phi(x,y)=0~\} \tag{1} L={ (x,y) ∣ ϕ(x,y)=0 }(1)其中， $(x,y)$ 为该图像轮廓的坐标点集。可以说，闭合曲线L，即当前图像轮廓由LSF的零水平集隐含的表达。LSM将图像轮廓，即二维闭合曲线描述为三维空间内LSF的演化。
       李纯明等人基于对LSM的思考，提出了一种基于距离正则化的不需要对水平集函数进行重新初始化的LSF模型，称为DRLSE模型。又由于对目标图像实现分割需要人为参与，帮助计算机进行选择，因此在本次实验报告中，我还设计了人机交互板块，对用户输入的特定信息进行分析并予以反馈。以下对DRLSE模型的主要迭代公式，以及LSF函数初始化的方法和图像轮廓的求解方式进行简单介绍。

二、DRLSE模型构建

(一) 水平集算法的主要迭代公式

       DRLSE模型是由能量公式 $\mathcal{E}(\varphi )=\mu R_p(\varphi )+{\mathcal{E}}_{ext}(\varphi )$ 演化而来，其中 ${\mathcal{E}}_{ext}(\varphi )$ 根据LSM的应用领域的不同而有不同的定义。在图像分割领域，DRLSE模型将其定义为
$$\mathcal{E}(\varphi )=\mu R_p(\varphi )+\lambda {\mathcal{L}}_g(\varphi )+\alpha {\mathcal{A}}_g(\varphi )\text{\hspace{1em}(2)}$$其中 $\mu >0$，$\lambda >0$ 并且$\alpha$为常数。在公式(2)中，其第二项 ${\mathcal{L}}_g$ 表示函数 $g$ (公式(11)) 沿着零水平集轮廓的线积分，当LSF的零水平集位于物体的边界时，该线积分最小。公式第三项 ${\mathcal{A}}_g$ 描述为当前零水平集，即二维闭合曲线内区域的加权面积。如果初始轮廓放置在物体外部，则 $\alpha >0$ ，以便在水平集演化过程中零水平轮廓会缩小；如果初始轮廓放置在物体内部，则 $\alpha <0$ 以展开轮廓。当零水平集靠近物体边界时，第三项控制曲线扩展或收缩的速率将减小，使得零水平集逐渐稳定下来（文献使用的参数为$\mu =0.2,\lambda =5,\alpha =-3$）。公式(3) 即对 公式(2) 的离散化。
$$\frac{\partial \varphi }{\partial t}=\mu \mathrm{div}\left(d_p(|\nabla \varphi |)\nabla \varphi \right)+\lambda {\delta }_{\epsilon }(\varphi )\mathrm{div}\left(g\frac{\nabla \varphi }{|\nabla \varphi |}\right)+\alpha g{\delta }_{\epsilon }(\varphi )\text{\hspace{1em}(3)}$$其中 $\delta (\varphi )$ 为狄拉克函数，$x$ 为空间变量。由于增加了距离正则化项，公式(3)中的所有空间导数都可以用中心差分格式来离散化。中心差分格式比传统水平集公式中常用的一阶迎风格式更精确、更有效。

       DRLSE模型等号右侧公式可以代替为
$$L\left({\varphi }_{i,j}^k\right)=\mu \mathrm{div}\left(d_p(|\nabla \varphi |)\nabla \varphi \right)+\lambda {\delta }_{\epsilon }(\varphi )\mathrm{div}\left(g\frac{\nabla \varphi }{|\nabla \varphi |}\right)+\alpha g{\delta }_{\epsilon }(\varphi )\text{\hspace{1em}(4)}$$在检查文献代码中，发现 公式(4) 中对散度的求解，即 $div(\cdot )$ 算子的计算方式与 公式(4) 描述有一些不同。借鉴文献代码，

（1）对第一项的求解实际上是按照
$$\mathrm{div}\left(d_p(|\nabla \varphi |)\nabla \varphi \right)=\mathrm{divergence}\left(dp(\nabla \varphi )\cdot {\varphi }_x-{\varphi }_x,dp(\nabla \varphi )\cdot {\varphi }_v-{\varphi }_v\right)+4\ast \mathrm{del}2(\varphi )\text{\hspace{1em}(5)}$$的方式来实现的。其中 $divergence$ 和 $del2$ 都是Matlab内置的函数。
（2）对第二项的求解实际上是按照
$$\mathrm{div}\left(g\frac{\nabla \varphi }{|\nabla \varphi |}\right)=g_x\cdot \frac{{\varphi }_x}{\nabla \varphi }+g_y\cdot \frac{{\varphi }_y}{\nabla \varphi }+g\cdot divergence\left(g_x,g_y\right)\text{\hspace{1em}(6)}$$代码披露如下

% Gradfx,Gradfy 为LSF的偏导数矩阵  
fa = divergence(dps.*Gradfx-Gradfx,dps.*Gradfy-Gradfy)+4*del2(faiLSF);  % 公式(5) 
% 需要防止分母出现为0的情况，否则将使得分式没意义出现NaN  
gx = Gradfx./G;   
gy = Gradfy./G;  
% g_Gfx,g_Gfy 为函数g的偏导数矩阵  
divMatrix = (g_Gfx.*gx+g_Gfy.*gy) + g.*divergence(gx,gy); % 公式(6)
1
2
3
4
5
6
7

函数左侧的偏微分方程可以通过有限差分格式来实现，即
$$\frac{\partial \varphi }{\partial t}=\left({\varphi }_{i,j}^{k+1}-{\varphi }_{i,j}^k\right)/\Delta t\text{\hspace{1em}(7)}$$如此一来，我们解决了离散公式的偏微分方程的求解。于是对于DRLSE模型则有迭代公式
$${\varphi }_{i,j}^{k+1}={\varphi }_{i,j}^k+\Delta tL\left({\varphi }_{i,j}^k\right),k=0,1,2,\dots \text{\hspace{1em}(8)}$$其中，$k$ 表示当前迭代的次数，$\Delta t$ 表示时间步长。我们用在笛卡尔坐标网格上定义的LSF，即 ${\varphi }_{i,j}$ 来表示水平集。对于 公式(8) 的理解，可以认为是：
$$之后的位置=当前位置+时间间隔\times 速度$$

(二) 水平集函数初始化

      通常使用二维阶跃函数作为初始LSF，因为它可以极其有效地生成。二维阶跃函数定义为
ϕ 0 ( x ) = { − c 0 ,  if  x ∈ R 0 c 0 ,  otherwise  (9) \phi_{0}(\mathbf{x})=\left\{

$$\begin{array}{ll} -c_{0}, & \text { if } \mathbf{x} \in R_{0} \\ c_{0}, & \text { otherwise } \end{array}$$ \right. \tag{9} ϕ0​(x)={−c0​,c0​,​ if x∈R0​ otherwise ​(9)其中$c_0$为一个正的常数，而$R_0$位于图像域当中。对于一张输入图像$I$，我们可以自定义初始的$R_0$范围。$x$是一个空间位置参量，实际上是图像内网格点的坐标值$(i,j)$。难点是判断该网格点是否在$R_0$内。然而使用二位阶跃函数的好处在于，当$(i,j)$位于$R_0$内时，LSF函数值为负；otherwise，为正。这样我们就很好判断当前点的位置了。
      此外，该区域$R_0$有时可以通过简单有效的初步分割步骤，如阈值分割，使其接近要分割的区域。因此，只需少量的迭代就可以将零水平集从的$R_0$边界移动到所需的图像内的对象的边界。

(三) 零水平集的计算-窄带构建

      窄带的构建过程，实际上是确定LSF零水平集位置的过程。在一维曲线当中，曲线某个区间存在零点的充要条件是

1）曲线连续；
2）区间两个端点的函数值的符号相反。

确定LSF的零交叉点也可以用相似的方法来确定。我们仍然使用掩模来完成。在当前处理像素点的领域$N_{i,j}^{(r)}$内，如果${\varphi }_{i-1,j}$和${\varphi }_{i+1,j}$是相反的符号，或者${\varphi }_{i,j-1}$和${\varphi }_{i,j+1}$是相反的符号，那么当前计算的网格点$(i,j)$称之为一个零交叉点。LSF的所有零交叉点的集合用$Z$来表示。然后，我们构造窄带为
$$B_r=\bigcup _{(i,j)\in Z}{N}_{i,j}^{(r)}\text{\hspace{1em}(10)}$$其中$N_{i,j}^{(r)}$是当前网格点的一个$(2r+1)\times (2r+1)$尺寸的邻域，这个就是掩模的尺寸大小。该过程可以称之为遍历过程。DRLSE的窄带实现包括以下步骤:

1. 初始化。确定一个LSF的零水平集函数。我们需要通过人机交互来选择初始零水平集轮廓，这在Section II当中已经说明。然后，依据此可以构造初始窄带$B_r^0={\bigcup }_{⟨i,j⟩\in Z^0}{N}_{i,j}^{(r)}$，其中$Z^0$是零水平集函数中零交叉点的集合。于是，我们应该开辟出一个新的存储区域，来存储这个集合。
2. 更新LSF函数。使用${\varphi }_{i,j}^{k+1}={\varphi }_{i,j}^k+\tau L\left({\varphi }_{i,j}^k\right)$，即公式(4)更新描述的窄带$B_r^k$。这一步需要遍历整个图像区域，完成对整个图像区域的新零水平集计算。
3. 更新窄带。确定窄带$B_r^k$内的${\varphi }_{i,j}^{k+1}$的所有零交叉点，将这些零交叉点的集合由$Z^{k+1}$来描述，然后通过赋值更新窄带的代数。这个对新一代的窄带的赋值关系可以描述为$B_r^{k+1}={\bigcup }_{(i,j)\in Z^{k+1}}{N}_{i,j}^{(r)}$。
4. 为窄带上的新像素赋值。位于窄带$B_r^{k+1}$内而不在窄带$B_r^k$内的像素点，若${\varphi }_{i,j}^{k+1}$>0则该点的像素值设置为h，否则设置为-h。h是一个常数，它可以设置为h=r+1并作为一个默认值。
5. 确定迭代的终止准则。如果连续迭代的零交叉点停止变化，或者超过规定的最大迭代次数，则停止迭代；否则，转到步骤2。每一次循环是在求出${\varphi }_{i,j}^{k+1}$之后。

文献代码中，对LSF零水平集位置的确定通过contour函数来实现，代码如下

% 在原始图像上绘制出此时的初始LSF的零水平集  
level = [0,0];  
contour(initial_faiLSF,level,'r','LineWidth',2);% contour思路借鉴
1
2
3

      DRLSE模型迭代的终止准则为，当连续迭代的零交叉点停止变化，或者超过规定的最大迭代次数，则模型停止迭代。

三、必要公式补充

​ 公式(3)中函数 $g$ 的表达式
$$g\triangleq \frac{1}{1+\nabla {\left|G_{\sigma }\ast I\right|}^2}\text{\hspace{1em}(11)}$$其中 $G_{\sigma }$ 是一个高斯核函数，$\sigma$ 为一个标准差。公式(11) 用来平滑图像当中的噪声，并求取平滑后的图像梯度场模的平方，再按照 公式(11) 的描述求解函数 $g$。但当前计算点位于图像边界时，此时梯度值较大， $g$ 值小，于是函数 $g$ 通常在对象边界处取较小的值。此外 $dp$ 函数被定义为
$$d_p(s)\triangleq \frac{p^{\prime }(s)}{s}\text{\hspace{1em}(12)}$$其中函数 $p$ 为势能函数，依据文献采用double-well即双井式势能函数 $p=p_2$ 来设计，$p_2$ 表达式为
p 2 ( s ) = { 1 ( 2 π ) 2 ( 1 − cos ⁡ ( 2 π s ) ) ,  if  s ≤ 1 1 2 ( s − 1 ) 2 ,  if  s ≥ 1 (13) p_{2}(s)=\left\{

$$\begin{array}{ll} \frac{1}{(2 \pi)^{2}}(1-\cos (2 \pi s)), & \text { if } s \leq 1 \\ \frac{1}{2}(s-1)^{2}, & \text { if } s \geq 1 \end{array}$$ \right. \tag{13} p2​(s)={(2π)21​(1−cos(2πs)),21​(s−1)2,​ if s≤1 if s≥1​(13)和一阶微分的表达式
$$p_2^{\prime }(s)=\{\begin{array}{ll}\frac{1}{2\pi }\sin (2\pi s),& \text{ if }s\le 1\\ s-1,& \text{ if }s\ge 1\end{array}\text{\hspace{1em}(14)}$$对于 $\delta (\varphi )$ 和 $H$ 分别为狄拉克函数和Heaviside函数，他们由两个函数 ${\delta }_{\epsilon }$ 和 $H_{\epsilon }$ 光滑近似，如
$${\delta }_{\epsilon }(x)=\{\begin{array}{ll}\frac{1}{2\epsilon }\left[1+\cos \left(\frac{\pi x}{\epsilon }\right)\right],& |x|\le \epsilon \\ 0,& |x|>\epsilon \end{array}\text{\hspace{1em}(15)}$$和
$$H_{\epsilon }(x)=\{\begin{array}{ll}\frac{1}{2}\left(1+\frac{x}{\epsilon }+\frac{1}{\pi }\sin \left(\frac{\pi x}{\epsilon }\right)\right),& |x|\le \epsilon \\ 1,& x>\epsilon \\ 0,& x<-\epsilon \end{array}\text{\hspace{1em}(16)}$$注意到 ${\delta }_{\epsilon }$ 为 $H_{\epsilon }$ 的导数。参数 $\epsilon$ 通常设置为1.5。