一种基于二维模糊聚类的超声图像分割方法转让专利
申请号 : CN200610116339.3
文献号 : CN100595789C
文献日 : 2010-03-24
发明人 : 汪源源 , 余锦华
申请人 : 复旦大学
摘要 :
本发明属于超声图像处理技术领域，具体为一种二维模糊聚类结合斑点噪声滤波器和亮度补偿的B型超声图像分割方法。该方法将基于图像亮度信息的传统模糊聚类方法扩展为同时在像素和像素邻域值进行模糊更新的二维模糊聚类；在二维模糊聚类目标函数中引入由各向异性扩散斑点噪声滤波器提供的邻域信息，加强二维模糊聚类对斑点噪声的鲁棒性；将亮度不均匀假设为乘性噪声的基础上，通过在二维模糊聚类目标函数中引入亮度补偿因子，利用图像中目标和背景存在的均匀性差异，加强二维模糊聚类对亮度不均匀噪声的鲁棒性。本发明能够在分割图像中同时抑制斑点噪声和亮度不均匀的影响，提供一种超声图像分割的新的有效方法。
权利要求 :
1、一种基于模糊聚类的超声图像分割方法，其特征在于具体步骤为：首先，提出2DFCM的概念，将邻域信息引入到聚类过程中，以提高算法灵活性和收敛速度；接着，针对斑点噪声的抑制，采用各向异性扩散滤波器SRAD为2DFCM提供邻域信息，提高算法对斑点噪声的鲁棒性；最后，将由信号衰减引起的亮度不均匀视为乘性噪声，根据超声图像获取过程中信号处理的特点，修改2DFCM的目标函数，引入亮度补偿因子，构造匀质化的二维模糊C均值聚类，该聚类记为2DHFCM。
2、根据权利要求1所述的基于模糊聚类的超声图像分割方法，其特征在于具体步骤为：(1) 对原始图像进行SRAD滤波在滤波前需采用EM算法对图像的概率分布进行拟合，并在此基础上计算图像目标背 景分割的门限，得到分割图，然后通过加窗的方法找到图像中包含背景最多的区域，根据 公式(13)计算初始扩散门限，再通过SRAD得到邻域约束数据i^(《1^1,2,•••，"};(2) 采用高斯混合分布模型估计2DFCM的初始参数首先根据原始数据1={^ /=1,2,'","}和滤波数据^={《| /=1,2,'","}计算二维直方图， 然后使用方差为0. 625的二维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要 高斯分布数设为聚类数尺，直方图的各个极大值点作为聚类中心；(3) 2DFCM迭代，迭代的目标函数为公式(4)，约束条件为^〜=1，执行以下步骤直至聚类中心在两次迭代前后的改变量小于一个很小的常数£=0.001,① 用公式(6)更新隶属度矩阵，a，② 用公式(7)更新聚类中心m,， S"③ 计算更新前后聚类中心的改变量llM,-M。wl12，这里M。w表示更新前的M，风， 表示更新后的M;(4) 2DHFCM初始化，这里2DHFCM为匀质化的二维模糊C均值聚类，将2DFCM收敛后的隶属度矩阵以及聚类中心作为2DHFCM的初始值，并用公式(20) 计算初始的亮度补偿值；以滤波数据f={《|;=1，2,〜，"}计算一维直方图，然后使用方差为 1.0的一维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为期 望的分类数c;(5) 2DHFCM迭代，迭代的目标函数为公式(17)，约束条件为力^ =1 ，执行以下步骤直至迭代后存在的不同聚类中心数量小于期望的分类数c,① 用公式(18)更新隶属度矩阵^，② 用公式(19)更新聚类中心^,^,③ 用公式(20)更新两个亮度补偿向量，A、 x,④ 如果聚类中心更新后，两聚类中心间的欧氏距离小于W，则将两聚类中心合并， 统计不同聚类中心的数量；
说明书 :
一种基于二维模糊聚类的超声图像分割方法
技术领域
本发明属于超声图像处理技术领域，具体涉及一种二维模糊聚类结合斑点噪声滤波器 和亮度补偿的超声图像分割方法。 技术背景
图像分割是指按某个特征如亮度或纹理，将图像划分为几个区域的过程。在医学图像 分析领域，图像分割是其它应用，如组织测量、解剖结构分析和组织定征的前提条件和关 键步骤。尽管已提出很多医学图像分割的方法，但其中大多数方法是对信噪比相对较高的 核磁共振图像或CT图像提出的。超声诊断因其安全、方便、无痛、价廉等优点己成为临 床不可缺少的诊断方法之一，但是由于超声图像中大量噪声的存在，使超声图像分割仍然 是一个悬而未决的难题。
在超声成像中，当人体组织的结构尺寸与入射超声波波长相近或小于波长时，超声束 发生散射。相位不同的散射回波相互干涉产生斑点噪声。作为超声图像中最为显著的噪声， 斑点噪声的存在降低了超声图像的质量，掩盖或模糊了图像中对诊断非常重要的细节信 息。因此，近年来众多学者陆续提出了维纳滤波[1]、自适应加权中值滤波[2]、小波软阈值 去噪[3]、各向异性扩散[4]等抑制斑点噪声的方法。除斑点噪声外，超声图像中由信号衰落 造成的亮度不均匀也是图像信噪比降低的主要原因。亮度不均匀使同一组织内部存在亮度 的变化，同时也降低了不同组织间的对比度，为超声图像的自动分割制造了更大的障碍。 因此，也提出过许多用于解决医学图像中的亮度不均匀问题[5'6]的方法，但很多方法仍然是 针对核磁共振图像提出的，为实现超声图像的自动分割，必须同时抑制斑点噪声和亮度不 均匀对分割结果的影响。 发明内容
本发明的目的是提出一种能同时抑制斑点超声和亮度不均匀对分割影响的超声图像自 动分割方法。
本发明提出的超声图像分割方法是一种基于模糊逻辑并融合斑点噪声抑制和亮度补偿 的超声图像分割方法，其具体步骤如下：首先，提出二维模糊C均值(2DFCM)的概念，在 将邻域信息引入到聚类过程中的同时，提高算法灵活性和收敛速度；接着，针对斑点噪声 的抑制，采用各向异性扩散滤波器SRAD为2DFCM提供邻域信息，提高算法对斑点噪声
的鲁棒性；最后，在将由信号衰减引起的亮度不均匀视为乘性噪声的基础上，根据超声图像获取过程中信号处理的特点，修改2DFCM的目标函数，引入亮度补偿因子，构造匀质 化的二维模糊C均值聚类(2DHFCM)。
下面对各步骤作进一步具体描述。
相关概念：模糊C均值聚类(FCM: Fuzzy C-Means)171。
若将图像像素值组成"个样本的集合义Hx,l /=1,2/-,"}，那么图像的分割问题就转化为 将这"个样本分成c个聚类的问题。设M^—户1,...,}为每个聚类的中心，〜是第/个样 本对第乂类的隶属度函数，则用隶属度函数定义的聚类损失函数JV^，可写成：
^鼎=tl>,i,『 （1)
其中6>1，是一个可控制聚类结果模糊程度的常数， 一般情况下6取2。图像的优化分割， 就是通过迭代寻找聚类中心,与隶属度值〜，使目标函数Ao/取最小。在t〜-i的限制
条件下，根据拉格朗日条件极值法，可得到目标函数/K^取条件极小值的必要条件：
(Zll ||2》/(*-"
=,广m』J (2)
1>a n、
〜=^- (3)
传统FCM算法在进行图像分割时，仅利用了像素的灰度信息，因此聚类结果对噪声 较为敏感。在FCM的目标函数中加入邻域约束是提高FCM抗噪性能的有效手段^。所谓 邻域约束就是在计算当前像素对各类的隶属情况时，将与目标像素相邻的各像素的隶属情 况作为约束条件，其依据是自然图像相邻像素存在的高相关性。邻域约束的加入使聚类结 果保持了图像空间上的连续性，抑制了噪声对聚类结果的影响。构造邻域约束的方法可参 考文献[8]。
1、 2DFCM
假设原始图像数据表示为^={^|/=1,2,'","}，原始图像经过滤波器后得到的数据表示为 i={?,h'=i,2，一,W。图像的分割问题就是要将二维样本《Z, f》分成c个聚类，每个聚类由
一个二维向量表征^={( ,，^)l产l,…,c)。这里^为f的聚类中心，产l,…,c由于原始图 像数据与其滤波数据之间的强相关性，可以用g的分类结果来影响其对应像素x,的分类结
果，即当x,对以".为中心的第y'类具有较高的隶属度时，《应该对以^为中心的第j'类具
有相似的高隶属度。3?,的分类对x,分类的影响程度可以通过一个加权常数a控制。根据经验，a—般取6-8。根据以上分析，本发明提出的2dfcm的目标函数可写为：formula see original document page 6 在约束条件^^=1下，目标函数的求解可采用Lagrange(拉格朗日)乘子法，构造新的 无约束条件的目标函数：
formula see original document page 6
入为由Lagrange乘子法引入的中间变量。将目标函数分别对〜，〜和^求一阶导数并赋 值为零，可得到目标函数取极小值的迭代公式：
formula see original document page 6
二维模糊聚类结构的提出较传统fcm具有以下优势：邻域信息引入到2dfcm的目标 函数中，增强了聚类对噪声的鲁棒性；在二维方向上寻找聚类中心，有潜力提高算法的收 敛速度；结构灵活，可以针对不同类型噪声，采用不同的滤波器提供邻域约束信息。
为抑制斑点噪声对聚类结果的影响，本发明采用各向异性扩散滤波器sradw来提供 2dfcm所需的邻域信息。各向异性扩散的基本原理是在抑制边缘附近区域平滑的同时， 加强平行于边缘方向的平滑。在srad中，滤波是一种扩散的过程，扩散的力度由下式决 定：formula see original document page 6
其中《(x，;;;0称为瞬时变化系数，它充当了扩散过程中的边缘检测器；的(^称为噪声比例函
数，它充当了扩散门限。扩散系数c(《)的构造应该是在均匀区域加强平滑，而在垂直于边
缘的方向抑制平滑。因此，srad在实际应用中的效果还取决于如下两个因素： 一是噪声
图像中边缘检测的有效性。即：《(x,少;/)能否有效区分由噪声或者纹理造成的小瞬时变化系
数和由边缘造成的大瞬时变化系数；二是扩散门限值估计的准确性。目卩：的W能否准确判
别何时加强或者抑制扩散过程。文献[4]提出的方法很好地解决了第一个问题，并在文献[9]
中进行了扩展，直接将其应用到超声图像的边缘检测中。对于第二个问题，文献[4]中仍然依靠经验选取初始扩散门限。
为给2DFCM提供准确的邻域约束信息，本发明使用二状态的瑞利、高斯混合分布对 超声图像灰度分布进行拟合，并采用EM算法实现混合分布的分解；根据分解结果预测图 像中斑点噪声均匀分布的区域；通过对均匀区域统计特性的分析获取各向异性扩散的扩散 参数。二状态的瑞利、高斯混合分布如下：
二 c^p(x I + o:2/3(x I Gotw."") (9)
其中a!和g2为混合系数，满足cq+a产l，且有：
1
2;rcr
exp
2-2 2a2
(10)
(11)
EM混合分布分解的方法可参考文献[IO]。当EM算法收敛后，根据似然分割准则[1()]，得 到目标区域、背景区域的判决门限：
2[1ogC57cr2) + C^/a2]/1/〃2 -l/o":
(a,(12)
/2[log(cr，） + «2Zai]/|lAT2 -1//?2
得到判决门限后将图像分割为目标区域和背景区域(噪声区域)。然后将图像划分为30x30 的小块，选择噪声分布最均匀的块，即包含噪声点最多的块，根据下式计算初始扩散门限：
va柳)]
(13)
其中z(O表示斑点噪声均匀分布的区域，vw[z(0]和^分别表示此区域的方差和均值； 2、 2DHFCM
在2DFCM中本发明通过引入各向异性扩散作为邻域约束来提高聚类对斑点噪声的鲁 棒性，为增强算法对实际超声图像的适应性，还需要考虑解决由信号衰减所造成的亮度不 均匀问题。
亮度不均匀是在各类医学图像中存在的一个普遍问题，如核磁共振图像。迄今为止解 决亮度不均匀问题较成功的方法是将其假设为一种乘性噪声[5'6]，然后利用对数压縮将其转 变为加性噪声以简化后继的处理过程。假设J^(x^2,…,x4， 1^{力,》,...具}分别表示理想 背散射信号数据和实际背散射信号数据，这里理想数据指不包含信号衰减的数据，则理想 数据和实际数据间的关系可表示如下：
"=JT, xS,.， ls^" (14)
超声图像获取时包含了对数压縮、低通滤波、插值等操作，所以通常将观测超声图像灰度信号和实际背散射信号的关系简化为对数压縮关系。根据图像灰度信号与实际背散射信号 的对数压縮关系，则可将理想图像像素和观测图像像素的关系表示如下：
V/e{l,2,..."} (15)
其中5={艮艮…,y^即为对观测像素P(",》,…,W的亮度补偿。因为各向异性扩散滤波 器对观测超声图像具有均值保留特性，因此对各向异性扩散滤波后数据的亮度补偿可以采 用与原始图像相同的方法。设经过各向异性扩散后的数据？ = Jj具有亮度补偿
^={",72, ...Jn}，则F和i?的关系表示如下：formula see original document page 8
在2DFCM的目标函数中引入亮度补偿，即将式(4)中的x,用乂-A代替，3f,用义.-》代替， 就得到2DHFCM的目标函数：
formula see original document page 8
在约束条件|;〜=1下，目标函数(17)的求解可采用与2DFCM相同的方法，得到目标函数
>1
取极小值的迭代公式：
formula see original document page 8
在2DHFCM算法迭代前，需要预先决定的初始参数为聚类数、初始聚类中心和初始亮 度补偿估计。由于2DHFCM在每次参数迭代时需要更新三组相互影响的参数，如果初始 参数不准确则会由误差的传递造成算法的错误迭代，因此2DHFCM的正确收敛依赖于初 始参数的准确性。这里利用自底向上松弛标注法和高斯混合分布模型为2DHFCM提供初 始参数估计和保证算法的正确收敛。自底向上松弛标注法是在聚类中引入空间纹理信息的 方法，该方法通常包括两个处理过程：首先通过一种聚类算法如K均值法来预测初始的标 注概率，然后通过重新调整分类将空间连续性引入最终的分类结果中。这里先用2DFCM对图像进行预分类，然后将其分类的结果作为2DHFCM的初始状态；接着再用2DHFCM 重新调整分类的结果，达到同时消除斑点噪声和亮度不均匀的目的。在执行2DFCM前同 样需要给定聚类数及初始聚类中心，本发明中采用高斯混合分布模型来估计初始参数。即 将原始数据Z和i构成的二维直方图视为多个二维高斯分布混合而成，混合分布的个数即 为初始聚类数，混合分布的中心位置坐标即为初始聚类中心，实现过程可参考文献[ll]。 附图说明
图1、对包含斑点噪声超声图像分割结果的比较。其中图l(a)为理想图像，图l(b)为 加入斑点噪声的模拟图像，图l(c)为SRAD的滤波结果，图l(d)为FCM分割的结果，图 1(e)为文献[8]分割的结果，图l(f)为2DFCM的分割结果，图l(g)为2DHFCM的分割结果。
图2、对包含斑点噪声和亮度不均匀超声图像分割结果的比较。其中图2(a)为加入斑 点噪声和亮度不均匀的模拟图像，图2(b)为亮度不均匀模型，图2(c)为FCM分割的结果， 图2(d)为文献[8]分割的结果，图2(e)为2DFCM的分割结果，图2(f)为2DHFCM的分割结 果，图2(g)为对原始图像的亮度补偿图，图2(h)为对SRAD滤波图像的亮度补偿图。
图3、对一幅胎儿肾积水超声图像分割结果的比较。其中图3(a)为原始图像，图3(b) 为SRAD滤波的结果，图3(c)为2DFCM的分割结果，图3(d)为2DHFCM的分割结果，图 3(e)为对原始图像的亮度补偿图，图3(f)为对SRAD滤波图像的亮度补偿图。
图4、对图3的聚类过程说明。图4(a)为聚类数为30的结果，图4(b)为聚类数为10 的结果，图4(c)为聚类数为6的结果，图4(d)为聚类数为3的结果，图4(e)为聚类数为2 的结果。
图5、对一幅胎儿主动脉弓超声图像分割结果的比较。其中图5(a)为原始图像，图5(b) 为SRAD滤波的结果，图5(c)为2DFCM的分割结果，图5(d)为2DHFCM的分割结果，图 5(e)为对原始图像的亮度补偿图，图5(f)为对SRAD滤波图像的亮度补偿图。
图6、对图5的聚类过程说明。图6(a)为聚类数为30的结果，图6(b)为聚类数为10 的结果，图6(c)为聚类数为6的结果，图6(d)为聚类数为4的结果，图6(e)为聚类数为2 的结果。
图7、对一幅乳腺肿瘤超声图像分割结果的比较。其中图7(a)为原始图像，图7(b)为 2DHFCM的分割结果，图7(c)为对SRAD滤波图像的亮度补偿图，图7(d)为对分割图像提 取边缘轮廓的结果。 具体实施方式
以下通过整个算法的实现步骤，进一步对发明中的各个组成部分加以描述 1对原始图像进行SRAD滤波在滤波前需采用EM算法对图像的概率分布进行拟合，并在此基础上计算图像目标背 景分割的门限，得到分割图，然后通过加窗的方法找到图像中包含背景最多的区域，根据 公式(13)计算初始扩散门限。再通过SRAD得到邻域约束数据f ={《U=l,2,*",/7}。
2采用高斯混合分布模型估计2DFCM的初始参数
首先以根据原始数据;=1,2,'","}和滤波数据》={3f,|片1,2,…,W计算二维直方 图，然后使用方差为0.625的二维高斯滤波器对直方图进行平滑，将平滑直方图中存在的 主要高斯分布数设为聚类数(设为幻，直方图的各个极大值点作为聚类中心。
3 2DFCM迭代，迭代的目标函数为公式（4),约束条件力^ =1 。
./='
执行以下步骤直至聚类中心在两次迭代前后的改变量小于一个很小的常数eK).OOl Stepl)用公式(6)更新隶属度矩阵;v Step2)用公式(7)更新聚类中心^.，
Step3)计算更新前后聚类中心的改变量llM^-M。/』2。这里M。/rf表示更新前的M， M,表 示更新后的M。
4 2DHFCM初始化
将2DFCM收敛后的隶属度矩阵以及聚类中心作为2DHFCM的初始值，并用公式(20) 计算初始的亮度补偿值。以滤波数据f ={《.|/=1,2,〜,"}计算一维直方图，然后使用方差为 1.0的一维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为期 望的分类数(设为c)。
5 2DHFCM迭代,迭代的目标函数为公式(4)，约束条件尤^ =1 。
执行以下步骤直至迭代后存在的不同聚类中心数量小于期望的分类数c。 Stepl)用公式(18)更新隶属度矩阵^。
Step2)用公式(19)更新聚类中心mj和。
Step3)用公式(20)更新两个亮度补偿向量/?,、 y,。
Step4)如果聚类中心更新后，两聚类中心间的欧氏距离小于V^，则将两聚类中心合并。统 计不同聚类中心的数量。 实验结果分析如下：
对本发明的上述方法，进行了仿真测试。仿真测试的条件如下： 仿真平台为主频为2.0GHz的PC机，编程工具为Visual C++6.0。为验证算法在各种情
况下的性能，分别采用了一幅仅受斑点噪声污染的模拟图，和一幅既有斑点噪声又有亮度不均匀的模拟图。产生模拟图的方法为采用Jensen等开发的仿真软件Fidd一II。对模拟图 分割的有效性用分割准确度来衡量，分割准确度定义如下：
^ = ^C謹CT xl00% (21)
J v ro皿
其中AU^表示正确分类的像素点、数andA^M，^刚象^^③戶;?剤象m点激。g l(a)为包
含一个圆形目标的理想图像，图l(b)为用Fidd一II产生的受斑点噪声污染的模拟图像，图 1(c)-(g)分别给出了 SRAD的滤波结果，以及FCM，文献[8]， 2DFCM和2DHFCM的分割 结果。表1分别给出了 FCM，文献[8]， 2DFCM和2DHFCM分割结果的SA值。由结果可 见2DFCM和2DHFCM在仅包含斑点噪声的情况下分割结果相似，且明显优于FCM和文 献[8]的分割结果。
表1四种方法对仅含斑点噪声的模拟图的分割准确度对比
FC KFCM一 FCM2D HFCM2D
M S2[8]
Fig.2 92.3 96.26 98.88 98.64
(c=2) 2
下面对图l(b)加入亮度不均匀的情况，如图2(a)所示。图2(b)为该不均匀亮度的模型， 图2(c)-(h)分别给出了 FCM，文献[8]， 2DFCM和2DHFCM的分割结果，以及对原始图像 和对SRAD滤波图像的亮度补偿图。表2分别给出了四种方法分割结果的SA值。由结果 可见2DFCM在含有亮度不均匀的情况下算法失效，2DHFCM获得了较高的分割准确度。
表2四种方法对既含斑点噪声又含亮度不均匀的模拟图的分割准确度对比
FC KFCM— FCM2D HFCM2D
M
Fig.2 66.1 69.41 71.79 98.56
(c=2) 7
下面给出算法在实际超声图像分割中的应用实例。图3(a)为一幅经腹部超声得到的胎 儿肾积水图像，图3(b)-(f)分别给出了 SRAD滤波的结果，2DFCM和2DHFCM的分割结 果，以及对原始图像及对SRAD滤波图像的亮度补偿图。为了说明2DHFCM的聚类过程，图4给出了 2DHFCM在迭代过程中不同聚类数情况下的中间结果。图4禾口图5给出了对 一幅具有较复杂形状的胎儿主动脉弓超声图像的分割结果对比及聚类过程说明。图6给出 了对一幅图像质量较差且目标背景间对比度较低的乳腺图像分割结果。
由对模拟超声图像和实际超声图像的分割结果可见，本发明能够同时抑制超声图像中 斑点噪声和亮度不均匀对分割结果的影响，且分割中不需要人为调整经验参数，是一种有 效的超声图像分割新方法。
参考文献
[I] R. C. Gonzalez, R. E. Woods. Digital image processing [M]. Beijing: Publishing House of Electronics Industry, 2003. 209.
[2] ML Karaman， M. A. Kutay， G. Bozdagi. An adaptive speckle suppression filter for medical
ultrasonic imaging [J]. IEEE Trans. Medical Image, 1995, 14(2) : 283-292. [3] L. D. David. De-noising by soft-thresholding [J]. IEEE Trans. Information Theory, 1995, 41(3): 613-627.
[4] YU Yong-jian and T. A. Scott. Speckle reducing anisotropic diffusion [J]. IEEE Trans. Image
Processing, 2002， 11(11): 1260-1270. [5] D. L. Pham, J. L. Prince, "Adaptive fuzzy segmentation of magnetic resonance images,"
鹏£ 7>認Me"廳g., vol. 18, no. 9, pp. 737-752， Sept 1999. [6] G. R Xiao， M. Brady, J. A. Noble, and Yongyue Zhang， "Segmentation of ultrasound B-mode
images with intensity inhomogeneity correction," Thms*. /應g., vol. 21, no. 1,
pp.48-57, Jan. 2002
[7] J, C. Bezdek, "A convergence theorem for the fuzzy ISODATA clustering algorithm,"
Trara. M?c/u7?e /她//.， vol. PAMI-2, pp. 1-8, 1980.
[8] Chen Song匿can， Zhang Dao-qiang, "Robust Image Segmentation Using FCM With Spatial
Constraints Based on New Kernel-Induced Distance Measure," rram*. S;ast^勝，Mm.
Q/Z^w" vol. 34, no.4, pp. 1907-1916, Apr.2004. [9] YU Yong-jian and T. A. Scott. Edge detection in ultrasound imagery using the instantaneous
coefficient of variation [J]. IEEE Trans. Image Processing, 2004， 13(12): 1640-1655. [10] M. S. Hassouna, A. A. Farag, S. Hushek and T. Moriarty, "Cerebrovascular segmentation
fromTOF using stochastic models," Medical Image Analysis., vol. 10, pp. 2-18, 2006
[II] F. Kurugollu, B. Sankur， A. E. Harmanci, "Color image segmentation using histogram multithresholding and fusion," /腿ge a/W版o" C鐘p她7?g., vol. 19， pp. 915-928, 2001.
图像分割是指按某个特征如亮度或纹理，将图像划分为几个区域的过程。在医学图像 分析领域，图像分割是其它应用，如组织测量、解剖结构分析和组织定征的前提条件和关 键步骤。尽管已提出很多医学图像分割的方法，但其中大多数方法是对信噪比相对较高的 核磁共振图像或CT图像提出的。超声诊断因其安全、方便、无痛、价廉等优点己成为临 床不可缺少的诊断方法之一，但是由于超声图像中大量噪声的存在，使超声图像分割仍然 是一个悬而未决的难题。
在超声成像中，当人体组织的结构尺寸与入射超声波波长相近或小于波长时，超声束 发生散射。相位不同的散射回波相互干涉产生斑点噪声。作为超声图像中最为显著的噪声， 斑点噪声的存在降低了超声图像的质量，掩盖或模糊了图像中对诊断非常重要的细节信 息。因此，近年来众多学者陆续提出了维纳滤波[1]、自适应加权中值滤波[2]、小波软阈值 去噪[3]、各向异性扩散[4]等抑制斑点噪声的方法。除斑点噪声外，超声图像中由信号衰落 造成的亮度不均匀也是图像信噪比降低的主要原因。亮度不均匀使同一组织内部存在亮度 的变化，同时也降低了不同组织间的对比度，为超声图像的自动分割制造了更大的障碍。 因此，也提出过许多用于解决医学图像中的亮度不均匀问题[5'6]的方法，但很多方法仍然是 针对核磁共振图像提出的，为实现超声图像的自动分割，必须同时抑制斑点噪声和亮度不 均匀对分割结果的影响。 发明内容
本发明的目的是提出一种能同时抑制斑点超声和亮度不均匀对分割影响的超声图像自 动分割方法。
本发明提出的超声图像分割方法是一种基于模糊逻辑并融合斑点噪声抑制和亮度补偿 的超声图像分割方法，其具体步骤如下：首先，提出二维模糊C均值(2DFCM)的概念，在 将邻域信息引入到聚类过程中的同时，提高算法灵活性和收敛速度；接着，针对斑点噪声 的抑制，采用各向异性扩散滤波器SRAD为2DFCM提供邻域信息，提高算法对斑点噪声
的鲁棒性；最后，在将由信号衰减引起的亮度不均匀视为乘性噪声的基础上，根据超声图像获取过程中信号处理的特点，修改2DFCM的目标函数，引入亮度补偿因子，构造匀质 化的二维模糊C均值聚类(2DHFCM)。
下面对各步骤作进一步具体描述。
相关概念：模糊C均值聚类(FCM: Fuzzy C-Means)171。
若将图像像素值组成"个样本的集合义Hx,l /=1,2/-,"}，那么图像的分割问题就转化为 将这"个样本分成c个聚类的问题。设M^—户1,...,}为每个聚类的中心，〜是第/个样 本对第乂类的隶属度函数，则用隶属度函数定义的聚类损失函数JV^，可写成：
^鼎=tl>,i,『 （1)
其中6>1，是一个可控制聚类结果模糊程度的常数， 一般情况下6取2。图像的优化分割， 就是通过迭代寻找聚类中心,与隶属度值〜，使目标函数Ao/取最小。在t〜-i的限制
条件下，根据拉格朗日条件极值法，可得到目标函数/K^取条件极小值的必要条件：
(Zll ||2》/(*-"
=,广m』J (2)
1>a n、
〜=^- (3)
传统FCM算法在进行图像分割时，仅利用了像素的灰度信息，因此聚类结果对噪声 较为敏感。在FCM的目标函数中加入邻域约束是提高FCM抗噪性能的有效手段^。所谓 邻域约束就是在计算当前像素对各类的隶属情况时，将与目标像素相邻的各像素的隶属情 况作为约束条件，其依据是自然图像相邻像素存在的高相关性。邻域约束的加入使聚类结 果保持了图像空间上的连续性，抑制了噪声对聚类结果的影响。构造邻域约束的方法可参 考文献[8]。
1、 2DFCM
假设原始图像数据表示为^={^|/=1,2,'","}，原始图像经过滤波器后得到的数据表示为 i={?,h'=i,2，一,W。图像的分割问题就是要将二维样本《Z, f》分成c个聚类，每个聚类由
一个二维向量表征^={( ,，^)l产l,…,c)。这里^为f的聚类中心，产l,…,c由于原始图 像数据与其滤波数据之间的强相关性，可以用g的分类结果来影响其对应像素x,的分类结
果，即当x,对以".为中心的第y'类具有较高的隶属度时，《应该对以^为中心的第j'类具
有相似的高隶属度。3?,的分类对x,分类的影响程度可以通过一个加权常数a控制。根据经验，a—般取6-8。根据以上分析，本发明提出的2dfcm的目标函数可写为：
入为由Lagrange乘子法引入的中间变量。将目标函数分别对〜，〜和^求一阶导数并赋 值为零，可得到目标函数取极小值的迭代公式：
二维模糊聚类结构的提出较传统fcm具有以下优势：邻域信息引入到2dfcm的目标 函数中，增强了聚类对噪声的鲁棒性；在二维方向上寻找聚类中心，有潜力提高算法的收 敛速度；结构灵活，可以针对不同类型噪声，采用不同的滤波器提供邻域约束信息。
为抑制斑点噪声对聚类结果的影响，本发明采用各向异性扩散滤波器sradw来提供 2dfcm所需的邻域信息。各向异性扩散的基本原理是在抑制边缘附近区域平滑的同时， 加强平行于边缘方向的平滑。在srad中，滤波是一种扩散的过程，扩散的力度由下式决 定：
其中《(x，;;;0称为瞬时变化系数，它充当了扩散过程中的边缘检测器；的(^称为噪声比例函
数，它充当了扩散门限。扩散系数c(《)的构造应该是在均匀区域加强平滑，而在垂直于边
缘的方向抑制平滑。因此，srad在实际应用中的效果还取决于如下两个因素： 一是噪声
图像中边缘检测的有效性。即：《(x,少;/)能否有效区分由噪声或者纹理造成的小瞬时变化系
数和由边缘造成的大瞬时变化系数；二是扩散门限值估计的准确性。目卩：的W能否准确判
别何时加强或者抑制扩散过程。文献[4]提出的方法很好地解决了第一个问题，并在文献[9]
中进行了扩展，直接将其应用到超声图像的边缘检测中。对于第二个问题，文献[4]中仍然依靠经验选取初始扩散门限。
为给2DFCM提供准确的邻域约束信息，本发明使用二状态的瑞利、高斯混合分布对 超声图像灰度分布进行拟合，并采用EM算法实现混合分布的分解；根据分解结果预测图 像中斑点噪声均匀分布的区域；通过对均匀区域统计特性的分析获取各向异性扩散的扩散 参数。二状态的瑞利、高斯混合分布如下：
二 c^p(x I + o:2/3(x I Gotw."") (9)
其中a!和g2为混合系数，满足cq+a产l，且有：
1
2;rcr
exp
2-2 2a2
(10)
(11)
EM混合分布分解的方法可参考文献[IO]。当EM算法收敛后，根据似然分割准则[1()]，得 到目标区域、背景区域的判决门限：
2[1ogC57cr2) + C^/a2]/1/〃2 -l/o":
(a,
/2[log(cr，） + «2Zai]/|lAT2 -1//?2
得到判决门限后将图像分割为目标区域和背景区域(噪声区域)。然后将图像划分为30x30 的小块，选择噪声分布最均匀的块，即包含噪声点最多的块，根据下式计算初始扩散门限：
va柳)]
(13)
其中z(O表示斑点噪声均匀分布的区域，vw[z(0]和^分别表示此区域的方差和均值； 2、 2DHFCM
在2DFCM中本发明通过引入各向异性扩散作为邻域约束来提高聚类对斑点噪声的鲁 棒性，为增强算法对实际超声图像的适应性，还需要考虑解决由信号衰减所造成的亮度不 均匀问题。
亮度不均匀是在各类医学图像中存在的一个普遍问题，如核磁共振图像。迄今为止解 决亮度不均匀问题较成功的方法是将其假设为一种乘性噪声[5'6]，然后利用对数压縮将其转 变为加性噪声以简化后继的处理过程。假设J^(x^2,…,x4， 1^{力,》,...具}分别表示理想 背散射信号数据和实际背散射信号数据，这里理想数据指不包含信号衰减的数据，则理想 数据和实际数据间的关系可表示如下：
"=JT, xS,.， ls^" (14)
超声图像获取时包含了对数压縮、低通滤波、插值等操作，所以通常将观测超声图像灰度信号和实际背散射信号的关系简化为对数压縮关系。根据图像灰度信号与实际背散射信号 的对数压縮关系，则可将理想图像像素和观测图像像素的关系表示如下：
V/e{l,2,..."} (15)
其中5={艮艮…,y^即为对观测像素P(",》,…,W的亮度补偿。因为各向异性扩散滤波 器对观测超声图像具有均值保留特性，因此对各向异性扩散滤波后数据的亮度补偿可以采 用与原始图像相同的方法。设经过各向异性扩散后的数据？ = Jj具有亮度补偿
^={",72, ...Jn}，则F和i?的关系表示如下：
在2DFCM的目标函数中引入亮度补偿，即将式(4)中的x,用乂-A代替，3f,用义.-》代替， 就得到2DHFCM的目标函数：
在约束条件|;〜=1下，目标函数(17)的求解可采用与2DFCM相同的方法，得到目标函数
>1
取极小值的迭代公式：
在2DHFCM算法迭代前，需要预先决定的初始参数为聚类数、初始聚类中心和初始亮 度补偿估计。由于2DHFCM在每次参数迭代时需要更新三组相互影响的参数，如果初始 参数不准确则会由误差的传递造成算法的错误迭代，因此2DHFCM的正确收敛依赖于初 始参数的准确性。这里利用自底向上松弛标注法和高斯混合分布模型为2DHFCM提供初 始参数估计和保证算法的正确收敛。自底向上松弛标注法是在聚类中引入空间纹理信息的 方法，该方法通常包括两个处理过程：首先通过一种聚类算法如K均值法来预测初始的标 注概率，然后通过重新调整分类将空间连续性引入最终的分类结果中。这里先用2DFCM对图像进行预分类，然后将其分类的结果作为2DHFCM的初始状态；接着再用2DHFCM 重新调整分类的结果，达到同时消除斑点噪声和亮度不均匀的目的。在执行2DFCM前同 样需要给定聚类数及初始聚类中心，本发明中采用高斯混合分布模型来估计初始参数。即 将原始数据Z和i构成的二维直方图视为多个二维高斯分布混合而成，混合分布的个数即 为初始聚类数，混合分布的中心位置坐标即为初始聚类中心，实现过程可参考文献[ll]。 附图说明
图1、对包含斑点噪声超声图像分割结果的比较。其中图l(a)为理想图像，图l(b)为 加入斑点噪声的模拟图像，图l(c)为SRAD的滤波结果，图l(d)为FCM分割的结果，图 1(e)为文献[8]分割的结果，图l(f)为2DFCM的分割结果，图l(g)为2DHFCM的分割结果。
图2、对包含斑点噪声和亮度不均匀超声图像分割结果的比较。其中图2(a)为加入斑 点噪声和亮度不均匀的模拟图像，图2(b)为亮度不均匀模型，图2(c)为FCM分割的结果， 图2(d)为文献[8]分割的结果，图2(e)为2DFCM的分割结果，图2(f)为2DHFCM的分割结 果，图2(g)为对原始图像的亮度补偿图，图2(h)为对SRAD滤波图像的亮度补偿图。
图3、对一幅胎儿肾积水超声图像分割结果的比较。其中图3(a)为原始图像，图3(b) 为SRAD滤波的结果，图3(c)为2DFCM的分割结果，图3(d)为2DHFCM的分割结果，图 3(e)为对原始图像的亮度补偿图，图3(f)为对SRAD滤波图像的亮度补偿图。
图4、对图3的聚类过程说明。图4(a)为聚类数为30的结果，图4(b)为聚类数为10 的结果，图4(c)为聚类数为6的结果，图4(d)为聚类数为3的结果，图4(e)为聚类数为2 的结果。
图5、对一幅胎儿主动脉弓超声图像分割结果的比较。其中图5(a)为原始图像，图5(b) 为SRAD滤波的结果，图5(c)为2DFCM的分割结果，图5(d)为2DHFCM的分割结果，图 5(e)为对原始图像的亮度补偿图，图5(f)为对SRAD滤波图像的亮度补偿图。
图6、对图5的聚类过程说明。图6(a)为聚类数为30的结果，图6(b)为聚类数为10 的结果，图6(c)为聚类数为6的结果，图6(d)为聚类数为4的结果，图6(e)为聚类数为2 的结果。
图7、对一幅乳腺肿瘤超声图像分割结果的比较。其中图7(a)为原始图像，图7(b)为 2DHFCM的分割结果，图7(c)为对SRAD滤波图像的亮度补偿图，图7(d)为对分割图像提 取边缘轮廓的结果。 具体实施方式
以下通过整个算法的实现步骤，进一步对发明中的各个组成部分加以描述 1对原始图像进行SRAD滤波在滤波前需采用EM算法对图像的概率分布进行拟合，并在此基础上计算图像目标背 景分割的门限，得到分割图，然后通过加窗的方法找到图像中包含背景最多的区域，根据 公式(13)计算初始扩散门限。再通过SRAD得到邻域约束数据f ={《U=l,2,*",/7}。
2采用高斯混合分布模型估计2DFCM的初始参数
首先以根据原始数据;=1,2,'","}和滤波数据》={3f,|片1,2,…,W计算二维直方 图，然后使用方差为0.625的二维高斯滤波器对直方图进行平滑，将平滑直方图中存在的 主要高斯分布数设为聚类数(设为幻，直方图的各个极大值点作为聚类中心。
3 2DFCM迭代，迭代的目标函数为公式（4),约束条件力^ =1 。
./='
执行以下步骤直至聚类中心在两次迭代前后的改变量小于一个很小的常数eK).OOl Stepl)用公式(6)更新隶属度矩阵;v Step2)用公式(7)更新聚类中心^.，
Step3)计算更新前后聚类中心的改变量llM^-M。/』2。这里M。/rf表示更新前的M， M,表 示更新后的M。
4 2DHFCM初始化
将2DFCM收敛后的隶属度矩阵以及聚类中心作为2DHFCM的初始值，并用公式(20) 计算初始的亮度补偿值。以滤波数据f ={《.|/=1,2,〜,"}计算一维直方图，然后使用方差为 1.0的一维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为期 望的分类数(设为c)。
5 2DHFCM迭代,迭代的目标函数为公式(4)，约束条件尤^ =1 。
执行以下步骤直至迭代后存在的不同聚类中心数量小于期望的分类数c。 Stepl)用公式(18)更新隶属度矩阵^。
Step2)用公式(19)更新聚类中心mj和。
Step3)用公式(20)更新两个亮度补偿向量/?,、 y,。
Step4)如果聚类中心更新后，两聚类中心间的欧氏距离小于V^，则将两聚类中心合并。统 计不同聚类中心的数量。 实验结果分析如下：
对本发明的上述方法，进行了仿真测试。仿真测试的条件如下： 仿真平台为主频为2.0GHz的PC机，编程工具为Visual C++6.0。为验证算法在各种情
况下的性能，分别采用了一幅仅受斑点噪声污染的模拟图，和一幅既有斑点噪声又有亮度不均匀的模拟图。产生模拟图的方法为采用Jensen等开发的仿真软件Fidd一II。对模拟图 分割的有效性用分割准确度来衡量，分割准确度定义如下：
^ = ^C謹CT xl00% (21)
J v ro皿
其中AU^表示正确分类的像素点、数andA^M，^刚象^^③戶;?剤象m点激。g l(a)为包
含一个圆形目标的理想图像，图l(b)为用Fidd一II产生的受斑点噪声污染的模拟图像，图 1(c)-(g)分别给出了 SRAD的滤波结果，以及FCM，文献[8]， 2DFCM和2DHFCM的分割 结果。表1分别给出了 FCM，文献[8]， 2DFCM和2DHFCM分割结果的SA值。由结果可 见2DFCM和2DHFCM在仅包含斑点噪声的情况下分割结果相似，且明显优于FCM和文 献[8]的分割结果。
表1四种方法对仅含斑点噪声的模拟图的分割准确度对比
FC KFCM一 FCM2D HFCM2D
M S2[8]
Fig.2 92.3 96.26 98.88 98.64
(c=2) 2
下面对图l(b)加入亮度不均匀的情况，如图2(a)所示。图2(b)为该不均匀亮度的模型， 图2(c)-(h)分别给出了 FCM，文献[8]， 2DFCM和2DHFCM的分割结果，以及对原始图像 和对SRAD滤波图像的亮度补偿图。表2分别给出了四种方法分割结果的SA值。由结果 可见2DFCM在含有亮度不均匀的情况下算法失效，2DHFCM获得了较高的分割准确度。
表2四种方法对既含斑点噪声又含亮度不均匀的模拟图的分割准确度对比
FC KFCM— FCM2D HFCM2D
M
Fig.2 66.1 69.41 71.79 98.56
(c=2) 7
下面给出算法在实际超声图像分割中的应用实例。图3(a)为一幅经腹部超声得到的胎 儿肾积水图像，图3(b)-(f)分别给出了 SRAD滤波的结果，2DFCM和2DHFCM的分割结 果，以及对原始图像及对SRAD滤波图像的亮度补偿图。为了说明2DHFCM的聚类过程，图4给出了 2DHFCM在迭代过程中不同聚类数情况下的中间结果。图4禾口图5给出了对 一幅具有较复杂形状的胎儿主动脉弓超声图像的分割结果对比及聚类过程说明。图6给出 了对一幅图像质量较差且目标背景间对比度较低的乳腺图像分割结果。
由对模拟超声图像和实际超声图像的分割结果可见，本发明能够同时抑制超声图像中 斑点噪声和亮度不均匀对分割结果的影响，且分割中不需要人为调整经验参数，是一种有 效的超声图像分割新方法。
参考文献
[I] R. C. Gonzalez, R. E. Woods. Digital image processing [M]. Beijing: Publishing House of Electronics Industry, 2003. 209.
[2] ML Karaman， M. A. Kutay， G. Bozdagi. An adaptive speckle suppression filter for medical
ultrasonic imaging [J]. IEEE Trans. Medical Image, 1995, 14(2) : 283-292. [3] L. D. David. De-noising by soft-thresholding [J]. IEEE Trans. Information Theory, 1995, 41(3): 613-627.
[4] YU Yong-jian and T. A. Scott. Speckle reducing anisotropic diffusion [J]. IEEE Trans. Image
Processing, 2002， 11(11): 1260-1270. [5] D. L. Pham, J. L. Prince, "Adaptive fuzzy segmentation of magnetic resonance images,"
鹏£ 7>認Me"廳g., vol. 18, no. 9, pp. 737-752， Sept 1999. [6] G. R Xiao， M. Brady, J. A. Noble, and Yongyue Zhang， "Segmentation of ultrasound B-mode
images with intensity inhomogeneity correction," Thms*. /應g., vol. 21, no. 1,
pp.48-57, Jan. 2002
[7] J, C. Bezdek, "A convergence theorem for the fuzzy ISODATA clustering algorithm,"
Trara. M?c/u7?e /她//.， vol. PAMI-2, pp. 1-8, 1980.
[8] Chen Song匿can， Zhang Dao-qiang, "Robust Image Segmentation Using FCM With Spatial
Constraints Based on New Kernel-Induced Distance Measure," rram*. S;ast^勝，Mm.
Q/Z^w" vol. 34, no.4, pp. 1907-1916, Apr.2004. [9] YU Yong-jian and T. A. Scott. Edge detection in ultrasound imagery using the instantaneous
coefficient of variation [J]. IEEE Trans. Image Processing, 2004， 13(12): 1640-1655. [10] M. S. Hassouna, A. A. Farag, S. Hushek and T. Moriarty, "Cerebrovascular segmentation
fromTOF using stochastic models," Medical Image Analysis., vol. 10, pp. 2-18, 2006
[II] F. Kurugollu, B. Sankur， A. E. Harmanci, "Color image segmentation using histogram multithresholding and fusion," /腿ge a/W版o" C鐘p她7?g., vol. 19， pp. 915-928, 2001.