CN1924926A - 一种基于二维模糊聚类的超声图像分割方法 - Google Patents

Abstract

本发明属于超声图像处理技术领域，具体为一种二维模糊聚类结合斑点噪声滤波器和亮度补偿的B型超声图像分割方法。该方法将基于图像亮度信息的传统模糊聚类方法扩展为同时在像素和像素邻域值进行模糊更新的二维模糊聚类；在二维模糊聚类目标函数中引入由各向异性扩散斑点噪声滤波器提供的邻域信息，加强二维模糊聚类对斑点噪声的鲁棒性；将亮度不均匀假设为乘性噪声的基础上，通过在二维模糊聚类目标函数中引入亮度补偿因子，利用图像中目标和背景存在的均匀性差异，加强二维模糊聚类对亮度不均匀噪声的鲁棒性。本发明能够在分割图像中同时抑制斑点噪声和亮度不均匀的影响，提供一种超声图像分割的新的有效方法。

Description

一种基于二维模糊聚类的超声图像分割方法

技术领域

本发明属于超声图像处理技术领域，具体涉及一种二维模糊聚类结合斑点噪声滤波器和亮度补偿的超声图像分割方法。

技术背景

图像分割是指按某个特征如亮度或纹理，将图像划分为几个区域的过程。在医学图像分析领域，图像分割是其它应用，如组织测量、解剖结构分析和组织定征的前提条件和关键步骤。尽管已提出很多医学图像分割的方法，但其中大多数方法是对信噪比相对较高的核磁共振图像或CT图像提出的。超声诊断因其安全、方便、无痛、价廉等优点已成为临床不可缺少的诊断方法之一，但是由于超声图像中大量噪声的存在，使超声图像分割仍然是一个悬而未决的难题。

在超声成像中，当人体组织的结构尺寸与入射超声波波长相近或小于波长时，超声束发生散射。相位不同的散射回波相互干涉产生斑点噪声。作为超声图像中最为显著的噪声，斑点噪声的存在降低了超声图像的质量，掩盖或模糊了图像中对诊断非常重要的细节信息。因此，近年来众多学者陆续提出了维纳滤波^[1]、自适应加权中值滤波^[2]、小波软阈值去噪^[3]、各向异性扩散^[4]等抑制斑点噪声的方法。除斑点噪声外，超声图像中由信号衰落造成的亮度不均匀也是图像信噪比降低的主要原因。亮度不均匀使同一组织内部存在亮度的变化，同时也降低了不同组织间的对比度，为超声图像的自动分割制造了更大的障碍。因此，也提出过许多用于解决医学图像中的亮度不均匀问题^[5，6]的方法，但很多方法仍然是针对核磁共振图像提出的，为实现超声图像的自动分割，必须同时抑制斑点噪声和亮度不均匀对分割结果的影响。

发明内容

本发明的目的是提出一种能同时抑制斑点超声和亮度不均匀对分割影响的超声图像自动分割方法。

本发明提出的超声图像分割方法是一种基于模糊逻辑并融合斑点噪声抑制和亮度补偿的超声图像分割方法，其具体步骤如下：首先，提出二维模糊C均值(2DFCM)的概念，在将邻域信息引入到聚类过程中的同时，提高算法灵活性和收敛速度；接着，针对斑点噪声的抑制，采用各向异性扩散滤波器SRAD为2DFCM提供邻域信息，提高算法对斑点噪声的鲁棒性；最后，在将由信号衰减引起的亮度不均匀视为乘性噪声的基础上，根据超声图像获取过程中信号处理的特点，修改2DFCM的目标函数，引入亮度补偿因子，构造匀质化的二维模糊C均值聚类(2DHFCM)。

下面对各步骤作进一步具体描述。

相关概念：模糊C均值聚类(FCM：Fuzzy C-Means)^[7]。

