基于边力函数能量泛函的图像分割算法

唐银敏 李自纳 胡 静

1(郑州大学信息工程学院 河南 郑州 450000) 2(南阳农业职业学院信息工程学院 河南 南阳 473000)

0 引 言

图像分割在医学、定位、计算机视觉等领域发挥了极其重要的作用[1]。然而,由于低分辨率、高噪声和模糊的边界,它具有极大的挑战性。随着网络时代的不断发展,实时性以及准确性要求使得实现快速高精度的图像分割成为研究发展的热点[2]。

活动轮廓模型(Active Contour Model,ACM)的基本思想是利用轮廓曲线提取目标,轮廓曲线包括两类:显式ACM和隐式ACM[3-4]。显式ACM在求解过程中使用参数方程显式表示演化曲线[5]。隐式采用有符号距离函数(Signed Distance Functions,LSF)代替参数方程来表示演化曲线[6]。实际上,噪声或弱边界往往存在于真实图像中。由于噪声的影响和灰度值的不连续性,标准的基于边缘的最小二乘法无法检测出图像的边界。因此,Li等[7]提出了一种距离正则化水平集演化方法,以保持LSF的正则性,并将等高线曲线推向所需的位置。基于区域的最小二乘法利用图像的区域统计信息来定位目标边缘。比较基于边缘的最小二乘法,这些模型对噪声不敏感,例如区域可缩放拟合模型[8]、局部图像拟合模型[9]和鲁棒噪声区域ACM[10]。特别的,Chan等[11]提出了一种基于区域的Chan-Vese模型,它能成功地分割出亮度均匀的图像。局部图像拟合模型是另一种基于区域的LSM模型,它以截断高斯函数作为局部区域约束来分割图像,可以分割出强度不均匀的图像,但对强噪声图像敏感,RSF模型也能分割出非均匀图像,但当图像是严重噪声或非均匀时,分割效果较差,而且该模型会使集函数过于平坦,不利于水平集函数的正则性。

为解决上述方法中存在的问题,提出一种基于有序统计滤波能量驱动的鲁棒主动轮廓模型。通过实验结果分析可知提出的方法能够更加快速、准确地分割图像,并且对初始轮廓和参数具有较强的鲁棒性。

1 基础理论

1.1 几何主动轮廓模型

几何主动轮廓(Geodesic Active Contour,GAC)模型是一种基于边缘指示函数的主动轮廓模型[12]。该模型提出了如下能量函数:

(1)

式中:Eint代表内部能量项,即长度项;Eext代表外部能量项,即面积项。

Eint(C)=α|C′(q)|2

(2)

Eext(C)=λg|▽I(C′(q))|2

(3)

式中:α和λ是常数;▽表示梯度算子;g是边缘指示函数。g定义如下:

(4)

式中:Gσ是具有标准差σ的高斯滤波函数;*表示卷积。可以得出梯度下降流方程:

(5)

1.2 距离正则化水平集演化模型

距离正则化水平集演化模型(Distance Regularized Level Set Evolution,DRLSE)的能量函数如下:

(6)

式中:μ、λ和α是常数,分别是距离正则项、长度项和面积项的系数。式(6)中的g表示边缘指示函数:

(7)

式中:Gσ是具有标准差σ的高斯滤波函数。式(6)中的p表示双阱势函数,可以写成:

(8)

Hε和δε函数表达式为:

(9)

(10)

可以通过变化的微积分获得以下梯度下降流方程:

αgδε(φ)

(11)

式中:dp(x)=p′(x)/x表示演化速度函数;p′(x)表示双阱势函数的导数函数。

1.3 CV模型

设I:Ω→R为输入图像,C为闭合曲线,提出了如下能量函数:

(12)

式中:λ1、λ2和v为常数;outside(C)代表轮廓C以外的区域;inside(C)代表轮廓C内的区域;c1和c2是两个常数,分别近似于区域outside(C)和inside(C)中的强度平均值。式(12)中的前两项是数据项,它们用于驱动曲线停止在目标边界处。第三项是长度项,用于平滑和缩短曲线。从式(12)可知,只有当轮廓C位于目标边界时,能量ECV才能达到最小值。

