自适应尺度参数的非均匀图像快速水平集分割
2021-04-15龚进慧张贵仓
龚进慧 张贵仓
(西北师范大学数学与统计学院 甘肃 兰州 730070)
0 引 言
图像分割在计算机视觉和图像处理领域仍然占有重要的地位。图像分割面临的主要挑战是如何分割具有噪声、强度不均匀和纹理复杂的图像。最近几十年来，学者们提出了许多图像分割方法。主动轮廓模型( Active Contour Model，ACM)[1]自提出以来一直受到人们的广泛关注。ACM算法能够提供平滑、封闭的轮廓线，以亚像素的精度覆盖目标边界。水平集方法是由Osher等[2]提出的活动轮廓模型族的成员之一。水平集方法用较高维函数的零水平集来表示轮廓，通常称为水平集函数，并且将轮廓的运动转换为水平集函数的演化。根据能量泛函中嵌入的图像性质，现有的水平集方法可分为两类：基于边缘的模型[3-5]和基于区域的模型[6-10]。
基于边缘的模型通常使用局部图像的梯度信息来迫使活动轮廓朝目标的期望边界移动，这对于分割具有清晰边界的图像特别有效。然而，这些模型对噪声和初始轮廓都很敏感，在梯度值比较小的弱边界上很容易发生边缘泄漏。基于区域的模型利用图像统计信息来推动活动轮廓覆盖目标边界。它们对噪声和轮廓初始化不太敏感，因为它们利用区域图像信息来驱动水平集演化。此外，它们还可以分割具有弱边界甚至无边界的图像。最流行的基于区域的模型之一是Chan-Vese(CV)模型[6]，这是众所周知的Mumford-Shah(MS)分割模型[11]的分段常量情况。文献[9，12-13]提出了一种基于区域的压力函数，并引入了一种快速的水平集演化公式，以使模型更快地收敛。随后，将局部灰度信息广泛应用于主动轮廓模型中，提高了噪声和强度不均匀性对图像分割的精度[14-18]。文献[14-15]将局部二元拟合(LBF)能量引入基于区域的活动轮廓模型的能量函数，以获得具有强度不均匀性的图像的满意分割结果。然而，基于局部区域的模型对初始轮廓也很敏感，容易陷入局部极小值。
强度不均匀性是现实世界图像中的常见现象，尤其是在磁共振(MR)图像中，其通常由成像装置的缺陷或照明变化引起。具有强度不均匀性的图像通常被描述为一个真实的分段常数图像乘以一个空间变化的光滑场[18]。这种空间变化的光滑场称为偏置场。偏置校正方法已被广泛研究，以处理具有强度不均匀的图像[19]。许多非均匀图像分割方法[20-23]都是基于偏置校正方法的。Li等[20]提出了一种新颖的水平集模型(LIC)，该模型考虑局部图像强度信息以演化水平集函数，并利用局部加权K均值聚类方法来评估偏置场。Zhang等[22]提出了一种局部统计活动轮廓模型(LSACM)，它引入了聚类方差信息。然而，这些模型对初始轮廓的位置也很敏感。此外，Li等[20]指出核函数的尺度参数应根据强度不均匀性的程度适当选择。然而，在LIC模型中却没有提到如何估计强度不均匀性的程度。Min等[24]提出了一种新的双极小化模型(DM)，该模型采用多层结构自适应地确定每个像素的尺度参数。
本文提出一种新的水平集非均匀图像分割方法。从非均匀图像模型出发，推导出最优分割平面，为每个像素提供最优分割。在平面上，提出一种新的基于区域的压力函数，并利用该函数在整个图像区域的水平集公式中定义了一个能量泛函。该方法通过对能量泛函进行最小化，在对非均匀图像进行分割的同时，对偏置场进行估计。此外，还引入了一个新的偏置场初始化，以提高初始轮廓的鲁棒性，这可以确保分割过程首先聚焦在有很大变化的区域，然后延伸到整个图像区域，而不管初始轮廓的位置如何。为了准确估计偏置场，针对核函数设计了一种新的自适应尺度参数。此外，在数值实现中采用了一种简单快速的水平集进化公式来提高效率。
1 背景模型
1.1 CV模型
Chan等[6]提出的CV模型假定物体和背景的强度属于两个可分离的区域。该模型的主要概念是找出图像的具体划分，将图像划分为两个聚类。对于域Ω中给定的图像I，通过最小化以下能量泛函来建立CV模型：
E=γ·length(C)+v·Area(inside(C))+
(1)
式中：γ≥0、v≥0、λ1>0、λ2>2是固定参数;length(C)是闭合轮廓线C的长度；Area是C的内部区域面积；c1和c2是轮廓C内部和外部的平均强度。在水平集方法中，等高线C、inside(C)和outside(C)的不相交区域表示为：
通过将式(1)相对于水平集函数φ最小化，相应的梯度流方程为：
v-λ1(I-c1)2+λ2(I-c2)2]
(2)
式中：δ是狄利克雷函数；▽是梯度算子，div()是散度算子。
同样，c1和c2是通过以下方式获得的：
(3)
(4)
式中：H(φ(x))是Heaviside函数。
CV模型对灰度均匀图像分割效果良好。此外，它对初始轮廓的位置不敏感。但是，当轮廓C内外的图像强度不均匀时，常数c1和c2可能不能准确地描述图像。因此，CV模型可能无法对灰度不均匀的图像进行分割。
1.2 张氏模型
张氏模型[9]将基于区域的模型与基于梯度的图像分割模型相结合，在该模型中引入基于区域的压力函数S。基于区域的模型的能量泛函被定义为：
(5)
式中：α是固定参数；μ是对象的灰度级下限。
下限μ可以作为图像的最优分割平面。因此，关于水平集函数φ最小化式(5)，并促使轮廓附着到灰度级大于μ的区域。相应的梯度流动方程为：
(6)
在CV模型中，当忽略前两项并且λ1=λ2=1时，能量泛函的梯度流动方程简化为：
(7)
1.3 LIC模型
为了正确分割强度不均匀的图像，Li等[20]提出了LIC模型。该方法基于一个常用的图像模型，其定义为：
I=bJ+n
(8)
式中：I是所测量的图像强度；b是偏置场；J是真实图像；n是加性噪声。在文献[20]中，关于真实图像J和偏置场b的两个假设如下：
(1) 偏置场b在图像域Ω中变化缓慢。
(2) 真实图像J可以近似地以在每个不相交区域中的常数来表示，即J(x)≈ci,x∈Ωi，其中{Ωi}i=1，2，…，N是图像域Ω的分区。
能量方程被定义为：
(9)
式中：K(·)是核函数。
当考虑LIC模型中的两相情况并引入水平集函数φ时，能量函数可以写成：
(10)
式中：M1(φ(x))=H(φ(x))，M2(φ(x))=1-H(φ(x))是每个集群的成员函数。
最佳的ci和b可通过以下方法获得：
(11)
(12)
式中：*代表卷积运算。
对于固定的c和b，通过求解以下梯度流方程，可以得到E相对于φ的极小值。
(13)
LIC模型采用加权K-均值聚类方法对偏置场进行评价。该方法通过偏置校正，可以分割出一些强度不均匀的图像。然而，LIC模型对初始轮廓位置的鲁棒性不强，核函数的尺度参数是经验性选择的，这是LIC模型的固有缺陷。另外，在LIC模型中采用有限差分格式进行数值实现。有限差分格式的时间步长T应满足CFL条件。为了保持数值稳定性，CFL条件将时间步长Δt限制为很小的长度，从而导致收敛缓慢。
2 本文方法
2.1 基于区域的压力函数设计
在LIC模型中，从式(13)导出以下方程：
(14)
卷积运算可以看作是一个获得加权平均值的过程。在假设(1)的基础上，可以得到：
|K*b2-(K*b)2|<ξ
式中：ξ是一个小的正常数。因此，可以假设K*b2≈(K*b)2，上面的方程可以表示为：
(15)
式中：α=2(K*b)(c1-c2)。
(16)
式中：K是具有标准差σb的高斯核函数。
在考虑N=2并引入水平集函数φ的情况下，可以将能量泛函写成：
(17)
类似于LIC模型，可以得到最优ci和b：
(18)
(19)
Heaviside函数H(φ)被正则化如下：
(20)
相应的梯度流动方程是：
(21)
2.2 新的偏置场初始化
当初始化偏置场b(x)=1,x∈Ω时，LIC模型在第一次迭代中退化为CV模型。在第一次迭代之后，零水平集可能远离目标的实际边界，这可能导致错误的偏置场估计和分割结果。因此，LIC模型对初始轮廓的位置敏感。为了克服这个缺点，引入了一个新的偏置场，它可以确保分割过程首先聚焦在物体真实边界附近的区域，然后扩散到整个图像区域。
由于整体图像强度变化的趋势与偏置场的趋势相似，因此通过保持这一趋势可以更好地估计偏置场。根据图像模型I=bJ+n和假设a，即在图像域Ω中，偏置场b变化缓慢，新的偏置场初始化定义为：
b0=K*(I/N0)
(22)
式中：N0是图像强度的平均值。
通过偏差校正，校正后的图像为：
(23)
S=I-K*(K*I)
(24)
由于目标边界位于图像强度变化较大的区域，新的偏置场初始化可以保证分割过程集中在真实边界附近的区域。不考虑初始轮廓的位置，因此可以先分割目标实际边界附近的图像。随着水平集的演化，远离目标真实边界的区域可以被分割。新的偏置场初始化可以提高初始轮廓定位的鲁棒性。
2.3 自适应尺度参数
核函数的尺度参数在偏置场估计中起着重要作用。小尺度参数可能使算法陷入局部极小值，而较大的参数可能无法正确估计偏置场。为了克服这一缺点，设计了一种新的自适应尺度参数，该参数依赖于局部图像的不均匀程度。根据图像模型I=bJ+n和假设a和b，可以看出，图像强度在图像边界附近的区域中变化很大，而在其他区域中变化缓慢。图像强度的变化是由图像在远离真实物体边界的平滑区域内的不均匀性所引起的。因此，可以使用局部区域中的方差来测量图像不均匀程度。可以根据平滑区域中局部图像不均匀性的程度来获得比例参数。
在桥台搭板进行设置期间，需要充分考虑到其坡度值，保证搭板长度的合理性。不仅要确保车辆行驶的安全性，更是要求桥体荷载力达到预期设计标准，能够承受更多的车辆负重，这样一来才可以将沉降控制在合理范围之内。虽然该项施工技术的应用有着较强的科学性，但还是不能够有效地解决桥头跳车问题，并且还可能导致路基压实度无法满足相关要求，严重的情况下还会导致搭板出现断裂的问题，引发严重的交通安全事故。所以，施工单位在施工之前需要充分的考虑沉降段的实际情况，确保搭板宽度和桥面宽度相适应，搭板边缘在缘石边缘的内部不能大于50cm，同时还要确保搭板与配筋设计方案的合理性，从而有效提高工程施工质量，给人们出行安全提供保障。
首先，采用中值滤波和双边滤波对图像I进行预处理，以减小噪声的影响，得到I1。然后，需要得到远离目标边界的光滑区域Ωh。利用拉普拉斯算子检测目标边界。由于大部分像素位于光滑区域，因此使用区域Ωh0粗略地表示远离目标边界的区域。
Ωh0=Ω(x)F(x) (25) 式中：F=K1*I1；T=τ·min(max(|F|),T0)是阈值；T0是恒定的阈值；K1是八领域的拉普拉斯算子；τ是时间步长。 局部区域中的图像强度变化越大，|F|的值越大。由于图像强度在远离物体真实边界的区域内变化缓慢，并且大部分像素位于光滑区域，所以F(x)在区域Ωh0附近，服从高斯分布。为了去除边界和噪声区域，在Ωh1中更精确地表示平滑区域。 Ωh1=Ω(x) max(cF-1.5σF,-T)< (26) 式中：cF、σF是F(x)在Ωh0区域的均值和标准差。为了使平滑区域更加精确，应该删除区域Ωh1的外部边缘，并获得平滑区域Ωh。 在获得平滑区域Ωh之后，得到以x为中心的(2ω+1)×(2ω+1)的正方形领域Nω(x)中的图像I1的局部方差D(x)。对于均匀图像，平滑区域的局部方差为零，尺度参数为无穷大。尺度参数与局部图像不均匀程度成反比。核函数的尺度参数可以设计为： (27) 式中：β是一个常数；cD是图像I1在区域Ωh中D(x)的平均值。从式(27)中可以清楚地看到，图像不均匀程度越大，尺度参数越小。可以根据图像自适应地选择适当的比例参数。对于σb>(length(I)+width(I))/4，(length(I)表示图像I的长度，width(I)表示图像I的宽度),图像强度是均匀的，可以用张氏模型对图像进行分割。 2.4.1算法流程图 本文算法流程如图1所示。 图1 本文算法流程 2.4.2算法过程描述 该算法的主要过程概括如下： 步骤1分别用式(22)和式(27)初始化偏置场b0和σb，并将水平集函数φ0初始化为以下二元函数： 式中:Ω0是图像阈Ω中的一个子集;∂Ω0是Ω0的分界线。 步骤2K=1∶m。 (1) 利用式(18)计算ci。 (2) 根据式(30)演化水平集函数φ。 (3) 用高斯滤波器对水平集函数进行正则化。 (4) 如果φ>0，则φ=1；否则，φ=-1。 步骤3利用式(19)计算b。 步骤4检查水平集函数是否满足停止条件，如果不，返回步骤2。 事实上，δ(φ)可以被|▽φ|取代[13]，那么式(21)的稳态等效于： (28) 基于观察，水平集函数在演化过程中必须保持为符号距离函数，为了防止它在界面附近太陡或者太平[26]，水平集函数必须满足条件|▽φ|=1。因此，式(28)可以简化为： (29) 为了保持数值实现的稳定性，在式(29)的迭代过程中，需要对水平集函数进行正则化。利用高斯滤波过程对水平集函数进行正则化[12]。使用加法运算符分裂(AOS)方案。水平集演化过程可以分为两部分，第一部分是求解式(29)，另一部分是对水平集函数进行正则化。采用带标准偏差σR的高斯滤波器对水平集函数进行正则化，并采用显式有限差分格式对式(29)进行离散，如下： (30) 式中：Δt表示时间步长；n表示迭代次数。 由于式(30)右边对于φ是独立的，可以选择方程的相对大的时间步长。因此，所提出的模型比LIC模型更快。实际上，每次迭代时对偏置场的估计是耗时的主要因素。当偏置场固定时，仅在一次迭代后，水平集不能达到稳定状态。因此，利用n次迭代来演化所提出的用于固定偏置场的方法中的水平集函数。 在实验中，当水平集函数的变化保持小于给定的阈值η时，迭代将自动停止，即，如果式(31)被满足，迭代将停止。 (31) //初始化水平集函数 initialLSF = c0*ones(size(Img)); initialLSF(30:90,50:90) = -c0; u=initialLSF; hold on; contour(u,[0 0],′r′); title(′Initial contour′); //初始轮廓 imagesc(Img,[0, 255]); colormap(gray); axis off; axis equal epsilon=1; b=ones(size(Img)); //初始化偏置场 K=fspecial(′gaussian′,round(2*sigma)*2+1,sigma); //高斯核 KI=conv2(Img,K,′same′); //图像卷积 KONE=conv2(ones(size(Img)),K,′same′); //高斯核与常值1函数的卷积，在计算LBF时会用到 [row,col]=size(Img); N=row*col; for n=1:iter_outer [u, b, C]= lse_bfe(u,Img, b, K,KONE, nu,timestep,mu,epsilon, iter_inner); //水平集迭代 if mod(n,2)==0 pause(0.001); imagesc(Img,[0, 255]); colormap(gray); axis off; axis equal; hold on; contour(u,[0 0],′r′);。 iterNum=[num2str(n), ′iterations′]; title(iterNum); hold off; end end 本文方法的性能在合成图像和真实图像上得到了广泛的评价。所有实验都是在MATLAB R2016a编程环境下进行的，使用的计算机型号为Lenovo G40，其配置为CPU 2.10 GHz、内存8 GB和Windows 10操作系统。实验中使用的图像像素各不相同，图2中第一幅像素为84×84,第二幅为98×100，第三幅为98×98，第四幅是MATLAB自带的一幅图像。图3中图像来自兰州第一人民医院眼科影像资源库，其像素为490×480。图4是在加入噪声的情况下的分割结果。图5(d)是合成图像，像素为103×93。经过多次实验，最终选择参数为：β=5.2，ω=1，T0=5，σR=0.9，Δt=1.0，η=1，m=5。其中时间步长可以取得比传统水平集方法大一些，本文实验过程中就分别取到了从0.1到100的宽幅范围。使用大的时间步长虽然可以加快迭代速度，但是选得太大可能会在边缘的地方引起错误。即大的时间步长和精确的边缘定位是矛盾的，通常情况下选择时间步长都小于10。 (a) 初始轮廓的原始图像 (b) 最终分割结果图2 本文方法在四幅强度不均匀图像上的实验结果 (a) 原始图像(b) 最终分割结果(c) 估计偏置场 (d) 偏差校正图像初始轮廓 图3 本文方法在两幅强度不均匀的真实图像中的应用 (a) 泊松噪声(b) 0.01的椒盐噪声(c) 0.03的椒盐噪声图4 加噪声图像分割结果 图5 三幅真实图像和一幅强度不均匀的合成图像的分割 在这一部分中，将本文模型与CV模型、LBF模型、LIC模型在强度不均匀的合成图像和真实图像上进行了比较。 图2为本文方法对四幅强度不均匀图像的分割。在这四幅图像中，强度不均匀性明显。结果表明，本文方法能够对灰度不均匀的图像进行分割，并取得了满意的分割效果。 与LIC模型相似，本文模型也能够分割灰度不均匀的图像，并同时估计偏置场。该方法应用于两幅强度不均匀的真实医学图像，如图3所示。可以看出，该模型得到了准确的分割结果，并且在偏置校正图像中，轮廓内部和外部区域的图像强度变得相当均匀。实验结果表明，该模型在图像分割和偏置场估计方面具有良好的性能。 图4是在加入噪声的情况下本文算法的图像分割结果。第1列的图像中加入了泊松噪声，第2列的图像中加入了密度为0.01的椒盐噪声，第3列的图像中加入了密度为0.03的椒盐噪声，第2行分别是各种噪声下的图像分割结果。实验结果表明，在加入上述各种噪声的情况下，本文算法仍可以较好地分割目标图像。 为了进一步验证本文算法的图像分割性能，并与其他算法进行比较，本文采用Jaccard Similarity(JS)等图像分割性能指标[25-28]，利用下面两个图像分割性能评价指标进行评价： 1) 面积重叠误差。 (32) 式中：A和M分别表示水平集算法分割区域以及手动分割区域的像素数量。 当面积重叠误差指标为零时，表示最精确分割，为100%时表示算法分割的结果与真实结果完全不重合，即错误分割。 2) 边界平均距离。设C是算法分割出的目标边界，对于边界C上的每一点Pi(i=1,2,…,N)计算Pi与真实目标边界S的距离dist(Pi,S)，然后求平均值。 (33) 如图5和表1分别给出了本文算法与其他经典算法的实验对比结果。在图5中给出了利用CV模型、LBF模型、LIC模型，以及本文方法对三幅真实图像和一幅强度不均匀合成图像的分割结果, 对于最后一幅图像，σR设置为0.5。CV模型采用全局强度信息的均值，如图5(b)所示，无法分割出强度不均匀的图像。LBF模型利用局部强度的均值来拟合被测图像，从而可以对灰度不均匀的图像进行分割。然而，该模型在能量泛函中引入了许多局部极小值，在某些情况下导致轮廓卡在不正确的位置上，如图5(c)所示。图5(d)和图5(e)分别是LIC模型与本文算法分割结果，因为它们在水平集演化过程中分割图像并同时估计了偏置场，因此比其他两种模型对于灰度不均匀图像的分割效果更好。但LIC模型对于灰度严重不均匀的合成图像的分割效果不是很理想，而本文模型不管是对真实图像还是对灰度严重不均匀的合成图像都取得了很好的分割效果。 表1 各算法的迭代次数、运算时间、面积重叠误差、边界平均距离对比 可以看出，本文算法在面积重叠误差、边界平均距离、迭代次数、运算时间上均优于其他几种算法。 为了进一步分析新偏置场初始化的影响，给出LIC模型与本文模型对于不同的初始轮廓分割结果的比较，如图6所示。结果表明，当初始轮廓位于物体内部时，LIC模型得到了正确的分割结果。然而，当初始轮廓同时包含目标和背景时，LIC模型无法对目标进行分割。而本文模型不管轮廓在什么位置，都取得了很好的分割结果，因此本文模型对初始轮廓具有较强的鲁棒性。 (a) 原始图像及初始轮廓 (b) LIC模型 分割结果 (c) 本文模型 分割结果 图6 LIC模型与本文模型对于初始轮廓位置的比较 本文提出一种新的水平集分割方法。在非均匀图像模型的基础上，推导出图像域的最优分割平面。在平面上，提出一种新的基于区域的压力函数，并在整个图像区域的水平集公式中定义了一个能量泛函。通过对能量泛函进行最小化来实现图像分割和偏置场估计。此外，引入了新的偏置场初始化以提高对初始轮廓的鲁棒性，并且为核函数设计了新的自适应尺度参数来准确地估计偏置场。在人工合成图像与真实图像上的仿真实验结果表明，该方法能够迅速、准确地对非均匀图像进行分割，并且对初始轮廓具有较强的鲁棒性。该方法分割结果能够辅助人工判读，用于辅助临床实习教学，对诊断提供一定参考信息。近几年许多学者基于水平集的图像分割方法做了大量的工作，该方法已经成功用于医学图像中病灶区域的分割，比如对于心脏的分割，大脑灰白质的分割，未来它还将会被应用于军事、航天、人脸识别等更多领域。
F(x)2.4 算法实现
2.5 部分实验编程
3 实验结果与分析
3.1 实验结果
3.2 实验分析
4 结 语
