复小波域混合概率图模型的超声医学图像分割

夏平 施宇 雷帮军 龚国强 胡蓉 师冬霞

医学影像已成为医生临床疾病辅助诊断最重要的手段之一[1−3];超声医学影像以其实时显示、对人体无伤害、适合软组织诊断等特点在眼部、心脏、腹部以及人体浅表部位等部位的疾病诊断中广泛应用.超声医学图像分割是医学影像分析的首要前提,是获取医学影像中病变区域与特定组织、器官信息的一种重要手段,其目的是从超声影像中分离出医生诊断感兴趣的区域,对辅助医生临床诊断、定位病变组织、规划手术治疗等方面显示出越来越重要的临床价值[4].

由于超声成像过程中的干涉现象及采集设备电子器件的随机扰动,使得超声图像存在大量不规则的斑点噪声,降低了图像中目标与背景区域的对比度,大大降低了超声图像的分辨率,使得图像中不同组织间的灰度差异性不明显,削弱了超声图像质量,从而影响了医生对病灶的识别,也增大了医生查出弥散性组织病变的难度.获取超声医学图像中特定组织器官或病变区域定量信息的过程中,超声医学图像分割是不可或缺的手段[5];超声影像的分割除传统分割手段[6]外,近年来随着统计学理论[7]、模糊集理论[8]、神经网络[9]、多尺度理论[10]的发展与应用出现了很多新的分割算法[11−14].文献[15]针对医学图像分割中区域主动轮廓模型易受目标和背景区域面积比的影响,且对初始位置敏感的问题,提出模糊C-均值聚类(Fuzzy C-means,FCM) 协作改进CV (Chan-Vese) 模型(FCM-CV) 的图像分割算法;该算法通过在CV 模型中增加能量权值函数消除面积比的影响后,用FCM 粗分割结果指导设定改进CV 模型零水平集的初始位置,从而提高了图像分割的稳定性和准确性.随着小波技术的发展与应用,基于小波域马尔科夫随机场(Markov random field,MRF) 分割方法的研究十分活跃;文献[16]提出了小波域中树结构化MRF 的图像分割算法,按照图像分类层次树的结构形式,该算法在小波域多分辨率中构建了序列MRF 模型,它们相互嵌套,并通过在小波域的相邻尺度间和同一分辨率内的两个层次上递归运算实现图像的有效分割;该算法在表示图像结构信息的同时,较好地描述了图像的非稳态性,得到较好的分割效果.但小波分析本身缺乏平移不变性,且只能在水平、垂直、以及45◦三个方向建模图像的特征信息,因而,基于小波的图像分割在获取目标边缘的连续性不理想,检测结果存在一定“空洞”现象.相对于小波分析,双树复小波变换[17]采用两路实离散小波分析实现复小波分析,改善了小波分析中对平移的敏感性,并实现了信息的多方向获取,在一定程度上克服了目标分割的“空洞”问题;文献[18]利用双树复小波矩形状信息对医学淋巴结进行精确定位,对已探测的淋巴结进行形状特征提取,判断其是否为真实的淋巴结;然后,应用水平集DRLSE(Distance regularized level set evolution) 模型对淋巴结进行轮廓描述,实现淋巴结完整的分割.

超声医学影像中,其成像过程所造成的斑点噪声以及相关组织纹理是影响图像分割准确性的重要因素,结合超声影像的特点,本文提出了复小波域混合概率图模型的超声医学图像分割方法.利用双树复小波变换(Dual tree-complex wavelet transform,DT-CWT) 的多分辨率分析具有的多方向性和近似的平移不变性,捕捉超声图像的局部弱特征信息[19];多分辨率的相邻层间构建的有向图模型确立父-子节点标号相互联系,构建MRF 模型来建立同层复小波系数标号与邻域间的联系;而噪声是统计独立的,从而,通过构建混合概率图模型建立了当前节点的标号与其祖先节点及同层邻域间节点之间联系,有效地获取超声影像中有用信息,实现超声影像的可靠分割.

1 超声影像的双树复小波域稀疏描述

