深度学习（4）——强度非均匀性下图像分割的水平集方法及其在MRI中的应用（上）

强度非均匀性下图像分割的水平集方法及其在MRI中的应用

---

摘要

强度不均匀性经常出现在真实世界的图像中，这对图像分割提出了相当大的挑战。最广泛使用的图像分割算法是基于区域并且通常依赖于感兴趣区域中图像强度的均匀性，这通常由于强度不均匀性而不能提供准确的分割结果。本文提出了一种新的基于区域的图像分割方法，能够处理分割中的强度不均匀性。首先，基于具有强度不均匀性的图像模型，我们导出图像强度的局部强度聚类特性，并且为每个点的邻域中的图像强度定义局部聚类准则函数。然后，相对于邻域中心对该局部聚类标准函数进行积分，以给出图像分割的全局标准。在水平集公式中，该标准根据表示图像域分区的水平集函数和考虑图像的强度不均匀性的偏置场来定义能量。因此，通过最小化该能量，我们的方法能够同时分割图像并估计偏置场，并且估计的偏置场可以用于强度不均匀性校正（或偏差校正）。我们的方法已经在合成图像和各种模态的真实图像上得到验证，在存在强度不均匀性的情况下具有期望的性能。实验表明，与众所周知的分段光滑模型相比，我们的方法对初始化更加鲁棒，更快速，更准确。作为一种应用，我们的方法已被用于磁共振（MR）图像的分割和偏差校正，具有良好的结果。

1、引言

由于各种因素，例如照明的空间变化和成像装置的缺陷，强度不均匀性经常发生在真实世界的图像中，这使图像处理和计算机视觉中的许多问题复杂化。特别地，由于要分割的区域中的强度范围之间的重叠，对于具有强度不均匀性的图像，图像分割可能相当困难。这使得不可能基于像素强度识别这些区域。那些广泛使用的图像分割算法[4]，[17]，[18]，[23]通常依赖于强度均匀性，因此不适用于强度不均匀性的图像。 通常，强度不均匀性在图像分割中一直是一个具有挑战性的难题。

水平集方法最初用作跟踪界面和形状的数值技术[14]，在过去十年中越来越多地应用于图像分割[2]，[4]，[5]，[8] - [12]，[15]。在水平集方法中，轮廓或曲面表示为更高维函数的零水平集，通常称为水平集函数。 利用水平集表示，可以基于完善的数学理论，包括变化微积分和偏微分方程（PDE），以原理方式公式化和解决图像分割问题。水平集方法的一个优点是涉及曲线和曲面的数值计算可以在固定的笛卡尔网格上执行，而不必参数化这些对象。此外，水平集方法能够表示具有复杂拓扑的轮廓/表面并以自然方式改变其拓扑。

用于图像分割的现有水平集方法可以分为两大类：基于区域的模型[4]，[10]，[17]，[18]，[20]，[22]和基于边缘的模型[3] ，[7]，[8]，[12]，[21]。基于区域的模型旨在通过使用特定区域描述符来指导每个感兴趣区域以指导活动轮廓的运动。然而，很难为具有强度不均匀性的图像定义区域描述符。大多数基于区域的模型[4]，[16] - [18]都是基于强度同质性的假设。典型的例子是[4]，[16]-[18]中提出的分段常数（PC）模型。在[20]，[22]中，水平集方法是基于Mumford和Shah[13]最初提出的一般分段光滑（PS）公式提出的。这些方法不假设图像强度的均匀性，因此能够分割具有强度不均匀性的图像。然而，这些方法在计算上太昂贵并且对轮廓的初始化非常敏感[10]，这极大地限制了它们的实用性。基于边缘的模型使用边缘信息进行图像分割。这些模型不假设图像强度的均匀性，因此可以应用于具有强度不均匀性的图像。然而，这种类型的方法通常对初始条件非常敏感，并且经常在具有弱物体边界的图像中遭受严重的边界泄漏问题。

在本文中，我们提出了一种新的基于区域的图像分割方法。从普遍接受的具有强度不均匀性的图像模型中，我们推导出局部强度聚类特性，因此为每个点的邻域中的强度定义局部聚类准则函数。该局部聚类标准被整合在邻域中心上以定义能量函数，其被转换为水平集公式。通过水平集演化的交错过程和偏置场的估计来实现对该能量的最小化。作为一种重要的应用，我们的方法可用于磁共振（MR）图像的分割和偏差校正。请注意，本文是我们在会议论文[9]中提出的初步工作的扩展版本。

