A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogen

332 阅读12分钟

水平集

Abstract

在现实世界中,图像的强度不均匀性是经常发生的,这给图像分割带来了很大的挑战。目前应用最广泛的图像分割算法是基于区域的,通常依赖于感兴趣区域内图像强度的均匀性,由于强度的不均匀性,往往无法提供准确的分割结果。提出了一种新的基于区域的图像分割方法,该方法能够有效地处理图像分割过程中的强度不均匀性。首先,基于具有强度不均匀性的图像模型,推导出图像强度的局部强度聚类性质,并定义了每个点邻域内图像强度的局部聚类准则函数。然后将该局部聚类准则函数与邻域中心相结合,给出了一种全局的图像分割准则。在水平集公式中,该判据根据水平集函数定义能量,水平集函数表示图像域的划分和偏置场,偏置场解释图像的强度不均匀性。因此,通过最小化该能量,我们的方法能够同时分割图像并估计出偏场,估计出的偏场可以用于强度不均匀性校正(或偏置校正)。我们的方法已经在各种模式的合成图像和真实图像上得到了验证,在存在强度不均匀性的情况下,性能良好。实验表明,该方法具有较强的初始化鲁棒性、较快的速度和较好的精度。作为一种应用,我们的方法已被用于磁共振(MR)图像的分割和偏置校正,取得了良好的效果。

INTRODUCTION

由于光照的空间变化和成像设备的不完善等因素,真实图像中往往会出现强度不均匀现象,这使得许多在计算机图像处理和计算机视觉的问题变得复杂。特别是对于强度不均匀的图像,由于待分割区域的强度范围存在重叠,分割难度较大。这使得基于像素强度来识别这些区域是不可能的。目前广泛使用的图像分割算法[4]、[17]、[18]、[23]通常依赖于灰度均匀性,因此不适用于灰度不均匀的图像。在图像分割中,强度不均匀性一直是一个具有挑战性的难题。

水平集方法最初是用来跟踪界面和形状[14]的数值技术,在过去的十年中,[2]、[4]、[5]、[8]-[12]、[15]等图像分割中得到了越来越多的应用。在水平集方法中,等值线或曲面表示为高维函数的零水平集,通常称为水平集函数。利用水平集表示,可以根据已有的数学理论,包括变分法和偏微分方程(PDE),有原则地构造和求解图像分割问题。水平集方法的一个优点是,涉及曲线和曲面的数值计算可以在一个固定的笛卡尔网格上进行,而无需对这些对象进行参数化。此外,水平集方法能够表示具有复杂拓扑结构的轮廓/曲面,并以自然的方式改变其拓扑结构。
现有的图像分割层次集方法主要分为两大类:基于区域的模型[4]、[10]、[17]、[18]、[20]、[22]和基于边缘的模型[3]、[7]、[18]0、[12]、[21]。基于区域的模型旨在通过使用特定的区域描述符来指导活动轮廓的运动来识别感兴趣的每个区域。然而,对于具有强度不均匀性的图像,很难定义区域描述符。基于区域的[4]、[16]-[18]模型大多基于强度均匀性假设。一个典型的例子是在[4]、[16]-[18]中提出的分段常数(PC)模型。在[20]、[22]中,水平集方法是基于Mumford和Shah[13]最初提出的一般分段光滑(PS)公式提出的。这些方法不假设图像强度的均匀性,因此能够分割具有强度不均匀性的图像。然而,这些方法计算开销太大,而且对轮廓[10]的初始化非常敏感,这极大地限制了它们的实用价值。基于边缘的模型利用边缘信息进行图像分割。这些模型不假设图像强度的均匀性,因此可以应用于具有强度不均匀性的图像。然而,这类方法一般对初始条件比较敏感,在目标边界较弱的图像中往往存在严重的边界泄漏问题
本文提出了一种新的基于区域的图像分割方法。从具有强度不均匀性的图像的一般模型出发,推导出局部强度聚类的性质,并由此定义了每个点邻域内强度的局部聚类准则函数。该局部聚类准则在邻域中心上进行集成,定义一个能量泛函,并将其转化为一个水平集公式。这种能量的最小化是通过水平集演化和偏置场估计的交错过程实现的。作为一种重要的应用,我们的方法可以用于磁共振(MR)图像的分割和偏置校正。注意,本文是我们在会议文件[9]中提出的初步工作的扩展版本。
本文组织如下。在第二节中,我们首先回顾了两个著名的基于区域的图像分割模型。在第三节中,我们提出了一种用于图像分割和偏场估计的能量最小化框架,然后将其转换为第四节中用于能量最小化的水平集公式。第五节给出了实验结果,然后讨论了我们的模型与第六节中的分段光滑Mumford-Shah模型和分段常数Chan-Vese模型之间的关系。