在水平集方法中,轮廓C⊂Ω由Lipschitz函数φ:Ω⊂R的零水平集表示,称为水平集函数,水平集函数如式(13)所示。当x点分别位于曲线C的外部和内部时,则水平集函数φ分别为正和负;否则,水平集函数φ等于零。

(13)

能量函数ECV重写如下:

(14)

式中:Hε(x)和δε(x)分别是近似Heaviside函数和近似Dirac函数[13]。

(15)

(16)

使用标准的梯度下降方法来最小化式(14)中的能量函数,可以得到以下梯度下降流方程:

(17)

其中c1和c2是:

(18)

最后,用迭代法φk+1=φk+Δt×(∂φ/∂t)求出水平集函数φ。c1和c2涉及图像强度的全局属性,它们不包含任何局部信息,因此,在对强度不均匀的图像进行分割时,分割结果将是错误的。

2 顺序统计滤波模型

2.1 二阶微分特性

GAC模型和DRLSE模型通常使用面积项αgδε(φ)来加速零水平集的演化速度。但由面积项驱动的演化过程是单向的,这意味着水平集函数不能自适应地选择演化方向,只能根据设定方向向内收缩或向外扩张,而不能根据图像信息调整演化方向。因此,本文的目的之一是寻找一种既能准确定位边界位置又能指导曲线演化方向的边缘信息。

以CV模型为例,CV模式下的数据拟合项为-[λ1|I(y)-c1(x)|2-λ2|I(y)-c2(x)|2],c1和c2分别代表轮廓C外部和内部区域中强度的平均值。并且设置λ1=λ2=1,则CV模型中的数据拟合项可以重写为FCV(φ,c1,c2)=-(λ1|I-c1|2-λ2|I-c2|2)=(I-c2)2-(I-c1)2。I(x,y)表示图像中的一个点,因此(I-c1)2表示c1和点I(x,y)之间的差的平方,而(I-c2)2表示c2和点I(x,y)之间的差的平方。(I-c1)2和(I-c2)2都代表一次微分计算。之后,(I-c2)2-(I-c1)2表示二次微分计算。因此,数据拟合项FCV(φ,c1,c2)也可以称为二阶微分数据项。二阶微分的零点是灰度值剧烈变化的地方,即目标边界的位置。

2.2 边力函数

在本文中,利用顺序统计滤波来构造具有二阶微分特性的边力函数(Edge Force Function,EFF)。令Ω→Rd为图像域,当d=1时,它代表灰度图像,而d=3时,它代表彩色图像。首先,定义一个以x为中心、宽度为ω的正方形邻域Ωx。

其次,通过使用MATLAB中的二维顺序统计滤波器,将局部窗口Ωx中的像素值从低到高排序,然后获取图像强度的最大值和图像强度的最小值分别为f1和f2。

f1(x)=max[I(y)|y∈Ωx]f2(x)=min[I(y)|y∈Ωx]

(19)

式中:I(y)表示局部窗口Ωx中的图像强度值。对于宽度为ω的某点x,可以直接根据式(19)计算f1(x)和f2(x)的值。

然后,介绍了以下能量函数:

(20)

式中:C是闭合轮廓;inside(C)和outside(C)分别表示轮廓C内部和外部的区域;λ1和λ2是常数。当式(20)中的能量函数被最小化时,轮廓C将位于边界上。

曲线C可以由Lipschitz函数φ的零水平集代替,因此,式(20)中的能量函数可以重新定义如下:

(21)

δε(φ)在式(22)中是Hε(φ)的导数。Hε(x)和δε(x)表示如下:

(22)

通过用最速下降法最小化式(21)中的能量函数,可以得到梯度下降流方程[14]:

(23)

式中:α是可以控制OSF模型分割速度的正常数;参数s表示为控制OSF模型分割范围而增加的图像的标准差。ei(x)表示如下:

(24)

然后,将EFF定义如下:

fEFF(x)=α(λ1e1(x)-λ2e2(x))

(25)

演化水平集函数:

φn+1=φn+Δt·(∂φ/∂t)

(26)

式中:Δt是时间步长。(Ai+1-Ai)/Ai<10-5是曲线演化的停止标准,其中Ai和Ai+1分别表示迭代前后轮廓所包围的区域。