本文的结构如下。我们首先回顾了第二节中两个著名的基于区域的图像分割模型。在第三节中，我们提出了一种用于图像分割和偏置场估计的能量最小化框架，然后在第四节中将其转换为能量最小化的水平集公式。实验结果在第五节中给出，接着讨论了我们的模型与分段光滑的Mumford-Shah和分段常数Chan-Vese模型在第六部分之间的关系。本文结论在第七节。

2.背景

设$\Omega$为图像域，$I:\Omega \to R$为灰度图像。在[13]中，图像$I$的分割是通过找到一个轮廓$C$来实现的，它将图像域分割为不相交的区域${\Omega }_1,\cdot \cdot ,{\Omega }_N$，以及一个分段光滑函数$u$，该函数近似于图像$I$，在每个区域${\Omega }_i$都是光滑的。这可以被表述为一个最小化Mumford-Shah函数的问题 。

$$F^{MS}(u,C)={\int }_{\Omega }(I-u{)}^2\mathrm{d}x+u{\int }_{\Omega \backslash C}\mid \nabla u{\mid }^2\mathrm{d}x+\nu \mid C\mid \text{\hspace{1em}(1)}$$

$\mid C\mid$是轮廓的长度。在上述公式的等号右侧，第一项是数据项，它强制$u$近似图像$I$，第二项是光滑项，它强制$u$在由轮廓$C$分隔的每个区域内光滑。第三项用于正规化轮廓$C$。

在$\Omega$内，${\Omega }_1,\cdot \cdot ,{\Omega }_N$是被轮廓$C$分隔的区域，即$\Omega \backslash C=U_{i=1}^N{\Omega }_i$。所以，能量$F^{MS}(u,C)$可以表示成

$$F^{MS}(u_1,\cdot \cdot ,u_N,{\Omega }_1,\cdot \cdot ,{\Omega }_N)=\sum _{i=1}^N\left({\int }_{{\Omega }_i}(I-u_i{)}^2\mathrm{d}x+u{\int }_{{\Omega }_i}\mid \nabla u{\mid }^2\mathrm{d}x+\nu \mid C_i\mid \right)$$

其中$u_i$是在区域${\Omega }_i$上定义的光滑函数。 旨在最小化这种能量的方法称为分段光滑（PS）模型。 在[20]，[22]中，水平集方法被提出作为用于图像分割的PS模型。

能量$F^{MS}$的变量包括N个不同的函数$u_1,\cdot \cdot ,u_N$。${\Omega }_i$中每个函数$u_i$的平滑度必须通过在函数$F^{MS}$中施加平滑项$u{\int }_{{\Omega }_i}\mid \nabla u{\mid }^2\mathrm{d}x$来确保。为了最小化能量$F^{MS}$，引入了用于求解与相应平滑项相关联的函数$u_1,\cdot \cdot ,u_N$的N个PDF，并且必须在轮廓$C$或区域${\Omega }_1,\cdot \cdot ,{\Omega }_N$的演化中的每个时间步骤处求解。该过程在计算上非常消耗资源。此外，PS模型对轮廓或区域的初始化敏感。从第V-A部分的一些实验结果可以看出这些困难。