DT-CWT 采用两路实小波对信号进行分解以实现复小波变换[17],其中,一路为CWT 的实部,另一路为其虚部;复小波函数:Ψ(t)=Ψh(t)+jΨg(t),其中,Ψh(t)、Ψg(t) 均为实小波.其分解过程用树A和树B 两个实小波分解树表述,对应树A、树B的滤波器分别为h0(k)、h1(k) 和g0(k)、g1(k);如图1(a).每层CWT 分解得到2 个低频分量和±15◦、±45◦、±75◦等6 个高频方向子图像.如图1(b) 所示.

树A、树B 对应的滤波器组除满足完全重构的条件外,还需使Ψ(t) 解析[17],即:

将树A、树B 两个实DWT 分别由两个方阵Hh和Hg表示,则DT-CWT:

超声图像I进行复小波分析,有:Ch=Hh·I和Cg=Hg·I;其中,Ch、Cg分别为DT-CWT 的实部和虚部系数,Ch+jCg为DT-CWT 稀疏描述的高频复值方向子带,从而得到6 个方向子带:

式中,每一子带的复小波实部与虚部由v={1,2}确定,mmm=(m1,m2) 为位置坐标,6 个方向子带用u ∈{1,2},i ∈{1,2,3} 标示.

图1 DT-CWT 变换及其子带方向Fig.1 DT-CWT and the sub-band direction

相对于实小波分析,DT-CWT 继承了小波分析的多尺度分析特性和稀疏表示能力,同时具有更好的方向选择性与近似的平移不变性.

2 复小波域超声图像的混合概率图模型

N × N格网S上的超声医学图像I,经J -1 级复小波分析后形成J幅不同尺度的高频子带图像,J幅图像大小集合为其中,S(0)为原始图像;S(j)对应第j级复小波分析的子带图像.本文基于格网集合上构建复小波域超声医学图像的混合概率图模型.如图2 所示.

2.1 标记场的混合概率图模型建模

定义复小波系数的标记场:

式中,标记FS(n)对应于分辨率n的复小波系数 CS(n) 的标记场,标记场F(n)=表示尺度n上空间(i,j) 的标记值.类似地,FFFS(J−1)对应最低分辨率复小波系数CCCS(J−1)的标记场,FFFS(0)对应于原始图像的标记场.

标记场描述了像素点特征值在“父-子”间的属性传递性以及同分辨率中局部的相关属性,不同分辨率尺度上的观测特征展现对应的区域特征,因而标记场的多分辨率描述体现了这两种属性共同作用的结果.

依据复小波域系数的持续性,相邻尺度间复小波系数的标记具有Markov 性,即多分辨率分析构成的层间Markov 链表征了层间的相关性,并且,将父节点标记属性传递给其子节点,尺度间的“父-子”节点联系形成复小波域中系数间的贝叶斯网络结构.如图3(a) 和(b) 所示.

图2 DT-CWT 域系数向量结构Fig.2 Coefficient vector structure of DT-CWT domain

标记场建模中,尺度间标记具有一阶Markov转移特性[20].多分辨率分析中,相邻尺度间的标记场服从一阶Markov 分布,用转移概率描述父、子复小波系数标号的依赖关系为:

式(5) 的一阶Markov 性反映了高分辨率尺度上的标记场仅受相邻的低分辨率尺度标记场影响,与其他尺度的标号无关.如图3(c) 所示.

图3 DT-CWT 域标记场贝叶斯网络模型Fig.3 Marking field Bayesian network model in DT-CWT domain

尺度间标记的传递属性仅考虑相邻层间“父-子”节点的相互作用,如图3(b) 中父节点“a”与子节点“b”等;节点标记传递属性数学建模:

式中,fi+1、fi分别表示父、子节点标记值.设医学图像分为K类,其传递属性用父节点标记值到其子节点标记值转移概率p来描述,则标记转移概率矩阵:

式中,ε ∈[0,1]为相邻尺度间“父-子”节点具有相同标记的概率.

式(8) 转移概率定义使式(7) 构成一强对称矩阵,此时,“父-子”节点标记在传递过程中信息容量达到最大.本文采用ε=1,即高分辨率尺度标记场的初始值由相邻低分辨率尺度标记场的实现来确定.

