融合双边滤波算子的多相水平集图像分割算法

史 娜， 孔慧华， 秦 鹏

(中北大学理学院数学系， 太原 030051)

超声检查是一种无辐射的安全检查，对人体没有创伤、没有损害，现已被广泛用于医学临床诊断。乳腺超声图像的准确快速分割对于妇科疾病的临床诊断具有重要意义，有助于医生提高诊断的正确率。乳腺超声图像具有低对比度的特点，其斑点噪声较多，肿瘤区域的边缘模糊且形状不规则。对于恶性乳腺肿瘤而言，待分割的目标边缘形状复杂、内部灰度较不均匀、目标与背景的对比度较低，因此，准确得到乳腺肿瘤病灶区域的边缘轮廓，是一项有难度且富有挑战的工作。

目前，针对医学超声图像分割问题的活动轮廓模型主要有以下3种类型:参数活动轮廓模型、几何活动轮廓模型、无边缘活动轮廓模型等。

根据曲线演化的不同，1988年，Kass等[1]提出了演化曲线参数化表示的活动轮廓模型，也称为蛇(snake)模型，该模型利用初始动态闭合曲线在感兴趣区域特征的约束下进行演化，当能量泛函达到极小值时，目标区域的边缘轮廓就是参数化表示的曲线位置。虽然蛇模型可以提供光滑闭合的边缘曲线，但无法自动处理运动界面的拓扑结构变化和局部形变。

1988年，Osher等[2]提出经典的几何活动轮廓模型，该算法将N维曲面或曲线的演化问题，转化为N+1维水平集函数的隐函数来求解，避免曲线的重新参数化表示，使曲线的演化更加灵活，具有自适应演化曲线拓扑结构变化的优势。2001年，Chan等[3]以水平集方法为基础，利用简化的Mumford-Shah泛函，建立分段常值的水平集分割模型(CV模型)；1997年，Caselles等[4]提出测地线活动轮廓模型；2002年，Vese等[5]在考虑全局区域信息的CV模型基础上，提出分片光滑和分段常值的多相水平集分割模型。

由于CV模型不能很好地解决灰度非均匀图像的分割问题，2007年，Li等[6]提出一种基于图像局部灰度信息拟合的活动轮廓(local binary fitting, LBF)模型，该模型以高斯函数为核函数对图像进行平滑处理；2009年，Wang等[7]结合全局能量泛函，将LBF模型改进为全局和局部能量信息相结合的图像分割模型。

近期，中外学者在超声图像分割领域，还提出基于多相水平集的活动轮廓模型[8-13]，王洁等[14]提出基于Vese-Chan模型的连续最大流方法，刘磊等[15]提出一种基于加权指数平均比率算子改进的CV模型，兰红等[16]利用改进的CV模型对磁共振图像进行分割。

针对乳腺超声图像中肿瘤区域的结构特征，提出一种融合双边滤波算子的多相水平集图像分割模型。改变以往将高斯函数作为核函数，引入双边滤波算子来进行边缘平滑，作为非线性滤波算子，它能有效保护边缘信息。结合超声图像的全局和局部灰度能量信息，以及多相水平集函数，使分割算法更适用于确定乳腺肿瘤区域的边界。

1 传统的无边缘活动轮廓模型

1.1 基于全局信息的Chan-Vase模型

传统的无边缘活动轮廓CV模型[3]，其算法以简化Mumford-Shah泛函为基础，利用初始演化曲线内外部图像灰度的分段常值，去逼近目标和背景区域。将曲线长度惩罚项和区域面积的惩罚项加入能量泛函，并用φ表示水平集函数，C表示演化曲线，对于待分割图像Ω(x，y)，所定义的能量泛函为

(1)

式(1)中:μ、ν、λ1、λ2分别为取值为正的参数，演化曲线C内部区域为Ω1，其平均灰度为c1；外部区域为Ω2，其平均灰度为c2，表达式为

(2)

(3)

式中：H(φ)通常为Heaviside函数；ε为很小的正数，即