在变分水平集公式[4]中，Chan和Vese将Mumford-Shah函数简化为以下能量公式：
$$\begin{array}{cc}{F}^{CV}(\varphi ,c_1,c_2)=& {\int }_{\Omega }\mid I(x)-c_1{\mid }^{2}H(\varphi (x))\mathrm{d}x\\ & +{\int }_{\Omega }\mid I(x)-c_2{\mid }^2(1-H(\varphi (x)))\mathrm{d}x\\ & +\nu {\int }_{\Omega }\mid \nabla (\varphi (x))\mid \mathrm{d}x\end{array}\text{\hspace{1em}(2)}$$
$H$为阶跃函数（Heaviside function），$\varphi$是水平集函数，零水平轮廓 C = { x : ϕ ( x ) = 0 } C = \{ x : \phi(x)=0 \} C={x:ϕ(x)=0}将图像域$\Omega$分割成两个不相交的区域 Ω 1 = { x : ϕ ( x ) &gt; 0 } \Omega_1 = \{ x : \phi(x)&gt;0 \} Ω1​={x:ϕ(x)>0}和 Ω 2 = { x : ϕ ( x ) &lt; 0 } \Omega_2 = \{ x : \phi(x)&lt;0 \} Ω2​={x:ϕ(x)<0}。上述公式中，前两项是数据拟合项，第三项是权重，其中$\nu >0$，用于正则化零水平轮廓。因此，通过找到水平集函数$\varphi$和使能量$F^{CV}$最小化的常数$c_1$和$c_2$来实现图像分割。该模型是分段常数（PC）模型，因为它假设图像$I$可以分别通过区域${\Omega }_1$和${\Omega }_2$中的常数$c_1$和$c_2$来近似。

3.联合分割和偏差场估计的变分框架

A.图像模型和问题公式

为了处理图像分割中的强度不均匀性，我们基于描述真实世界图像的组成的图像模型来制定我们的方法，其中强度不均匀性归因于图像的分量。在本文中，我们考虑以下强度不均匀性的乘法模型。根据各种模态（例如摄像机和MRI）的成像物理特性，可以将观察到的图像建模为:
$$I=bJ+n\text{\hspace{1em}(3)}$$
其中$J$是真实图像，$b$是考虑强度不均匀性的分量，$n$是加性噪声。分量$b$被称为偏置场（或阴影图像）。真实图像$J$测量被成像物体的固有物理特性，因此假设它是分段（近似）恒定的。假设偏置场$b$缓慢变化。加性噪声$n$可以假设为零均值高斯噪声。

在本文中，我们将图像$I$视为在连续域$\Omega$上定义的函数$I:\Omega \to R$。关于真实图像$J$和偏差场$b$的假设可以更具体地说明如下：

• (A1)、偏置场$b$缓慢变化，这意味着$b$可以很好地近似于图像域中每个点的邻域中的常数。
• (A2)、真实图像$J$大致分别在不相交的区域${\Omega }_1,\cdot \cdot ,{\Omega }_N$中分别取N个不同的常数值$c_1,\cdot \cdot ,c_N$,其中， { Ω i } i = 1 N \{ \Omega_i \}_{i=1}^N {Ωi​}i=1N​,即$\Omega ={\bigcup }_{i=1}^N{\Omega }_i$,当$i̸=j$时，${\Omega }_i\bigcap {\Omega }_j≠0$。

基于公式$I=bJ+n$和(A1)、(A2)，我们提出了一种估计区域， { Ω i } i = 1 N \{ \Omega_i \}_{i=1}^N {Ωi​}i=1N​，常数 { c i } i = 1 N \{ c_i \}_{i=1}^N {ci​}i=1N​和偏置场$b$的方法。所获得的估计值分别由 { Ω ^ i } i = 1 N \{ \hat{\Omega}_i \}_{i=1}^N {Ω^i​}i=1N​，常数 { c ^ i } i = 1 N \{ \hat{c}_i \}_{i=1}^N {c^i​}i=1N​和偏置场$b$表示。所获得的偏置场$\hat{b}$应该缓慢变化，并且区域${\hat{\Omega }}_1,\cdot \cdot ,{\hat{\Omega }}_N$应该满足一定的规律性，以避免由图像噪声引起的伪分割结果。我们将基于上述图像模型和假设A1和A2来定义用于寻找此类估计的标准。该标准将根据区域${\Omega }_i$，常数$c_i$和函数$b$来定义，作为变分框架中的能量，其被最小化以找到最佳区域 { Ω ^ i } i = 1 N \{ \hat{\Omega}_i \}_{i=1}^N {Ω^i​}i=1N​，常数 { c ^ i } i = 1 N \{ \hat{c}_i \}_{i=1}^N {c^i​}i=1N​和偏置场$\hat{b}$，并同时完成图像分割和偏置场估计。

B.局部强度聚类属性