最后,使用优化的正则化函数φL和优化的长度函数φC分别对水平集函数进行正则化并平滑和缩短曲线。

(27)

φC=Gk*φL

(28)

式中:φL是优化的正则化函数,其目的是改善过零区域的斜率,并限制两端的斜率,以确保水平集函数具有很大的正则性;参数β用于控制过零区域的斜率,以确保分割效率;Gk是具有标准偏差k且大小为round(2×k)×2+1的高斯核。这两个函数的加入不仅减少了计算量,而且消除了长度和惩罚项的系数。

该模型的优势体现在以下几个方面:首先,由于f1(x)和f2(x)与大小受Ωx控制的局部区域中的图像强度有关,因此,OSF模型可以分割强度不均匀的图像。其次,由于EFF的二阶微分特性,曲线可以在演化过程中自适应地选择演化方向。第三,EFF是在迭代之前计算的,因此,OSF模型仅需要在迭代期间更新δε(φ),即可大大提高分割速度。第四,OSF模型大大提高了初始化和参数的鲁棒性。最后,由于添加了优化的正则化函数φL和优化的长度函数φC,极大地降低了计算复杂度,并且消除了长度和惩罚项的系数。

3 实验与结果分析

3.1 实验设置

OSF模型的步骤如下:

步骤1设置各种参数并初始化水平集函数。

步骤2根据式(19)计算f1(x)和f2(x)。

步骤3根据式(24)和式(25)计算EFF。

步骤4根据式(26)更新水平集函数。

步骤5分别根据式(27)和式(28)计算φL和φC。

步骤6如果达到迭代次数,则停止,并输出分割结果,否则返回步骤4。

在OSF模型中使用以下参数:Δt=0.1,ε=1,k=2.5,β=10,α=1.5,ω=7,λ1=λ2=1,C0=1。初始水平集函数φ0初始化为二进制步进函数,在零水平集外部和内部分别取1和-1。OSF模型是在2.6 GHz Inter Core i5个人计算机上的MATLAB R2017b中实现的。

3.2 初始轮廓的鲁棒性

在本节中,将选择五幅强度不均匀或边界较弱的不同图像,以证明OSF模型对初始轮廓具有鲁棒性。对于每幅图像,用三种不同的形状设置初始轮廓,分别是矩形、圆形和三角形。每个形状的初始轮廓在大小和位置上都不同。

初始轮廓和分割结果如图1所示。因为EFF是在迭代之前计算的,它与轮廓无关。因此,可以选择不同的初始轮廓,在不同的位置嵌入、包含或远离目标。

图1 不同形状、大小和位置的初始轮廓的分割结果

在图1中,白虚线代表初始轮廓,可以看到初始轮廓具有不同的位置、形状和大小。灰线代表最终的分割结果,可以看到,在初始轮廓具有不同形状、大小和位置的情况下,OSF模型可以精确地分割每幅图像。因此,可以得出一个结论,即OSF模型对初始轮廓具有鲁棒性,能清楚显示出不同形状、大小位置的轮廓,分割结果明确,OSF模型的鲁棒性效果明显。

3.3 分割实验

在本节中,将利用十幅强度不均匀、边界弱、背景不均匀或对比度低的典型图像来进行OSF模型的分割实验。初始轮廓、演化过程和最终分割结果如图2所示。在图像d和j中,设置β=9,在图像g中,设置α=1.0只是为了获得更好的分割效果。

图2 图像的分割结果

在图2中,规则的线代表初始轮廓,不规则的线代表演化曲线。在每个组中,第一列代表原始图像和初始轮廓。第二和第三列表示演化过程,最后一列表示准确且令人满意的最终分割结果。因此,可以得出结论,OSF模型可以正确地分割强度不均匀或边界较弱的图像。

3.4 对比实验

在本节中,将比较OSF模型的分割结果与一些经典的基于区域模型的分割结果,包括RSF模型、LIF模型、RSF&LoG模型和LPF模型。选择图3中的图像k-图像r来证明OSF模型的良好分割速度和分割精度。图3显示了具有相同初始轮廓的OSF模型、RSF模型、LIF模型、RSF&LoG模型和LPF模型之间的比较。根据花费的时间和每个模型的分割精度,将评估每个模型的性能。并根据戴斯相似系数(Dice Similarity Coefficient,DSC)[15]估算结果:

图3 对比实验结果

(29)

式中:ST表示通过手动获得的真实区域;SE表示通过实验获得的目标区域。DSC的值越接近1,则分割精度越高。分析结果列于表1。

表1 每个图像的迭代次数、花费的时间和DSC

在图3中,第一列代表原始图像和初始轮廓。第二到最后一列分别代表RSF模型、LIF模型、RSF&LoG模型、LPF模型和OSF模型的分割结果。表1显示了每幅图像的迭代次数、花费的时间和DSC。

在表1中,可以注意到OSF模型花费的时间比其他模型少得多,因为EFF是在曲线演化以及优化的正则项和长度项之前计算的。RSF模型和RSF&LoG模型花费的时间相对较长,因为它们每次迭代都具有四个卷积。LPF模型的分割时间也很长,这是因为要计算预拟合函数fs(x)和fl(x)。由于每次迭代的计算成本低,LIF模型的分割时间相对较短。此外,所提出模型的DSC在图像k、l、m、o、p、q、r中排名第一,在图像n中排名第二。因此,可以得出一个结论:与其他四个基于区域的经典模型相比,OSF模型可以更准确、更快地分割图像。

4 讨 论

4.1 边力函数

图4(a)显示了图像的EFF的计算结果,可以看出图中有明显的黑白环,黑环和白环相交的位置是目标边界所在的位置;图4(b)显示了横截面上具有二阶微分特性的边缘力能量EEFF的值。可以注意到,边缘力能量的值在目标边界外部和附近的区域为正,而边缘力能量的值在目标边界内部和附近的区域为负。下行中所示的二阶微分的零点是点A、B、C、D,它们表示目标边界所在的位置。因此,通过以上分析可以知道,可以通过计算EFF的二阶微分零点来检测目标边界。

图4 EFF的计算结果和边缘力能量的值

4.2 自适应演化曲线

EFF在本文中的另一个作用是吸引演化曲线以自适应地演化到目标边界。图5显示了一次迭代后的合成图像的分割结果。规则的线表示初始轮廓,不规则的线表示一次迭代后的轮廓。如图5所示,点a在目标内部,位于初始轮廓上。因此,根据图5可以知道点a的边缘力能量值为负,并且根据梯度下降流方程,点a应该移动到能量较高的区域。并且,由于能量值在初始轮廓之外的区域中为正,因此,点a向外扩展到点a1,移动速度为|∂φ/∂t|Δt,其中Δt是时间步长。

图5 曲线的演化方向

点b也在初始轮廓上,但是它位于目标边界之外。因此,根据图5,可以知道点b的边缘力能量值为正,因此根据梯度下降流方程,点b应该移动到能量较低的区域。并且,因为能量值在初始轮廓内的区域中为负,所以点b向内收缩到点b1,移动速度为|(∂φ/∂t)|Δt。

点c位于目标边界上,因此根据图5,可以知道点c的边缘力能量的值等于零,即|(∂φ/∂t)|=0,则运动停止。

4.3 参数β、k、α和ω的影响

4.3.1 关于参数β

式(27)对本文中的水平集函数进行了正则化。利用式(27)中的参数β控制过零区域的斜率,以保证水平集函数的正则性。首先选择图6中的图像s、t、u、v、w、x来显示参数β的有效性。图6中的第一行代表原始图像和初始轮廓;第二行表示没有参数β的分割结果;第三行表示有参数β的最终分割结果。可以明显地注意到,没有参数β的最终分割结果是错误的,而具有参数β的分割结果是正确的。

图6 在没有参数β和有参数β的情况下获得的分割结果

在分割实验中,参数β在图像d和j中进行了微调,以实现更好的分割。因此,在本节中,将选择图6中的图像v来讨论参数β的鲁棒性。图7(a)显示了原始图像和初始轮廓。图7(b)-图7(h)表示当参数β等于7到13时获得的最终分割结果。当参数β等于8到12时,最终的分割结果令人满意。因此,可知OSF模型对参数β具有较强鲁棒性。