II. BACKGROUND

设Ω是图像域,I:Ω→R是一个灰度图像。图像I的分割是通过找到一个轮廓C,他将图像域Ω分割成不相交的区域\Omega_{1}, \cdots,\Omega_{N},分段光滑函数u是图像I的近似并且是每个地区\Omega_{i}的内部平滑。这可以表示为一个最小化以下Mumford-Shah函数的问题

\mathcal{F}^{M S}(u, C)=\int_{\Omega}(I-u)^{2} d \mathbf{x}+\mu \int_{\Omega \backslash C}|\nabla u|^{2} d \mathbf{x}+\nu|C|        (1)

|C|在这里是轮廓C的长度,(1)右手边的第一项是数据项,他使得u接近于图片I,第二项是平滑项,他使u在每个被C分割的区域内保持平滑 .引入第三项对轮廓C进行正则化。
\Omega_{1}, \cdots, \Omega_{N}成为\Omega中被C分割,即\Omega \backslash C=\cup_{i=1}^{N} \Omega_{i}。然后,轮廓C可以表示为区域边界的并集,表示为C_{1}, \cdots, C_{N}, i.e. C=\cup_{i=1}^{N} C_{i}因此,上述能量\mathcal{F}^{M S}(u, C)可以等价地写成

\mathcal{F}^{M S}\left(u_{1}, \cdots, u_{N}, \Omega_{1}, \cdots, \Omega_{N}\right)
= \sum_{i=1}^{N}\left(\int_{\Omega_{i}}\left(I-u_{i}\right)^{2} d \mathbf{x}+\mu \int_{\Omega_{i}}\left|\nabla u_{i}\right|^{2} d \mathbf{x}+\nu\left|C_{i}\right|\right)

这里u_{i}是一个定义在\Omega_{i}区域上的光滑函数。旨在最小化这种能量的方法称为分段光滑(PS)模型。在[20]、[22]中,提出了层次集方法作为图像分割的PS模型。 能量的变量\mathcal{F}^{M S}包括N个不同的函数u_{1}, \cdots, u_{N},在\mathcal{F}^{M S}中每个函数u_{i}\Omega_{i}中的平滑,必须通过引入平滑项\mu \int_{\Omega_{i}}\left|\nabla u_{i}\right|^{2} d \mathbf{x}来保证,最小化这种能量,N个 PDE求解函u_{1}, \cdots,u_{N}相应平滑项的引入了平滑项和在轮廓C的演化过程中,每个时间步长都必须求解C或区域\Omega_{1}, \cdots,\Omega_{N}。这个过程在计算上很费劲。此外,PS模型对轮廓C或Ω1,Ωn的初始化很敏感。这些困难可以从V-A部分的一些实验结果中看出。
在变分水平集公式[4]中,Chan和Vese将Mumford-Shah函数简化为:

\begin{aligned} \mathcal{F}^{C V}\left(\phi, c_{1}, c_{2}\right)=& \int_{\Omega}\left|I(\mathbf{x})-c_{1}\right|^{2} H(\phi(\mathbf{x})) d \mathbf{x} \\ &+\int_{\Omega}\left|I(\mathbf{x})-c_{2}\right|^{2}(1-H(\phi(\mathbf{x}))) d \mathbf{x} \\ &+\nu \int_{\Omega}|\nabla H(\phi(\mathbf{x}))| d \mathbf{x} \end{aligned}          (2)