基于区域的图像分割方法通常依赖于要分割的每个区域中的强度的特定区域描述符（例如，强度均值或高斯分布）。然而，对于具有强度不均匀性的图像，难以给出这样的区域描述符。此外，强度不均匀性经常导致区域${\Omega }_1,\cdot \cdot ,{\Omega }_N$中强度分布之间的重叠。因此，不可能基于像素强度直接分割这些区域。然而，局部强度的性质很简单，可以有效地利用我们的图像分割方法的制定，同时估计偏置场。

基于公式$I=bJ+n$和(A1)、(A2)，我们能够得到局部强度的有用特性，这被称为局部强度聚类特性，如下所述和证明的那样。 具体而言，我们考虑一个圆形邻域，其半径为$\rho$，以每个点$y\in \Omega$为中心,由 O y ≜ { x : ∣ x − y ∣ ≤ ρ } O_y \triangleq \{ x:\mid x-y\mid \leq \rho \} Oy​≜{x:∣x−y∣≤ρ}定义。整个域$\Omega$的分区 { Ω i } i = 1 N \{ \Omega_i \}_{i=1}^N {Ωi​}i=1N​引起邻域$O_y$的分区，即 { O y ⋂ Ω i } i = 1 N \{ O_y\bigcap\Omega_i \}_{i=1}^N {Oy​⋂Ωi​}i=1N​ 形成 $O_y$的分区。 对于缓慢变化的偏置场$b$，圆形邻域$O_y$中所有$x$的值，$b(x)$接近$b(y)$，即
$$b(x)\approx b(y)\text{for x ∈ Oy}\text{\hspace{1em}(4)}$$
。因此，每个子区域$O_y\bigcap {\Omega }_i$中的强度$b(x)J(x)$接近常数$b(y)c_i$,即
$$b(x)J(x)\approx b(y)c_i\text{ for x∈Oy⋂Ωi}\text{\hspace{1em}(5)}$$

然后，鉴于公式(3)中的图像模型，我们有:$$I(x)\approx b(y)c_i+n(x)\text{for x∈Oy⋂Ωi}$$其中$n(x)$是加性零均值高斯噪声。因此，集合的强度:

I y i = { I ( x ) : x ∈ O y ⋂ Ω i } I_y^i = \{ I(x):x\in O_y \bigcap \Omega_i\} Iyi​={I(x):x∈Oy​⋂Ωi​}

形成具有聚类中心$m_i\approx b(y)c_i$的聚类，其可以被认为是从具有平均$m_i$的高斯分布中抽取的样本。显然，$N$个群集$I_y^1,\cdot \cdot \cdot ,I_y^N$是完全分离的，具有不同的聚类中心$m_i\approx b(y)c_i,i=1,\cdot \cdot \cdot ,N$（因为常数$c_1,\cdot \cdot \cdot ,c_N$是不同的，并且假设高斯噪声的方差$n$相对较小）。该局部强度聚类特性用于制定所提出的用于图像分割和偏置场估计的方法，如下所述。

C. 能量配方

上述局部强度聚类特性表明邻域$O_y$中的强度可以被分类为具有中心$m_i\approx b(y)c_i,i=1,\cdot \cdot \cdot ,N$的$N$个聚类。这允许我们应用标准的K均值聚类来对这些局部强度进行分类。具体来说，对于邻域$O_y$中的强度$I(x)$，K均值算法是一个迭代过程，以最小化聚类标准[19]，它可以连续的形式写成:
$$F_y=\sum _{i=1}^N{\int }_{O_y}\mid I(x)-m_i{\mid }^2{u}_i(x)\mathrm{d}x\text{\hspace{1em}(6)}$$

其中$m_i$是第$i$个聚类的聚类中心，$u_i$是要确定的区域${\Omega }_i$的隶属函数，即$$u_i(x)=1\text{, x∈Ωi}$$和 $$u_i(x)=0,\text{ x̸∈Ωi}$$由于$u_i$是区域${\Omega }_i$的隶属函数，我们可以将$F$重写为
$$F_y=\sum _{i=1}^N{\int }_{{\Omega }_i\bigcup O_y}\mid I(x)-m_i{\mid }^2\mathrm{d}x\text{\hspace{1em}(7)}$$