若将图像像素值组成n个样本的集合X＝{x_i|i＝1，2，...，n}，那么图像的分割问题就转化为将这n个样本分成c个聚类的问题。设M＝{m_j|j＝1，...，c}为每个聚类的中心，μ_ij是第i个样本对第j类的隶属度函数，则用隶属度函数定义的聚类损失函数J_FCM，可写成：

$J_{\mathrm{FCM}}=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b{\left|\right|x_i-m_j\left|\right|}^2---\left(1\right)$

其中b＞1，是一个可控制聚类结果模糊程度的常数，一般情况下b取2。图像的优化分割，就是通过迭代寻找聚类中心m_j与隶属度值μ_ij，使目标函数J_FCM取最小。在 $\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}=1$ 的限制条件下，根据拉格朗日条件极值法，可得到目标函数J_FCM取条件极小值的必要条件：

${\mathrm{\μ}}_{\mathrm{ij}}=\frac{{(1/{\left|\right|x_i-m_j\left|\right|}^2)}^{1/(b-1)}}{\underset{k=1}{\overset{c}{\mathrm{\Σ}}}{(1/{\left|\right|x_i-m_k\left|\right|}^2)}^{1/(b-1)}}---\left(2\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b}---\left(3\right)$

传统FCM算法在进行图像分割时，仅利用了像素的灰度信息，因此聚类结果对噪声较为敏感。在FCM的目标函数中加入邻域约束是提高FCM抗噪性能的有效手段^[8]。所谓邻域约束就是在计算当前像素对各类的隶属情况时，将与目标像素相邻的各像素的隶属情况作为约束条件，其依据是自然图像相邻像素存在的高相关性。邻域约束的加入使聚类结果保持了图像空间上的连续性，抑制了噪声对聚类结果的影响。构造邻域约束的方法可参考文献[8]。

1、2DFCM

假设原始图像数据表示为X＝{x_i|i＝1，2，...，n}，原始图像经过滤波器后得到的数据表示为 $\tilde{X}=\left\{{\tilde{X}}_i\right|i=\mathrm{1,2},\·\·\·,n\}.$ 图像的分割问题就是要将二维样本

分成c个聚类，每个聚类由一个二维向量表征 $\stackrel{\‾}{M}=\left\{(m_j,{\tilde{m}}_j)\right|j=1,\·\·\·,c\}.$ 这里

为

的聚类中心，j＝1，...，c由于原始图像数据与其滤波数据之间的强相关性，可以用

的分类结果来影响其对应像素x_i的分类结果，即当x_i对以m_j为中心的第j类具有较高的隶属度时，

应该对以

为中心的第j类具有相似的高隶属度。 的分类对x_i分类的影响程度可以通过一个加权常数α控制。根据经验，α一般取6-8。根据以上分析，本发明提出的2DFCM的目标函数可写为：

$J_{2\mathrm{DFCM}}=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]---\left(4\right)$

在约束条件 $\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}=1$ 下，目标函数的求解可采用Lagrange(拉格朗日)乘子法，构造新的无约束条件的目标函数：

$F=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]+\mathrm{\λ}(1-\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}})---\left(5\right)$

λ为由Lagrange乘子法引入的中间变量。将目标函数分别对μ_ij，m_j和

求一阶导数并赋值为零，可得到目标函数取极小值的迭代公式：