图7 具有不同参数β值的分割结果

4.3.2 关于参数k

式(28)用于平滑和缩短演化过程中的曲线。但是,高斯滤波会模糊图像的边界,并会导致丢失少量信息。因此,式(28)中的参数k应根据图像信息设置。在本节中,将设置k=2.0、k=2.5、k=3.0,以证明在OSF模型中选择的k=2.5是最合适的。在图8中,选择图6中的图像u、v、w、x和图像y进行比较实验。第一行表示原始图像和初始轮廓。第二行到最后一行分别代表k=2.0、k=2.5、k=3.0的最终的分割结果。表2记录了每幅图像的迭代次数和花费的时间。

图8 使用参数k的不同值进行分割的结果

表2 不同参数k下图像u到y的迭代次数和花费的时间

从图8和表2中可以发现,当k=2.0时,图像v和y均未成功分割。同时,当k=3.0时,图像v和w分割失败。可以使用参数k的不同值成功分割图像u和x。从表2可以看出,当k=2.5时,迭代次数和花费的时间最少。当k=2.0和k=2.5时,可以成功分割图像w。从表2可以看出,当k=2.5时,迭代次数和花费的时间都更少。类似的,当k=2.5和k=3.0时,可以将图像y分割,当k=2.5时,迭代次数和花费的时间会更少。因此,可以得出结论,在OSF模型中选择k=2.5是最有利的。

4.3.3 关于参数α

式(23)中使用的参数α可以控制OSF模型的分割速度。在分割实验中,为了更好地分割,在图像g中对参数α进行了微调。因此,选择图8中的图像x来证明OSF模型对于参数α是鲁棒的。

如图9所示,第一列表示原始图像和初始轮廓,而第二列到最后一列显示最终的分割结果,其中α=0.5、α=1.5、α=2.5、α=3.5、α=4.5和α=5.5。可以看到,只有当α=0.5时,OSF模型才不能正确分割图像。相反,当α=1.5、α=2.5、α=3.5、α=4.5和α=5.5时,最终的分割结果都是正确的。因此,可以得出一个结论,即OSF模型对参数α具有鲁棒性。

图9 不同参数α值的分割结果

4.3.4 关于参数ω

局部窗口Ωx中的参数ω决定了所提出的模型一次要处理的像素数。因此,参数ω可以在一定程度上控制所提出模型的分割速度。以图像y为例,通过更改参数ω的值并保持其他参数固定来说明参数ω的有效性。如图10所示,当ω=4时,最终的分割结果没有捕获到完整的目标。相反,当ω=8时,最终分割结果包括非目标部分。但是,当ω=5、ω=6和ω=7时,最终结果都是正确的,这表明所提出的模型对参数ω具有鲁棒性。此外,当ω等于5、6和7时,表3中记录了所提出模型所需的迭代次数和分割时间。从表3中可以发现ω=5、ω=6和ω=7都可以得到理想的结果,但是ω的值越大,所提出的模型所需的迭代次数和分割时间就越少,这可以证明所提方法中使用的参数ω可以控制所提模型的分割速度。因此,通过以上分析,可以得出关于参数ω的以下结论:(1) 所提出的模型在一定范围内对参数ω具有鲁棒性。(2) 所提出的方法中使用的参数ω会影响所提出模型的分割速度,这是因为参数ω可以一次控制所提出的模型要处理的像素数量,因此ω的值越大,则一次处理更多像素。

图10 不同参数ω值的分割结果

表3 不同参数ω下图像y的迭代次数和花费的时间

5 结 语

针对传统图像分割算法中存在的图像分割速度慢以及初始轮廓和参数鲁棒性差等问题,提出一种基于边力函数能量泛函的图像分割算法。通过实验结果分析可得如下结论:

(1) 提出的OSF模型能够快速、准确地获取目标边界,实现精确的图像分割。

(2) 由于能量函数中的f1和f2仅与图像强度有关,所以OSF模型对初始轮廓具有较强的鲁棒性。

(3) 正则项的引入大大地减少计算量,并消除了长度系数和惩罚项系数,使得算法对参数选择也具有一定的鲁棒性。