鉴于上式中的聚类准则和$m_i\approx b(y)c_i$对聚类中心的近似，我们定义了一个聚类标准，用于将$O_y$中的强度分类为:
$${\epsilon }_y=\sum _{i=1}^N{\int }_{{\Omega }_i\bigcup O_y}K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}x\text{\hspace{1em}(8)}$$
，其中$K(y-x)$作为非负窗函数引入，也称为核函数，例如$$K(y-x)=0\text{,x̸∈Oy}$$
使用窗口功能，可以将聚类标准函数${\epsilon }_y$重写为:
$${\epsilon }_y=\sum _{i=1}^N{\int }_{{\Omega }_i}K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}x\text{\hspace{1em}(9)}$$
该局部聚类准则函数是我们方法的基本元素。

局部聚类标准函数${\epsilon }_y$评估由$O_y$的分区 { O y ⋂ Ω i } i = 1 N \{O_y \bigcap\Omega_i\}_{i=1}^N {Oy​⋂Ωi​}i=1N​给出的邻域$O_y$中的强度的分类.${\epsilon }_y$的值越小，分类越好。当然，我们将整个域$\Omega$的最佳分区 { Ω i } i = 1 N \{\Omega_i\}_{i=1}^N {Ωi​}i=1N​定义为使得局部聚类准则函数${\epsilon }_y$对于$\Omega$中的所有$y$最小化。因此，我们需要联合最小化$\Omega$中所有$y$的${\epsilon }_y$.这可以通过在图像域$\Omega$上最小化${\epsilon }_y$相对于$y$的积分来实现。因此，我们定义能量$\epsilon \triangleq \int {\epsilon }_y\mathrm{d}y$,即：
$$\epsilon \triangleq \int \left(\sum _{i=1}^N{\int }_{{\Omega }_i}K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}x\right)\mathrm{d}y\text{\hspace{1em}(10)}$$
在本文中，如果积分在整个域$\Omega$上，我们省略积分符号的下标中的域$\Omega$（如上面的第一个积分中）。图像分割和偏置场估计可以通过最小化这个能量来执行。区域${\Omega }_1,\cdot \cdot \cdot ,{\Omega }_N$，常数$c_1,\cdot \cdot \cdot ,c_N$和偏置区域$b$。

核函数K的选择是灵活的。例如，它可以是截断的均匀函数，定义为$$K(u)=a\text{, ∣u∣≤ρ}$$和$$K(u)=0\text{, ∣u∣≥ρ}$$，其中$a$是正常数，使得$\int K(u)=1$.在本文中，核函数$K$被选择为由此定义的截断高斯函数。
$$K(u)=\{\begin{array}{ll}\frac{1}{a}{e}^{-\mid u{\mid }^2/2{\sigma }^2},& \text{for∣u∣≤ρ}\\ 0,& \text{otherwise}\end{array}\text{\hspace{1em}(11)}$$
其中$a$是归一化常数，使得$\int K(u)=1$，$\sigma$是高斯函数的标准偏差（或尺度参数），并且$\rho$是邻域$O_y$的半径。

注意，应该根据强度不均匀性的程度适当地选择邻域$O_y$的半径$\rho$。对于更局部化的强度不均匀性，偏置场$b$变化更快，因此上式中的近似仅在较小的邻域中有效。在这种情况下，应该使用较小的$\rho$作为邻域$O_y$的半径，对于上式中的截断高斯函数，缩放参数$\sigma$也应该更小。

4.水平集制定和能量最小化

我们在（10）中提出的能量${\epsilon }_y$用区域${\Omega }_1,\cdot \cdot \cdot ,{\Omega }_N$表示。很难从这种表达式推导出能量最小化问题的解决方案。 在本节中，通过用多个水平集函数表示不相交区域${\Omega }_1,\cdot \cdot \cdot ,{\Omega }_N$，将能量${\epsilon }_y$转换为水平集公式，在这些水平集函数上具有正则化项。在水平集配方中，能量最小化可以通过使用完善的变分方法来解决[6]。

在水平集方法中，水平集函数是带有正号和负号的函数，其可以用于将域$\Omega$的分区表示为两个不相交的区域${\Omega }_1$和${\Omega }_2$。设$\varphi :\Omega \to R$为水平集函数，然后其符号定义两个不相交的区域
$${\Omega }_1=x:\varphi (x)>0\text{and}{\Omega }_2=x:\varphi (x)<0\text{\hspace{1em}(12)}$$
形成域$\Omega$的分区。对于N>2的情况，可以使用两个或更多个水平集函数来表示$N$个区域${\Omega }_1,\cdot \cdot \cdot ,{\Omega }_N$。 对于$N=2$和$N>2$的情况，能量${\epsilon }_y$的水平集配方，分别称为两相和多相配方，将在接下来的两个小节中给出。

A.两相水平集配方

我们首先考虑两相情况：图像域$\Omega$被分割成两个不相交的区域${\Omega }_1$和${\Omega }_2$。 在这种情况下，水平集函数$\varphi$用于表示由（12）给出的两个区域${\Omega }_1$和${\Omega }_2$。区域${\Omega }_1$和${\Omega }_2$可以用它们分别由$M_1(\varphi )=H(\varphi )$和$M_2(\varphi )=1-H(\varphi )$定义的隶属函数表示，其中H是Heaviside函数。 因此，对于$N=2$的情况，(10)中的能量可以表示为以下水平集公式：
$$\epsilon =\int \left(\sum _{i=1}^N\int K(y-x)\mid I(x)-b(y)c_i{\mid }^2{M}_i(\varphi (x))\mathrm{d}x\right)\mathrm{d}y\text{\hspace{1em}(13)}$$
通过交换集成的顺序，我们有
$$\epsilon =\int \sum _{i=1}^N\left(\int K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}x\right)M_i(\varphi (x))\mathrm{d}x\text{\hspace{1em}(14)}$$