(4)

δ(φ)=H′(φ)

(5)

通过变分法最小化能量泛函，得到水平集演化的偏微分方程，采用有限差分法进行数值近似计算，梯度下降流为

(6)

式(6)中：面积参数ν≥0；尺度参数μ>0。

1.2 基于局部信息的LBF模型

2007年，针对CV模型无法正确分割灰度非均匀图像的问题， Li等[6]提出一种局部二值能量拟合的活动轮廓(local binary fitting model，LBF)模型，该模型提供了有效的解决方案。

对于定义域为Ω的灰度图像I(x，y)，给定点x，y∈Ω，待分割图像的局部灰度信息的能量泛函为

|I(y)-f1(x)|2H[φ(y)]dy+λ2×

{1-H[φ(y)]}dy

(7)

式(7)中:H(·)为Heaviside函数;λ1和λ2分别为取值为正的参数;f1和f2这两个平滑函数分别为轮廓曲线C内部和外部的图像局部结构信息，而Gσ为标准差为σ高斯核函数，表达式为

(8)

f2(x)]dx+νL(C)

(9)

根据变分水平集方法，求解能量泛函式(9)的极小值，可得f1和f2的表达式为

(10)

最终，LBF模型的曲线演化方程为

(11)

式(11)中:κ为演化曲线的曲率;e1和e2满足:

(12)

2 融合双边滤波算子的多相水平集分割模型

2.1 双边滤波算子

双边滤波算子(bilateral filtering)是一种依赖像素点的空间邻近相似度和灰度差异度的非线性滤波器[17]。双边滤波在高斯滤波的基础上，引入像素的灰度值进行局部加权平均，针对图像的斑点噪声进行平滑时，能够较好地保持边缘结构特征。

双边滤波算子的滤波核权重ωi在空间近邻函数Gs基础上，增加灰度相似函数Gr，表达式为

(13)

(14)

(15)

(16)

(17)

式中:i为像素索引；j为窗口ωi覆盖区域的像素索引；x为像素空间位置坐标；I为像素灰度值；参数σs为空间近邻因子；参数σr为灰度相似因子；σs越大平滑效果越好，σr越小保边效果越好。

2.2 多相水平集方法

2002年，针对分段常值和分片光滑图像，Chan等[5]采用N个水平集函数划分2N个目标区域，建立了多相水平集的图像分割模型。

首先，引入两个水平集函数φ1和φ2，待分割区域被划分为四个互不相交、互不重叠的部分(图1):

图1 2个水平集函数划分4个图像分割区域

Heaviside函数和Dirac函数的表达式如下:

针对2个水平集函数划分4个图像区域的特征函数分别为M1、M2、M3、M4，表达式如下:

2.3 基于全局或局部的活动轮廓模型

针对乳腺超声图像边界模糊、感兴趣区域对比度较低、斑点噪声较大的情况，在CV模型和局部二值拟合模型的基础上，综合考虑乳腺超声图像的全局和局部区域的灰度能量信息，建立全局和局部二值拟合模型(local and global binary fitting model， LGBF)，其能量泛函的水平集函数表达式为

ELGBF(φ，f1，f2，c1，c2)=(1-α)EL(φ，f1，f2)+

αEG(φ，c1，c2)+ER

(18)

H[φ(y)]}dy]dx

(19)

H[φ(x)]}dx

(20)

式中:输入图像为I:Ω→R；点x，y∈Ω；C为闭合演化曲线；λ1、λ2≥0；图像区域Ω被曲线C分为目标(C内部)和背景(C外部)两个区域;c1为曲线C内部平均灰度值;c2为曲线C外部的平均灰度值;φ为水平集函数;H为Heaviside函数。局部能量泛函为EL(φ，f1，f2);全局能量泛函为EG(φ，c1，c2);α为图像的全局信息参数;GBF为双边滤波算子。