其中H为Heaviside函数,ϕ是一个水平集函数,其零水平轮廓为C=\{\mathbf{x} : \phi(\mathbf{x})=0\}划分区域\Omega成不相交区域\Omega_{1}=\{\mathbf{x} : \phi(\mathbf{x})>0\} and \Omega_{2}=\{\mathbf{x} : \phi(\mathbf{x})<0\},(2)中的前两项为数据拟合项 ,第三项带有权重\nu>0,正则化的0水平集轮廓。图像分割是通过找到ϕ的水平集函数和常量c1和c2最小化能量\mathcal{F}^{C V}。这个模型是分段常数(PC)模型,它假定图像I在区域Ω1和Ω2可以用常数 c1和c2,分别近似。 为了解决图像分割中的强度不均匀性问题,提出了一种基于图像模型的图像分割方法。本文考虑了以下强度不均匀性的乘法模型。从物理成像的各种形式(如相机和MRI),观察到的图像I建模为
I=b J+n               (3)
其中。j为真实图像,b为强度不均匀性的分量,n为额外噪声。元素b称为bias field(或阴影图像)。真实图像J测量被成像对象的内在物理特性,因此假定它是分段的(近似的)常数。假设偏场b是缓慢变化的。加性噪声n可以假设为零均值高斯噪声。
本文我们认为图像I定义在连续域\Omega上的函数I:\OmegaR。关于真像J和偏场b的假设可以更具体地表述为:
(A1)偏置场b是缓慢变化的,这意味着b可以很好地被图像域中每个点的邻域内的一个常数近似。
(A2)真实图像J近似取N个不同的常量C_1,…,C_N在不相交区域\Omega_1,...\Omega_N其中{构成图像域的一个划分 基于区域的图像分割方法通常依赖于被分割区域的特定区域描述符(如强度均值或高斯分布)的强度。然而,对于具有强度不均匀性的图像,很难给出这样的区域描述。此外,强度不均匀经常导致不均匀区域\Omega_1,...\Omega_N之间分布的重叠。因此,直接基于像素强度分割这些区域是不可能的。然而,局部强度的性质是简单的,这可以有效地利用在制定我们方法,方法用于图像分割的同时估计偏置场
根据(3)中的图像模型和假设(A1)和(A2),我们提出一个方法来估算区域\left\{\Omega_{i}\right\}_{i=1}^{N},常数\left\{c_{i}\right\}_{i=1}^{N}和偏场b。获得的估计是\left\{\hat{\Omega}_{i}\right\}_{i-1}^{N},常数\left\{\hat{c}_{i}\right\}_{i-1}^{N},和偏差场\hat{b},获得的偏差场\hat{b}应该缓慢变换的,区域\left\{\Omega_{i}\right\}_{i=1}^{N}应该满足一定的正则属性,以避免图像噪声引起的虚假分割结果。我们将根据上述图像模型和假设A1和A2定义一个寻找此类估计的准则。这一标准将定义的地区Ωi常数ci,和函数b上,作为一个能量变分框架,为了找到最小化最优区域\left\{\hat{\Omega}_{i}\right\}_{i-1}^{N},常数\left\{\hat{c}_{i}\right\}_{i-1}^{N},和偏差场\hat{b}。同时实现了图像分割和偏置场估计 。
根据(3)中的图像模型和假设A1和A2,我们可以得到一个有用的局部强度属性,称为局部强度聚类属性,如下所述,并证明其正确性。具体地说,我们考虑一个圆形邻域以ρ为半径,每个点中心\mathbf{y} \in \Omega,定义为\mathcal{O}_{\mathbf{y}} \triangleq\{\mathbf{x} :|\mathbf{x}-\mathbf{y}| \leq \rho\},完整域\Omega的分区\left\{\Omega_{i}\right\}_{i=1}^{N}引入了邻域\mathcal{O}_{\mathbf{y}}的分割,如\left\{\mathcal{O}_{\mathbf{y}} \cap \Omega_{i}\right\}_{i=1}^{N}得到了\mathcal{O}_{\mathbf{y}},对于慢变偏场b,圆邻域\mathcal{O}_{\mathbf{y}}中所有x的值b(x)都接近b(y)

b(\mathbf{x}) \approx b(\mathbf{y}) \quad for \quad \mathbf{x} \in \mathcal{O}_{\mathbf{y}}    (4)

因此,强度b (x) J (x)在每个次区域\mathcal{O}_{\mathbf{y}} \cap \Omega_{i}接近常数b (y) ci

b(\mathbf{x}) J(\mathbf{x}) \approx b(\mathbf{y}) c_{i} \quad for \quad \mathbf{x} \in \mathcal{O}_{\mathbf{y}} \cap \Omega_{i}   (5)

然后,针对(3)中的图像模型,我们得到

I(\mathbf{x}) \approx b(\mathbf{y}) c_{i}+n(\mathbf{x}) \quad for \quad \mathbf{x} \in \mathcal{O}_{\mathbf{y}} \cap \Omega_{i}

\mathbf{I}_{\mathbf{y}}^{i}=\left\{I(\mathbf{x}) : \mathbf{x} \in \mathcal{O}_{\mathbf{y}} \cap \Omega_{i}\right\}形成一个以mi≈b(y)ci为聚类中心的聚类,可以看作是样本来自均值为mi为的高斯分布,显然N个簇\mathbf{I}_{\mathbf{y}}^{1}, \cdots, \mathbf{I}_{\mathbf{y}}^{N}被很好的划分。随着聚类中心的不同,mi≈b(y)ci, i=1,⋯⋯N(由于常数c1,⋯⋯cN是不同的,假设高斯噪声N的方差相对较小)。利用这种局部强度聚类特性,提出了一种图像分割和偏置场估计的方法。
上述局部强度聚类特性表明邻域\mathcal{O}_{\mathbf{y}}中的强度可以划分为N个聚类,中心为mi≈b(y)ci, i=1,⋯⋯N。这允许我们应用标准的K-means聚类来对这些局部强度进行分类。具体来说,对于邻域\mathcal{O}_{\mathbf{y}}的强度I(\mathbf{x}),K-means算法是一个迭代过程,它将聚类准则[19]最小化,可以写成连续形式

F_{\mathbf{y}}=\sum_{i=1}^{N} \int_{\mathcal{O}_{\mathbf{y}}}\left|I(\mathbf{x})-m_{i}\right|^{2} u_{i}(\mathbf{x}) d \mathbf{x}   (6)

mi是第i个聚类中心,ui属于区域Ωi的成员函数u_{i}(\mathbf{x})=1 for \mathbf{x} \in \Omega_{i} and u_{i}(\mathbf{x})=0 for \mathbf{x} \notin \Omega_{i},因为ui是该地区Ωi的隶属函数,我们可以重写F_{\mathbf{y}}

F_{\mathbf{y}}=\sum_{i=1}^{N} \int_{\Omega_{n}\cap \mathcal{O}_{\mathbf{y}}}\left|I(\mathbf{x})-m_{i}\right|^{2} d \mathbf{x}   (7)

针对式(7)中的聚类准则,利用mi≈b(y)ci近似聚类中心,定义了\mathcal{O}_{\mathbf{y}}中强度分类的聚类准则为

\mathcal{E}_{\mathbf{y}}=\sum_{i=1}^{N} \int_{\Omega_{n} \sigma_{y}} K(\mathbf{y}-\mathbf{x})\left|I(\mathbf{x})-b(\mathbf{y}) c_{i}\right|^{2} d \mathbf{x}(8)

其中K(y - x)作为非负窗口函数引入,也称为核函数,使得K(\mathbf{y}-\mathbf{x})=0 for \mathbf{x} \notin \mathcal{O}_{\mathbf{y}}。利用窗口函数,可以将聚类准则函数\mathcal{E}_{\mathbf{y}}重写为

\mathcal{E}_{\mathbf{y}}=\sum_{i=1}^{N} \int_{\Omega_{i}} K(\mathbf{y}-\mathbf{x})\left|I(\mathbf{x})-b(\mathbf{y}) c_{i}\right|^{2} d \mathbf{x}   (9)

该局部聚类准则函数是本文方法的基本组成部分。
局部聚类准则函数\mathcal{E}_{\mathbf{y}},评估邻域\mathcal{O}_{\mathbf{y}}在给定\left\{\mathcal{O}_{\mathbf{y}} \cap \Omega_{i}\right\}_{i=1}^{N} of \mathcal{O}_{\mathbf{y}}分类的强度。\mathcal{E}_{\mathbf{y}}值越小,分类效果越好。自然,我们定义整个区域\Omega的最优划分\left\{\Omega_{i}\right\}_{i-1}^{N},这样的局部聚类标准函数\mathcal{E}_{\mathbf{y}}\Omega的域内所有y的最小化。这可以通过对y在图像域Ω实现\mathcal{E}_{\mathbf{y}}的最小化的积分。因此,我们定义了能量\mathcal{E} \triangleq \int \mathcal{E}_{y} d \mathbf{y}

\mathcal{E} \triangleq \int\left(\sum_{i=1}^{N} \int_{\Omega_{i}} K(\mathbf{y}-\mathbf{x})\left|I(\mathbf{x})-b(\mathbf{y}) c_{i}\right|^{2} d \mathbf{x}\right) d \mathbf{y}   (10)

在这里,我们忽略域Ω的下标积分符号(如上面的第一个积分)如果集成在整个域Ω。图像分割和偏置场估计可以通过执行最小化能量关于regions \Omega_{1}, \cdots, \Omega_{N}, constants c_{1}, \cdots, c_{N}, and偏置场b
核函数K的选择是灵活的,例如,它可以是一个截断统一函数,定义为K(\mathbf{u})=a for |\mathbf{u}| \leq \rho and K(\mathbf{u})=0 for |\mathbf{u}|>\rho,a是一个是正的常数,这样∫K (u) = 1,本文选取核函数K为截断高斯函数,由

K(\mathbf{u})=\left\{\begin{array}{ll}{\frac{1}{a} e^{-|\mathbf{u}|^{2} / 2 \sigma^{2}},} & {\text { for }|\mathbf{u}| \leq \rho} \\ {0,} & {\text { otherwise }}\end{array}\right.    (11)

这里a是一个归一化常数,这样∫K (u) = 1,σ是标准偏差(或尺度参数)的高斯函数,ρ是\mathcal{O}_y邻域的半径。
我们把常数c_{1}, \cdots, c_{N}用向量\mathbf{c}=\left(c_{1}, \cdots, c_{N}\right)来表示,所以原能量函数可写为

\mathcal{E}(\phi, \mathbf{c}, b)=\int \sum_{i=1}^{N} e_{i}(\mathbf{x}) M_{i}(\phi(\mathbf{x})) d \mathbf{x}

其中

e_{i}(\mathrm{x})=\int K(\mathbf{y}-\mathbf{x})\left|I(\mathbf{x})-b(\mathbf{y}) c_{i}\right|^{2} d \mathbf{y}

所以

\mathcal{F}(\phi, \mathbf{c}, b)=\mathcal{E}(\phi, \mathbf{c}, b)+\nu \mathcal{L}(\phi)+\mu \mathcal{R}_{p}(\phi)

其中\mathcal{L}(\phi)\mathcal{R}_{p}(\phi)如下:

\mathcal{L}(\phi)=\int|\nabla H(\phi)| d \mathbf{x}
\mathcal{R}_{p}(\phi)=\int p(|\nabla \phi|) d \mathbf{x}

对于二分类N=2,M_{1}(\phi)=H(\phi)M_{2}(\phi)=1-H(\phi)

H_{\epsilon}(x)=\frac{1}{2}\left[1+\frac{2}{\pi} \arctan \left(\frac{x}{\epsilon}\right)\right]

最终求得

\begin{aligned} \frac{\partial \phi}{\partial t}=-\delta(\phi)\left(e_{1}-e_{2}\right)+\nu \delta(\phi) \operatorname{div}\left(\frac{\nabla \phi}{|\nabla \phi|}\right) & \\+\mu \operatorname{div}\left(d_{p}(|\nabla \phi|) \nabla \phi\right) \end{aligned}
d_{p}(s) \triangleq \frac{p^{\prime}(s)}{s}
\hat{c}_{i}=\frac{\int(b * K) I u_{i} d \mathbf{y}}{\int\left(b^{2} * K\right) u_{i} d \mathbf{y}}, \quad i=1, \cdots, N
u_{i}(\mathbf{y})=M_{i}(\phi(\mathbf{y}))
\hat{b}=\frac{\left(I J^{(1)}\right) * K}{J^{(2)} * K}
J^{(1)}= \sum_{i=1}^{N} c_{i} u_{i}\    J^{(2)}=\sum_{i=1}^{N} c_{i}^{2} u_{i}

SECTION IV.

Level Set Formulation and Energy Minimization

我们提出在(10)能量中的\varepsilon表达的\Omega_{1}, \cdots, \Omega_{N}。很难解决能量最小化问题来自这个表达式的大肠在本节中,能量E转换为水平集配方Ω1代表不相交的区域,⋯,ΩN水平集函数的数量,这些水平集函数的正则化项。在水平集公式中,能量最小化问题可以用已有的变分方法求解