为方便起见，我们用向量$\vec{c}=(c_1,\cdot \cdot \cdot ,c_N)$表示常量$c_1,\cdot \cdot \cdot ,c_N$。因此，水平集函数$\varphi$，向量$\vec{c}$和偏置字段$b$是能量$\epsilon$的变量，因此可以写为$\epsilon (\varphi ,\vec{c},b)$。从上式开始，我们可以用以下形式重写能量$\epsilon (\varphi ,\vec{c},b)$：
$$\epsilon (\varphi ,\vec{c},b)=\int \sum _{i=1}^N{e}_i(x)M_i(\varphi (x))\mathrm{d}x\text{\hspace{1em}(15)}$$
函数$e_i$定义为：
$$e_i(x)=\int K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}y\text{\hspace{1em}(16)}$$
可以使用以下等效表达式计算函数$e_i$：
$$e_i(x)=I^2{1}_K-2C_{i}I(b\ast K)+c_i^2(b^2\ast K)\text{\hspace{1em}(17)}$$
其中$\ast$是卷积运算，$1_K$是由$1_K(x)=\int K(y-x)\mathrm{d}y$定义的函数，它在除了图像域$\Omega$的边界附近以外的任何地方等于常数1。

上面定义的能量$\epsilon (\varphi ,\vec{c},b)$用作所提出的变分水平集公式的能量中的数据项，其由下式定义：
$$F(\varphi ,\vec{c},b)=\epsilon (\varphi ,\vec{c},b)+\nu L(\varphi )+\mu R_p(\varphi )\text{\hspace{1em}(18)}$$
$L(\varphi )$和$R_p(\varphi )$是如下定义的正则化项。

能量项$L(\varphi )$由下式定义:
$$L(\varphi )=\int \mid ▽H(\varphi )\mid \mathrm{d}x\text{\hspace{1em}(19)}$$
它计算了$\varphi$的零水平轮廓的弧长，因此通过惩罚其弧长来平滑轮廓[4]，[10]。

能量项$R_p(\varphi )$定义为
$$R_p(\varphi )=\int p(\mid ▽\varphi \mid )\mathrm{d}x\text{\hspace{1em}(20)}$$

具有电势（能量密度）函数$p:[0,\infty )\to R$使得所有$s$的$p(s)\ge p(1)$，即$s=1$是$p$的最小点。在本文中，我们使用$p(s)=\frac{1}{2}(s-1{)}^2$定义的势函数$p$。显然，当具有电势$p$时，能量$R_p(\varphi )$在$\mid ▽\varphi \mid =1$时被最小化，这是有符号距离函数的特征，称为有符号距离属性。因此，正则化项R称为距离正则化项，由Li等人引入。 [11]在更一般的变分水平集公式中称为距离正则化水平集演化（DRLSE）公式。读者可以参考[11]了解在DRLSE中维持水平集函数的有符号距离属性的必要性和机制。