为保持乳腺超声图像演化曲线的光滑性，同时，避免重新初始化水平集函数，在总的能量泛函ELGBF中引入惩罚项ER，包括:长度惩罚项L(φ)和能量惩罚项P(φ)。当总能量泛函ELGBF达到最小值时，长度约束项L(φ)使得演化曲线C保持尽可能短且光滑，能量约束项P(φ)使得水平集函数φ在曲线演化过程中保持为符号距离函数。惩罚项ER的表达式为

(21)

式(21)中:长度惩罚项参数ν>0;能量惩罚项参数μ>0。

为更准确地获取乳腺肿瘤病灶区域的细节信息，可将ELGBF模型推广到多相水平集情况，满足乳腺肿瘤超声图像的多目标分割需求。具体的能量泛函定义为

ELGBF(φ1，φ2，f1，f2，f3，f4，c1，c2，c3，c4)=

ν[Lε(φ1)+Lε(φ2)]+μ[P(φ1)+P(φ2)]

(22)

(23)

式(23)中：f1、f2、f3和f4分别为演化曲线C1和C2所分割的四个图像区域的局部灰度二值拟合能量函数;c1、c2、c3和c4分别为演化曲线C1和C2所分割的4个图像区域中的全局平均灰度值。

由变分法，固定水平集函数φ1和φ2，最小化总的能量泛函ELGBF，得到满足Euler-Lagrange方程的局部灰度二值拟合的能量函数f1、f2、f3和f4，以及全局灰度平均值c1、c2、c3和c4，具体表达式为

(24)

(25)

对应最小化φ1的梯度下降流方程为

(26)

式(26)中:div(·)为散度算子；δε(·)为正则化的Dirac函数，拟合能量为F13=Hε(φ2)(e3-e1);F24=[1-Hε(φ2)](e4-e2)。

同理，可得φ2的梯度下降流方程为

(27)

在乳腺肿瘤超声图像的分割过程中，演化曲线在全局和局部区域能量泛函的相互作用下，内外力相互补充，当演化曲线接近病灶区域的边界时，反映乳腺超声图像局部结构信息的能量函数占主导地位；当演化曲线远离病灶区域时，表征乳腺超声图像全局结构信息的能量函数起主要作用。

3 实验结果与分析

3.1 实验数据集

使用的真实病例的乳腺超声图像进行实验，包括67幅恶性乳腺肿瘤的超声图像，图像的长宽均在300～400像素。该数据集中的所有乳腺肿瘤超声图像，均配有医生手工标记的病灶区域。在实际图像分割实验时，对每幅图像边缘的机器信息和图像参数信息进行剪裁。

3.2 实验结果的评价函数

乳腺肿瘤超声图像具有以下结构特征:①病灶区域的灰度值比周围正常组织小很多，并且灰度值相对均匀，但是病灶边缘不规则且拓扑结构复杂；②由于超声成像机制所限，可检测到的乳腺肿瘤区域的直径通常大于5 mm，若小于该阈值，可视为背景区域；③在对病人进行超声检查时，可将病灶区域限制在超声图像中心区域的窗口内。

为了定量化评价各模型的分割效果，采用的评价指标分别是:Dice相似系数、Jaccard相似度、假阳性率RFP及假阴性率RFN。具体计算公式为

(28)

式(28)中:N(·)为区域中像素的个数；Sg为医生手工标注的肿瘤区域；Sm为分割算法所获得的分割区域；O为Sg和Sm的公共区域。Dice相似系数和Jaccard相似度都能够表示算法分割结果Sm与医生手动标记Sg结果的重叠率，Dice相似系数和Jaccard相似度的值越接近于1，分割越准确。假阳性率RFP的值表示医生手工标记的病灶区域去除公共区域的占比，而假阴性率RFN的值表示分割算法所得分割区域去除公共区域的占比，RFP和RFN的值越接近于0，分割效果越好。

3.3 结果与分析