标记场建模中,尺度内标记具有Markov 性.即尺度n中复小波系数的标记场是以∆ij为邻域的MRF,即:

式中,Pf(n)>0,∆ij为尺度n内的(i,j) 位置的邻域系统.

同一尺度标记场是以邻域系统的结构作为条件,MRF 建模中常采用1 阶或2 阶邻域系统,本文采用2 阶邻域系统描述同一尺度标记场的结构间联系.如图4 所示.文献[20]已证明,邻域系统基团势函数参数αm、βi、γi、δ等与标记场取值无关,因此,本文在复小波域MRF 建模中,对指定阶数的所有邻域系统均采用相同的势函数参数.

图4 DT-CWT 域标记场的2 阶邻域系统Fig.4 The second order neighborhood system in DT-CW

据Hammersley-Clifford 定理[21],MRF 和Gibbs 分布具有等效性,使用MLL (Multilevel logistic model) 模型来表示标记场的先验分布[7,12].本文采用各向同性的双点MLL 模型描述像素空间邻域的作用关系,仅考虑其二元势函数来计算标记之间的局部约束条件.

本文中,定义随机场∆={∆ij|(i,j)∈S=N×N}是S的邻域系统,∆ij是(i,j) 的邻域点集,用各向同性的双点MLL 模型来计算标记之间的局部约束条件,其势函数Uc与标记场的局部条件概率P分别用式(10) 和式(11) 描述.

式中,Z为概率分布归一化因子,即分配函数,c是二维图像基团,β是对应的势团参数,δ(·,·) 为信息采样函数,具有如下特性:是位置(i,j) 的邻域中不等于的相邻位置的数目.

尺度间与尺度内两种标记场体现了多分辨率分析中不同层间、同一层内的复小波系数特征相互作用与影响后表现出的特征属性,尺度间作用构成了一个有向图结构,在此结构中,父节点标记属性是以式(7) 转移概率矩阵形式传递给子节点,因此,多层传递中形成了一个贝叶斯网络结构:尺度内联系构成了无向MRF 的图模型,每个节点标记与邻域∆ij节点标记建立联系,从而,在尺度内构成了无向图模型:二者共同构建了复小波域空间标记场的混合概率图模型.

2.2 观测特征场建模

超声图像I做J-1 层复小波分解后高频系数为CCC,每一尺度复小波用对应层号n(1≤n ≤J -1)表示.将每一分辨率上对应位置不同频带方向的复小波系数构成向量,形成了对应的向量图像,则此向量图像表示了该分辨率中某一方向图像的观测特征,层号n=0 对应原始图像.则:

分辨率尺度n、位置(i,j) 的复小波系数的特征向量:

经J -1 层复小波分解后,定义图像I的格网S结构:

由式(15) 可知,J -1 分解后格网中共包含6J+1 个子带,尺度上的格网:

各尺度上对应的复小波系数C的向量场描述的特征场的位置在(i,j)∈S(n)所表示的范围内.不同位置集上相应的复小波系数形成观测特征场表示为:

式中,C(J−1)为最粗尺度的复小波系数,只有一个节点,C(0)对应原始图像.

特征场建模中,考虑到任意节点上复小波系数值的概率由每一观测特征由所属标号和所在分辨率的观测特征值确定,不同标号下观测特征相互统计独立.即:

图5 给出了前列腺囊肿图像复小波分析后第1、第2 层各方向高频图像(图5(a)) 及系数统计图(图5(b)),统计分析表明,信号经DT-CWT 后子带系数的边缘概率分布呈现尖峰、长拖尾的非高斯分布,如图5 所示,反映了DT-CWT 分析在尺度内对信号的“聚集”能力和稀疏表示能力.由此可知,子带高频系数值大部分位于“0”值附近,少数系数值较大,为此,用值“大”或“小”两状态描述子带系数的分布规律;基于此,文献[22]采用两个不同方差的高斯混合模型来近似逼近子带系数总体分布.本文中尺度n的观测特征值采用高斯混合模型建模,该分辨率中每一标号的观测特征值采用高斯模型拟合,则:

图5 医学图像的DT-CWT 分解及其高频系统统计Fig.5 Medical image of DT-CWT and Statistical characteristics of the coefficient

式中,分别表示尺度n上标号的均值向量和协方差矩阵;Q为特征向量维数.

2.3 误分割率最小的超声医学图像分割

分割的目标是对复小波域中每个复小波系数观测特征与其分配的标记一一对应,记:(Cs,Fs);用条件概率P(Cs|Fs) 表示复小波系数对标记的依赖关系,用边缘分布P(Fs) 表示标号的先验概率分布,在给定复小波系数Cs下标记Fs的概率:

分割过程中标记的获取与其观测特征应相适应,其误分割率Pe表示为:

最佳分割的目标是分配一组最佳标记Fs使图像的误分割率达到最小,即在已知观测复小波系数Cs的基础上,求使得P(Fs|Cs) 最大的分割

式(22) 等价于最小化能量问题.能量函数为标记场能量U和特征场能量U之和.即:

由式(20)～(24),采用基于局部条件概率的确定性算法,通过逐点更新图像标记完成图像分割.

2.4 分割算法实现

超声图像I大小为N ×N,对其进行J -1 层复小波分解,采用高斯混合模型建模各尺度特征场,每一标号的观测特征场采用高斯模型拟合,用MLL模型建模各个尺度的标记场,使用EM 算法完成模型中参数估计.分割算法步骤:

1) 对超声图像进行复小波多分辨率分解,获得观测场,获取最粗尺度上的初始分割.本文采用Kmean 聚类算法获取尺度n=J -1 上的初始分割结果;

2) 令尺度n ←n-1;

3)E步.用MPL 方法估计式(25)、式(26) 的模型参数[8]:

式中,n表示多分辨率分解的尺度,l表示迭代次数,m表示标记类别;

4)M步.用式(25) 和式(26) 估计的参数,采用ICM 求解式(23) 的最小能量,获取该尺度上标记场结果;

5) 尺度内重复步骤(3) 和步骤(4),直到满足迭代条件(∆U=Ui-Ui−1≤10−3) 即停止迭代,得尺度n的标记场;

6) 尺度间迭代:将尺度n的标记场结果映射到最邻近的精细尺度n-1 上,作为尺度n-1 的初始分割,重复步骤2),直至获取最细尺度的标记场;

7) 进行形态学处理,实现超声医学图像分割.

3 实验仿真与性能分析

3.1 实验参数与定量评价指标

采用囊肿、乳腺癌、输尿管末端囊肿、左肾肿物、前列腺等超声医学图像进行实验,医学图像源由市中心人民医院影像科提供.对比实验时采用复小波域K-means 算法、空域MRF 算法、小波域MRF 算法以及本文算法对图像进行分割.

参数选取:复小波分析采用Q-shift DT-CWT,进行3 层DT-CWT 分解;K-means 算法中,分为目标与背景2 类;MRF 中,势团参数选β=0.8;迭代运算条件:前后两次迭代能量之差∆U ≤10−3.

定量评价指标采用概率Rand 指数(Probabilistic rand index,PRI)[23]、信息变化指数(Variation of information,VoI)[24]、全局一致性误差指数(Global consistency error,GCE)[25]、边界偏移误差指数(Boundary displacement error,BDE)[26]等客观评价指标,从不同的角度反映分割结果的优劣以及对算法的比较与评价.

1) 概率Rand 指数(PRI) 用于度量待评测结果与参考结果之间的属性共生的一致性.待评测聚类标记场FFF,图像中的每个像素点i标记号fi,用信息采样函数δ(·,·) 指示两个像素点i、j是否属同一类.

参考分割结果记为ℜ,用集合ℜ={r11,r12,···,rMN} 描述,则待评测分割结果与参考结果之间PRI定义[23]:

PRI ∈[0,1],PRI值越大两个聚类结果的一致性越高,分割结果越好.

2)V oI从信息论的角度度量不同聚类之间的距离,反映分割信息量的丢失.定义[24]:

式中,I(F,ℜ) 是F和ℜ所共有的互信息.V oI ∈[0,∞),V oI值越小两个聚类结果的一致性越高,分割效果越好.