${\mathrm{\μ}}_{\mathrm{ij}}=\frac{{[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{c}{\mathrm{\Σ}}}{[{(x_k-m_j)}^2+\mathrm{\α}{({\tilde{x}}_k-{\tilde{m}}_j)}^2]}^{-\frac{1}{(b-1)}}}---\left(6\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b},{\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{\tilde{x}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b}---\left(7\right)$

二维模糊聚类结构的提出较传统FCM具有以下优势：邻域信息引入到2DFCM的目标函数中，增强了聚类对噪声的鲁棒性；在二维方向上寻找聚类中心，有潜力提高算法的收敛速度；结构灵活，可以针对不同类型噪声，采用不同的滤波器提供邻域约束信息。

为抑制斑点噪声对聚类结果的影响，本发明采用各向异性扩散滤波器SRAD^[4]来提供2DFCM所需的邻域信息。各向异性扩散的基本原理是在抑制边缘附近区域平滑的同时，加强平行于边缘方向的平滑。在SRAD中，滤波是一种扩散的过程，扩散的力度由下式决定：

$c\left(q\right)=\frac{1}{1+[q^2(x,y;t)-q_0^2\left(t\right)]/\left[q_0^2\left(t\right)(1+q_0^2\left(t\right))\right]}---\left(8\right)$

其中q(x，y；t)称为瞬时变化系数，它充当了扩散过程中的边缘检测器；q₀(t)称为噪声比例函数，它充当了扩散门限。扩散系数c(q)的构造应该是在均匀区域加强平滑，而在垂直于边缘的方向抑制平滑。因此，SRAD在实际应用中的效果还取决于如下两个因素：一是噪声图像中边缘检测的有效性。即：q(x，y；t)能否有效区分由噪声或者纹理造成的小瞬时变化系数和由边缘造成的大瞬时变化系数；二是扩散门限值估计的准确性。即：q₀(t)能否准确判别何时加强或者抑制扩散过程。文献[4]提出的方法很好地解决了第一个问题，并在文献[9]中进行了扩展，直接将其应用到超声图像的边缘检测中。对于第二个问题，文献[4]中仍然依靠经验选取初始扩散门限。

为给2DFCM提供准确的邻域约束信息，本发明使用二状态的瑞利、高斯混合分布对超声图像灰度分布进行拟合，并采用EM算法实现混合分布的分解；根据分解结果预测图像中斑点噪声均匀分布的区域；通过对均匀区域统计特性的分析获取各向异性扩散的扩散参数。二状态的瑞利、高斯混合分布如下：

              p(x)＝α₁p(x|Rayleigh)+α₂p(x|Gaussian)                (9)

其中α₁和α₂为混合系数，满足α₁+α₂＝1，且有：

$p\left(x\right|\mathrm{Rayleigh})=\frac{x}{{\mathrm{\β}}^2}\exp (-\frac{x^2}{{2\mathrm{\β}}^2})---\left(10\right)$

$p\left(x\right|\mathrm{Gaussian})=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \left[\frac{-{(x-\mathrm{\μ})}^2}{{2\mathrm{\σ}}^2}\right]---\left(11\right)$

EM混合分布分解的方法可参考文献[10]。当EM算法收敛后，根据似然分割准则^[10]，得到目标区域、背景区域的判决门限：

$T=\left\{\begin{array}{cc}\sqrt{2[\log ({\mathrm{\&beta;}}^2/{\mathrm{\&sigma;}}^2)+{\mathrm{\&alpha;}}_1/{\mathrm{\&alpha;}}_2]/|1/{\mathrm{\&beta;}}^2-1/{\mathrm{\&sigma;}}^2|}& ({\mathrm{\&alpha;}}_1\&GreaterEqual;{\mathrm{\&alpha;}}_2)\\ \sqrt{2[\log ({\mathrm{\&sigma;}}^2/{\mathrm{\&beta;}}^2)+{\mathrm{\&alpha;}}_2/{\mathrm{\&alpha;}}_1]/|1/{\mathrm{\&sigma;}}^2-1/{\mathrm{\&beta;}}^2|}& ({\mathrm{\&alpha;}}_1<{\mathrm{\&alpha;}}_2)\end{array}\right.---\left(12\right)$

得到判决门限后将图像分割为目标区域和背景区域(噪声区域)。然后将图像划分为30×30的小块，选择噪声分布最均匀的块，即包含噪声点最多的块，根据下式计算初始扩散门限：

$q_0\left(t\right)=\frac{\sqrt{\mathrm{var}\left[z\left(t\right)\right]}}{\stackrel{\&OverBar;}{z\left(t\right)}}---\left(13\right)$

其中z(t)表示斑点噪声均匀分布的区域，var[z(t)]和 z(t)分别表示此区域的方差和均值；

2、2DHFCM

在2DFCM中本发明通过引入各向异性扩散作为邻域约束来提高聚类对斑点噪声的鲁棒性，为增强算法对实际超声图像的适应性，还需要考虑解决由信号衰减所造成的亮度不均匀问题。

亮度不均匀是在各类医学图像中存在的一个普遍问题，如核磁共振图像。迄今为止解决亮度不均匀问题较成功的方法是将其假设为一种乘性噪声^[5，6]，然后利用对数压缩将其转变为加性噪声以简化后继的处理过程。假设X＝{x₁，x₂，...，x_n}，Y＝{y₁，y₂，...，y_n}分别表示理想背散射信号数据和实际背散射信号数据，这里理想数据指不包含信号衰减的数据，则理想数据和实际数据间的关系可表示如下：

                           Y_i＝X_i×B_i，1≤i≤n                (14)

超声图像获取时包含了对数压缩、低通滤波、插值等操作，所以通常将观测超声图像灰度信号和实际背散射信号的关系简化为对数压缩关系。根据图像灰度信号与实际背散射信号的对数压缩关系，则可将理想图像像素和观测图像像素的关系表示如下：

               y_i＝x_i+β_i   i∈{1，2，...n}      (15)

其中B＝{β₁，β₂，...，β_n}即为对观测像素Y＝{y₁，y₂，...，y_n}的亮度补偿。因为各向异性扩散滤波器对观测超声图像具有均值保留特性，因此对各向异性扩散滤波后数据的亮度补偿可以采用与原始图像相同的方法。设经过各向异性扩散后的数据 $\tilde{Y}=\{{\tilde{y}}_1,{\tilde{y}}_2,\&CenterDot;\&CenterDot;\&CenterDot;{\tilde{y}}_N\}$ 具有亮度补偿R＝{γ₁，γ₂，...，γ_n}，则

和R的关系表示如下：

${\tilde{y}}_i={\tilde{x}}_i+{\mathrm{\&gamma;}}_i,\&ForAll;i\&Element;\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;n\}---\left(16\right)$

在2DFCM的目标函数中引入亮度补偿，即将式(4)中的x_i用y_i-β_i代替， 用

代替，就得到2DHFCM的目标函数：

$J_{2\mathrm{DHFCM}}=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]---\left(17\right)$

在约束条件 $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1$ 下，目标函数(17)的求解可采用与2DFCM相同的方法，得到目标函数取极小值的迭代公式：

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{-\frac{1}{(b-1)}}}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{[{(y_k-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_k-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(18\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b(y_i-{\mathrm{\&beta;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},{\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b({\tilde{y}}_i-{\mathrm{\&gamma;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(19\right)$

${\mathrm{\&beta;}}_i=y_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{m}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b},{\mathrm{\&gamma;}}_i={\tilde{y}}_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{\tilde{m}}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b}---\left(20\right)$

在2DHFCM算法迭代前，需要预先决定的初始参数为聚类数、初始聚类中心和初始亮度补偿估计。由于2DHFCM在每次参数迭代时需要更新三组相互影响的参数，如果初始参数不准确则会由误差的传递造成算法的错误迭代，因此2DHFCM的正确收敛依赖于初始参数的准确性。这里利用自底向上松弛标注法和高斯混合分布模型为2DHFCM提供初始参数估计和保证算法的正确收敛。自底向上松弛标注法是在聚类中引入空间纹理信息的方法，该方法通常包括两个处理过程：首先通过一种聚类算法如K均值法来预测初始的标注概率，然后通过重新调整分类将空间连续性引入最终的分类结果中。这里先用2DFCM对图像进行预分类，然后将其分类的结果作为2DHFCM的初始状态；接着再用2DHFCM重新调整分类的结果，达到同时消除斑点噪声和亮度不均匀的目的。在执行2DFCM前同样需要给定聚类数及初始聚类中心，本发明中采用高斯混合分布模型来估计初始参数。即将原始数据X和 构成的二维直方图视为多个二维高斯分布混合而成，混合分布的个数即为初始聚类数，混合分布的中心位置坐标即为初始聚类中心，实现过程可参考文献[11]。

附图说明

图1、对包含斑点噪声超声图像分割结果的比较。其中图1(a)为理想图像，图1(b)为加入斑点噪声的模拟图像，图1(c)为SRAD的滤波结果，图1(d)为FCM分割的结果，图1(e)为文献[8]分割的结果，图1(f)为2DFCM的分割结果，图1(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得到邻域约束数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}.$

2 采用高斯混合分布模型估计2DFCM的初始参数

首先以根据原始数据X＝{x_i|i＝1，2，…，n}和滤波数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ 计算二维直方图，然后使用方差为0.625的二维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为聚类数(设为K)，直方图的各个极大值点作为聚类中心。

3 2DFCM迭代，迭代的目标函数为公式(4)，约束条件 $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1.$

执行以下步骤直至聚类中心在两次迭代前后的改变量小于一个很小的常数ε＝0.001

Step1)用公式(6)更新隶属度矩阵μ_ij，

Step2)用公式(7)更新聚类中心m_j，

Step3)计算更新前后聚类中心的改变量||M_new-M_old||²。这里M_old表示更新前的M，M_new表示更新后的M。

4 2DHFCM初始化

将2DFCM收敛后的隶属度矩阵以及聚类中心作为2DHFCM的初始值，并用公式(20)计算初始的亮度补偿值。以滤波数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ 计算一维直方图，然后使用方差为1.0的一维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为期望的分类数(设为c)。

5 2DHFCM迭代，迭代的目标函数为公式(4)，约束条件 $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1.$

执行以下步骤直至迭代后存在的不同聚类中心数量小于期望的分类数c。

Step1)用公式(18)更新隶属度矩阵μ_ij。

Step2)用公式(19)更新聚类中心mj和

Step3)用公式(20)更新两个亮度补偿向量β_i、γ_i。

Step4)如果聚类中心更新后，两聚类中心间的欧氏距离小于，则将两聚类中心合并。统计不同聚类中心的数量。

实验结果分析如下：

对本发明的上述方法，进行了仿真测试。仿真测试的条件如下：

仿真平台为主频为2.0GHz的PC机，编程工具为Visual c++6.0。为验证算法在各种情况下的性能，分别采用了一幅仅受斑点噪声污染的模拟图，和一幅既有斑点噪声又有亮度不均匀的模拟图。产生模拟图的方法为采用Jensen等开发的仿真软件Field_II。对模拟图分割的有效性用分割准确度来衡量，分割准确度定义如下：

$\mathrm{SA}=\frac{N_{\mathrm{CORRECT}}}{N_{\mathrm{TOTAL}}}\&times;100\%---\left(21\right)$

其中N_CORRECT表示正确分类的像素点数andN_TOTAL表示图像包含的所有像素点数。图1(a)为包含一个圆形目标的理想图像，图1(b)为用Field_II产生的受斑点噪声污染的模拟图像，图1(c)-(g)分别给出了SRAD的滤波结果，以及FCM，文献[8]，2DFCM和2DHFCM的分割结果。表1分别给出了FCM，文献[8]，2DFCM和2DHFCM分割结果的SA值。由结果可见2DFCM和2DHFCM在仅包含斑点噪声的情况下分割结果相似，且明显优于FCM和文献[8]的分割结果。

   表1四种方法对仅含斑点噪声的模拟图的分割准确度对比

	FCM	KFCM_S2[8]	FCM2D	HFCM2D
					Fig.2(c＝2)	92.32	96.26	98.88	98.64

下面对图1(b)加入亮度不均匀的情况，如图2(a)所示。图2(b)为该不均匀亮度的模型，图2(c)-(h)分别给出了FCM，文献[8]，2DFCM和2DHFCM的分割结果，以及对原始图像和对SRAD滤波图像的亮度补偿图。表2分别给出了四种方法分割结果的SA值。由结果可见2DFCM在含有亮度不均匀的情况下算法失效，2DHFCM获得了较高的分割准确度。

表2四种方法对既含斑点噪声又含亮度不均匀的模拟图的分割准确度对比

	FCM	KFCM_S2[8]	FCM2D	HFCM2D
					Fig.2(c＝2)	66.17	69.41	71.79	98.56

下面给出算法在实际超声图像分割中的应用实例。图3(a)为一幅经腹部超声得到的胎儿肾积水图像，图3(b)-(f)分别给出了SRAD滤波的结果，2DFCM和2DHFCM的分割结果，以及对原始图像及对SRAD滤波图像的亮度补偿图。为了说明2DHFCM的聚类过程，图4给出了2DHFCM在迭代过程中不同聚类数情况下的中间结果。图4和图5给出了对一幅具有较复杂形状的胎儿主动脉弓超声图像的分割结果对比及聚类过程说明。图6给出了对一幅图像质量较差且目标背景间对比度较低的乳腺图像分割结果。

由对模拟超声图像和实际超声图像的分割结果可见，本发明能够同时抑制超声图像中斑点噪声和亮度不均匀对分割结果的影响，且分割中不需要人为调整经验参数，是一种有效的超声图像分割新方法。

参考文献

[1]R.C.Gonzalez，R.E.Woods.Digital image processing[M].Beijing：Publishing House ofElectronics Industry，2003.209.

[2]M.Karaman，M.A.Kutay，G.Bozdagi.An adaptive speckle suppression filter for medicalultrasonic 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.ImageProcessing，2002，11(11)：1260-1270.

[5]D.L.Pham，J.L.Prince，“Adaptive fuzzy segmentation of magnetic resonance images，”IEEE Trans Med Imag.，vol.18，no.9，pp.737-752，Sept.1999.

[6]G.F.Xiao，M.Brady，J.A.Noble，and Yongyue Zhang，“Segmentation of ultrasound B-modeimages with intensity inhomogeneity correction，”IEEE Trans.Med Imag.，vol.21，no.1，pp.48-57，Jan.2002

[7]J.C.Bezdek，“A convergence theorem for the fuzzy ISODATA clustering algorithm，”IEEETrans.Pattern Anal.Machine Intell.，vol.PAMI-2，pp.1-8，1980.

[8]Chen Song-can，Zhang Dao-qiang，“Robust Image Segmentation Using FCM With SpatialConstraints Based on New Kernel-Induced Distance Measure，”IEEE Trans.Systems，Man.Cybern.，vol.34，no.4，pp.1907-1916，Apr.2004.

[9]YU Yong-jian and T.A.Scott.Edge detection in ultrasound imagery using the instantaneouscoefficient 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 segmentationfrom TOF using stochastic models，”Medical Image Analysis.，vol.10，pp.2-18，2006

[11]F.Kurugollu，B.Sankur，A.E.Harmanci，“Color image segmentation using histogrammultithresholding and fusion，”Image and Vision Computing.，vol.19，pp.915-928，2001.

Claims (2)

1、一种基于模糊聚类的超声图像分割方法，其特征在于具体步骤为：首先，提出2DFCM的概念，将邻域信息引入到聚类过程中，以提高算法灵活性和收敛速度；接着，针对斑点噪声的抑制，采用各向异性扩散滤波器SRAD为2DFCM提供邻域信息，提高算法对斑点噪声的鲁棒性；最后，将由信号衰减引起的亮度不均匀视为乘性噪声，根据超声图像获取过程中信号处理的特点，修改2DFCM的目标函数，引入亮度补偿因子，构造匀质化的二维模糊C均值聚类，这里2DHFCM为二维模糊C均值。

2、根据权利要求1所述的基于模糊聚类的超声图像分割方法，其特征在于具体步骤为：

(1)对原始图像进行SRAD滤波

在滤波前需采用EM算法对图像的概率分布进行拟合，并在此基础上计算图像目标背景分割的门限，得到分割图，然后通过加窗的方法找到图像中包含背景最多的区域，根据公式(13)计算初始扩散门限，再通过SRAD得到邻域约束数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\};$

(2)采用高斯混合分布模型估计2DFCM的初始参数

首先根据原始数据X＝{x_i|i＝1，2，…，n}和滤波数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ 计算二维直方图，然后使用方差为0.625的二维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为聚类数K，直方图的各个极大值点作为聚类中心；

(3)2DFCM迭代，迭代的目标函数为公式(4)，约束条件为 $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1,$

执行以下步骤直至聚类中心在两次迭代前后的改变量小于一个很小的常数ε＝0.001，

(1)用公式(6)更新隶属度矩阵，μ_ij，

(2)用公式(7)更新聚类中心m_j，

(3)计算更新前后聚类中心的改变量‖M_new-M_old‖²，这里M_old表示更新前的M，M_new表示更新后的M；

(4)2DHFCM初始化，这里2DHFCM为匀质化的二维模糊C均值聚类，

将2DFCM收敛后的隶属度矩阵以及聚类中心作为2DHFCM的初始值，并用公式(20)计算初始的亮度补偿值；以滤波数据 $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ 计算一维直方图，然后使用方差为1.0的一维高斯滤波器对直方图进行平滑，将平滑直方图中存在的主要高斯分布数设为期望的分类数c；

(5)2DHFCM迭代，迭代的目标函数为公式(17)，约束条件为 $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1,$

执行以下步骤直至迭代后存在的不同聚类中心数量小于期望的分类数c，

(1)用公式(18)更新隶属度矩阵μ_ij，

(2)用公式(19)更新聚类中心m_j，

(3)用公式(20)更新两个亮度补偿向量，β_i、γ_i，

(4)如果聚类中心更新后，两聚类中心间的欧氏距离小于，则将两聚类中心合并，统计不同聚类中心的数量；

$J_{2\mathrm{DFCM}}=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]---\left(4\right)$

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(x_i-m_j)}^2+{\mathrm{\&alpha;}({\tilde{x}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{[{(x_k-m_j)}^2+\mathrm{\&alpha;}{({\tilde{x}}_k-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(6\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},$ ${\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b{\tilde{x}}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(7\right)$

$q_0\left(t\right)=\frac{\sqrt{\mathrm{var}\left[z\left(t\right)\right]}}{\stackrel{\&OverBar;}{z\left(t\right)}}---\left(13\right)$

$J_{2\mathrm{DHFCM}}=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]---\left(17\right)$

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{[{(y_k-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_k-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(18\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b(y_i-{\mathrm{\&beta;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},$ ${\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b({\tilde{y}}_i-{\mathrm{\&gamma;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(19\right)$

${\mathrm{\&beta;}}_i=y_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{m}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b},$ ${\mathrm{\&gamma;}}_i={\tilde{y}}_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{\tilde{m}}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b}---\left(20\right)$

这里，b＝2，α取6-8；{β₁，β₂，…，β_n}＝B为对观测像素Y＝{y₁，y₂，…，y_n}的亮度补偿， $\tilde{Y}=\{{\tilde{y}}_1,{\tilde{y}}_2,\&CenterDot;\&CenterDot;\&CenterDot;{\tilde{y}}_N\}$ 为各向异性扩散后的数据，{γ₁，γ₂，…，γ_n}＝R为

的亮度补偿，其中z(t)表示斑点噪声均匀分布的区域，var[z(t)]和 z(t)分别表示此区域的方差和均值。

Images