为验证本文算法的可行性，首先对仿真图像进行分割实验。将本文中提出的LGBF模型与CV模型和LBF模型的分割结果进行对比。然后，对乳腺超声图像数据集进行分割实验，利用评价指标对各算法的分割效果做出定量评价。最后，将融合双边滤波算子的LGBF模型推广到多相水平集情况，能够得到乳腺肿瘤内部的病灶细节信息。

实验环境为Intel Core i7-9700CPU@ 3.00 GHz，内存为8 GB，64 bit操作系统，实验软件为MATLABR2014b。实验中常用的参数设置为ν=0.001×2552，μ=1，λ1=1，λ2=1，ε=1，Δt=0.1。

在仿真图像的分割实验中，参数设置如下:全局结构信息的权重参数α=0.4，双边滤波算子的窗宽ω=4，空间近邻因子σs=4，灰度相似因子σr=25。将传统的CV模型、LBF模型与本文LGBF模型进行对比分割实验[分辨率如图2(a)左下方所示]。活动轮廓模型的初始演化曲线如图2(b)所示。由图2(c)可见，CV模型只考虑全局灰度信息，故无法正确分割灰度不均匀的图像。针对第1、2幅图像，LBF模型能够有效分割灰度不均匀图像，然而，在第3幅图像的分割实验中，由于只注重局部灰度信息，在成功分割染色体的同时，出现对染色体过度分割的现象(如图2中第3行第4列所示)。

本研究提出的LGBF模型充分考虑了图像的局部和全局灰度能量信息，采用双边滤波算子对图像进行预处理，能够有效分割出目标区域，分割效果良好，如图2(e)所示。

图2 灰度仿真图像的分割结果

然后，针对乳腺恶性肿瘤的超声图像数据集进行分割实验。在数据集中，病灶区域的边缘轮廓可以分为两类，一类是肿瘤区域的边界模糊,如图3中第1、2幅图像所示;另一类是肿瘤区域的边界较光滑，如图3中第3、4幅图像所示。

本文的融合双边滤波算子的LGBF模型，在进行乳腺超声图像进行分割实验时，根据乳腺超声图像的结构特征，实验参数设置如下:全局结构信息的权重参数α=0.3，双边滤波算子的窗宽ω=9，空间近邻因子σs=4，灰度相似因子σr=15。

图像分辨率和初始演化曲线如图3(a)和图3(b)所示。从图3可以看出，CV模型在对乳腺病灶区域进行分割时，演化曲线受全局灰度能量的作用，将图像灰度信息较一致的区域进行分割，导致病灶区域和正常组织区域混淆在一起，不能得到正确的病灶区域，如图3(c)所示。而传统的LBF模型过分注重超声图像的局部能量信息，在分割出病灶区域的同时，目标区域和背景区域出现了过度分割的现象，如图3(d)所示。

本文融合双边滤波算子的LGBF模型，在曲线演化过程中，利用双边滤波算子对超声图像进行平滑处理后，演化曲线受乳腺肿瘤超声图像全局和局部灰度能量结构信息的综合作用，能够准确逼近医生手工标记的病灶区域边界,如图3(e)所示。最后，将本文所提算法推广到多相水平集情况，得到乳腺超声图像病灶区域的细节分割效果,如图3(f)所示。利用评价函数，对不同算法的实验结果进行定量分析，分别对比Jaccard相似度、Dice系数、假阳性率RFP以及假阴性率RFN，发现本文算法的指标分别为0.938 7、0.945 1、0.074 2、0.083 3，均优于传统的CV模型和LBF模型，如表1所示。

表1 三种算法的分割结果对比

图3 乳腺肿瘤超声图像的分割结果

4 结论

融合双边滤波算子对乳腺肿瘤超声图像进行平滑处理，同时，根据乳腺超声图像中肿瘤区域与其他人体组织结构的灰度信息有相似性特点，所建立的能量泛函要综合考虑超声图像的全局和局部灰度能量信息。因此，本文多相LGBF模型优于传统的CV模型和LBF模型，能够准确得到病灶区域的边缘轮廓，分割曲线的平滑性优于其他算法，分割结果与医生手动标记的病灶区域相吻合。