3)GCE用于度量不同分割结果之间的一致性误差[25].GCE ∈[0,1],GCE值越小表示细化误差越小,分割效果越好.

4)BDE待评测结果和参考结果边缘像素间平均距离来度量分割结果[26].BDE ∈[0,∞),值越小表示两者边界差异越小,分割效果越好.

3.2 实验结果与性能分析

结果如图6、图7 所示.由图6、图7,从清晰度角度看,基于空域MRF、小波域MRF 以及本文算法分割的图像均较清晰;从分割区域效果看,本文算法很好地将目标区域分割出来,其次,表现较好的是小波域MRF 算法,而其他算法表现差强人意;其主要原因在于,复小波域K-means 分割中并没有很好地利用复小波分析优良的特征提取特性;而空域MRF 算法中,在图像分割区域内部构建MRF,充分考虑了区域内部的一致性问题,但忽视了不同目标区域边界的各向异性,从而,导致应用该算法进行图像分割时易产生目标边缘模糊;小波域MRF 算法相比于本文算法,主要缺乏平移不变性以及方向特征信息表示不足.4 种分割算法的定量比较见图6,通过PRI、V oI、GCE、及BDE等指标的评价,本文算法的PRI指标高于其他算法0.014 以上;V oI指标低0.09 以上;GCE指标低0.03 以上;BDE指标低1.10 以上.本文算法的分割结果从边缘、轮廓的清晰度、细节信息捕捉等方面均好于其他三种算法.

图6 超声医学图像4 种算法分割结果比较Fig.6 Comparison of segmentation results of four algorithms for ultrasonic medical images

图7 医学超声图像4 种算法定量评价指标比较Fig.7 Comparison of quantitative evaluation indexes of four algorithms for medical ultrasound images

由图6 和图7,从清晰度角度看,基于空域MRF、小波域MRF 以及本文算法分割的图像均较清晰;从分割区域效果看,本文算法很好地将目标区域分割出来,其次,表现较好的是小波域MRF算法,而其他算法表现差强人意;其主要原因在于,复小波域K-means 分割中并没有很好地利用复小波分析优良的特征提取特性;而空域MRF 算法中,在图像分割区域内部构建MRF,充分考虑了区域内部的一致性问题,但忽视了不同目标区域边界的各向异性,从而,导致应用该算法进行图像分割时易产生目标边缘模糊;小波域MRF 算法相比于本文算法,主要缺乏平移不变性以及方向特征信息表示不足.4 种分割算法的定量比较见图6,通过PRI、V oI、GCE、及BDE等指标的评价,本文算法的PRI指标高于其他算法0.014 以上;V oI指标低0.09 以上;GCE指标低0.03 以上;BDE指标低1.10 以上.本文算法的分割结果从边缘、轮廓的清晰度、细节信息捕捉等方面均好于其他三种算法.

4 结论

本文提出了复小波域多分辨率MRF 模型的超声医学图像分割算法.复小波分析中每一分辨率提供了6 个方向的高频特征信息,并具有近似的平移不变性,适应了图像的非平稳性,从而能较好地捕捉超声图像弱特征信息;误分割率最小算法充分考虑复小波域中特征信息的层间、层内信息的相关性,利用了复小波域中每一尺度内同标号区域的特征信息分布规律,弥补了空域MRF 分割中参数估计复杂、小波域MRF 分割中缺乏平移不变性和特征提取方向性差的不足,得到了更好的分割效果.

理论推导和实验仿真验证了本文方法应用于超声医学图像目标检测的有效性以及实现超声医学图像目标检测所具有良好的检测性能.

未来及目前正着手做的工作包括:1) 本文算法在实现误分割率最小的能量函数最优解的求解中采用ICM 算法,正着手引入图割、置信度传播等算法实现能量最优解,探讨通过优化求解能量最小化模型,提高算法的分割效率;2) 探讨其他算法,如深度学习方法,应用于细胞、视网膜等医学图像分割;3)针对肺部CT 影像,在三维测量基础上,结合本文算法研究检测肺结节的方法预测和分析肺部肿瘤恶性情况,为临床诊断与治疗提供辅助手段.