通过最小化该能量，我们获得由水平集函数$\varphi$给出的图像分割的结果和偏置场$b$的估计。通过迭代过程实现能量最小化：在每次迭代中，我们最小化关于其每个变量$\varphi$，$\vec{c}$和$b$的能量$F(\varphi ,\vec{c},b)$，给定在先前迭代中更新的其他两个变量。我们给出了关于每个变量的能量最小化的解决方案，如下所述。

1）关于$\varphi$的能量最小化

对于固定的$\vec{c}$和$b$，通过使用标准梯度下降法，即求解梯度流方程，可以实现$F(\varphi ,\vec{c},b)$相对于$\varphi$的最小化
$$\frac{\partial \varphi }{\partial t}=-\frac{\partial F}{\partial \varphi }\text{\hspace{1em}(21)}$$
其中$\frac{\partial F}{\partial \varphi }$是能量$F$的Gateaux导数

通过微分变换[1]，我们可以计算出Gateaux导数$\frac{\partial F}{\partial \varphi }$并将相应的梯度流方程表示为:
$$\begin{array}{cc}\frac{\partial \varphi }{\partial t}=& -\delta (\varphi )(e_1-e_2)+\nu \delta (\varphi )div\left(\frac{\nabla \varphi }{\mid \nabla \varphi \mid }\right)\\ & +\mu div(d_p(\mid \nabla \varphi \mid )\nabla \varphi )\end{array}\text{\hspace{1em}(22)}$$
其中$\nabla$是梯度算子，$div(\cdot )$是散度算子，函数$d_p$定义为
$$d_p(s)\triangleq \frac{p^{{}^{\prime }}(s)}{s}$$

如[11]中所述，实现DRLSE的相同有限差分方案可用于水平集演化（22）。在根据（22）的电平集函数的演变期间，$\epsilon (\varphi ,\vec{c},b)$中的常数$c_1$和$c_2$以及偏置字段$b$分别通过最小化相对于$\vec{c}$和$b$的能量来更新，这将在下面描述。

2）关于$\vec{c}$的能量最小化

对于固定的$\varphi$和$b$，最小化能量$\epsilon (\varphi ,\vec{c},b)$ 的最佳$\vec{c}$，由$\hat{\vec{c}}=(\widehat{c_i},\cdot \cdot \cdot ,\widehat{c_N})$表示，由下式给出：
$$\widehat{c_i}=\frac{\int (b\ast K)I{u}_i\mathrm{d}y}{\int (b^2\ast K)u_i\mathrm{d}y}\text{,i=1,⋅⋅⋅,N}\text{\hspace{1em}(23)}$$
其中，$u_i(y)=M_i(\varphi (y))$

3)关于$b$的能量最小化

对于固定的$\varphi$和$\vec{c}$，最小化能量$\epsilon (\varphi ,\vec{c},b)$的最佳$b$，由$\hat{b}$表示，由下式给出：
$$\hat{b}=\frac{I{J}^{(1)}\ast K}{J^{(2)}\ast K}\text{\hspace{1em}(24)}$$
其中$J^{(1)}={\sum }_{i=1}^N{c}_i{u}_i,J^{(2)}={\sum }_{i=1}^N{c}_i^2{u}_i$。注意，上式中具有核函数$K$的卷积证实了偏置场导出的最优估计$\hat{b}$的缓慢变化的性质。

B.多相水平集配方

对于$n>3$的情况，我们可以使用两个或更多水平集函数${\varphi }_1,\cdot \cdot \cdot {\varphi }_k$定义区域$${\Omega }_i\text{,i=1,⋅⋅⋅,N}$$的$N$个隶属函数$M_i$:

$$M_i({\varphi }_1(y),\cdot \cdot \cdot {\varphi }_k(y))=\{\begin{array}{ll}1,& \text{ y∈Ωi }\\ 0,& \text{else}\end{array}$$
例如， $N=3$,我们使用两个水平集函数 ${\varphi }_1$和 ${\varphi }_2$定义 $M_1({\varphi }_1,{\varphi }_2)=H({\varphi }_1)H({\varphi }_2)$, $M_2({\varphi }_1,{\varphi }_2)=H({\varphi }_1)(1-H({\varphi }_2))$和 $M_3({\varphi }_1,{\varphi }_2)=1-H({\varphi }_1)$,给出我们方法的三阶段水平集。对于四阶情形 $N=4$， $M_i$的定义可以定义为 $M_1({\varphi }_1,{\varphi }_2)=H({\varphi }_1)H({\varphi }_2)$, $M_2({\varphi }_1,{\varphi }_2)=H({\varphi }_1)(1-H({\varphi }_2))$, $M_3({\varphi }_1,{\varphi }_2)=(1-H({\varphi }_1))H({\varphi }_2)$和 $M_4=(1-H({\varphi }_1))(1-H({\varphi }_2))$。

为了简化符号，我们用向量值函数$\Phi =({\varphi }_1,\cdot \cdot \cdot ,{\varphi }_k)$表示这些水平集函数${\varphi }_1,\cdot \cdot \cdot ,{\varphi }_k$。因此，隶属函数$M_i({\varphi }_1(y),\cdot \cdot \cdot ,{\varphi }_k(y))$可以写成$M_i(\Phi )$.公式$\epsilon \triangleq \int \left({\sum }_{i=1}^N{\int }_{{\Omega }_i}K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}x\right)\mathrm{d}y$中的能量$\epsilon$可以转换为多相水平集配方：
$$\epsilon (\Phi ,\vec{c},b)=\int \sum _{i=1}^N{e}_i(x)M_i(\Phi (x))\mathrm{d}x$$
$e_i$在前文公式中已给出：$e_i(x)=\int K(y-x)\mid I(x)-b(y)c_i{\mid }^2\mathrm{d}y$

对于函数$\Phi =({\varphi }_1,\cdot \cdot \cdot ,{\varphi }_k)$,我们定义正则化项$L(\Phi )={\sum }_{j=1}^{k}L({\varphi }_i)$和$R_p(\Phi )={\sum }_{j=1}^k{R}_p({\varphi }_i)$,其中$L(\Phi )$和$R_p(\Phi )$分别由$L(\varphi )=\int \mid ▽H(\varphi )\mid \mathrm{d}x$和$R_p(\varphi )=\int p(\mid ▽\varphi \mid )\mathrm{d}x$定义，用于每个水平集合函数$\varphi$。我们的多相水平集配方中的能量函数$F$由下式定义:
$$F(\Phi ,\vec{c},b)\triangleq \epsilon (\Phi ,\vec{c},b)+R_p(\Phi )\text{\hspace{1em}(25)}$$
可以通过求解以下梯度流动方程来执行上式中的能量$F(\Phi ,\vec{c},b)$相对于变量$\Phi =({\varphi }_1,\cdot \cdot \cdot ,{\varphi }_k)$的最小化:
$$\begin{array}{cc}\frac{\partial {\varphi }_1}{\partial l}=& -{\sum }_{i=1}^N\frac{\partial M_i(\varphi )}{\partial {\varphi }_1}{e}_i+\nu \delta ({\varphi }_1)div\left(\frac{\nabla {\varphi }_1}{\mid \nabla {\varphi }_1\mid }\right)\\ & +\mu div(d_p(\mid \nabla {\varphi }_1\mid )\nabla {\varphi }_1)\\ ⋮& \\ \frac{\partial {\varphi }_k}{\partial t}=& -{\sum }_{i=1}^N\frac{\partial M_i(\varphi )}{\partial {\varphi }_k}{e}_i+\nu \delta ({\varphi }_k)div\left(\frac{\nabla {\varphi }_k}{\mid \nabla {\varphi }_k\mid }\right)\\ & +\mu div(d_p(\mid \nabla {\varphi }_k\mid )\nabla {\varphi }_k)\end{array}\text{\hspace{1em}(26)}$$
能量$\epsilon (\Phi ,\vec{c},b)$的最小化可以通过与两相情况相同的过程来实现。并且很容易证明最小化能量$\epsilon (\Phi ,\vec{c},b)$的最佳$\vec{c}$和$b$由公式（23）和（24）给出，其中
$$u_i=M_i(\Phi )\text{,i=1,⋅⋅⋅,N}$$