CN109410214B - 一种具有抗噪性和偏场校正的医学图像分割方法 - Google Patents

Abstract

本发明公开了一种具有抗噪性和偏场校正的医学图像分割方法，该方法包括以下步骤：初始化隶属度矩阵和多项式拟合曲面的系数，并设置聚类个数、隶属度权重因子以及终止阈值；通过目标函数分别对偏场、拟合曲面系数以及隶属度求偏导，得到偏场、拟合曲面系数以及隶属度的更新计算公式，通过计算更新偏场、拟合曲面系数以及隶属度求得最优解；根据隶属度最大的准则分类标记像素点，并将能量损失函数加入初步分割的图像中，完成带偏场的医学图像的偏场校正和最终分割。本发明的方法很大程度地消除了偏场带来的影响，能够取得令人满意的分割效果，同时加入了抗噪模糊因子并采用多项式拟合曲面代替聚类中心，能够更好的抵抗噪声影响和保持更精确的图像细节。

Description

一种具有抗噪性和偏场校正的医学图像分割方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种具有抗噪性和偏场校正的医学数字图像分割方法。

背景技术

图像是人类视觉的基础，是自然景物的客观反映，也是人类信息收集和交流的有效途径之一，随着计算机、互联网等技术的发展，数字图像处理技术在工业、医疗、军事、交通等行业中发挥着日益重要的作用。

为了有效地提取和利用数字图像中所包含的信息，需要对图像进行分割。图像分割是图像识别和计算机视觉至关重要的预处理，没有正确的分割就不可能有正确的识别。其意义是将一副图像分割成一组互不相交的子区域，并且在同一区域内部具有相同或相似的特性，这里的特征可以是灰度、颜色、纹理等。图像分割是模式识别、计算机视觉、图像理解等领域的最重要的基础环节，它是由图像处理到图像分析的关键步骤，也是一种基本的计算机视觉技术。但是，进行分割仅有的依据是图像中像素的亮度及颜色，由计算机自动处理分割时，将会遇到各种困难。例如，光照不均匀、噪声的影响、图像中存在不清晰的部分，以及阴影等，常常发生分割错误。因此图像分割是需要进一步研究的技术。

模糊C均值(Fuzzy C-means)算法简称FCM算法，是一种基于目标函数的模糊聚类算法，主要用于数据的聚类分析。该算法理论成熟，应用广泛，是一种优秀的聚类算法。该算法由两部分组成，一是隶属度，简单理解就是一种相似程度，是一个样本相似于不同结果的程度指标。二是某样本点到各个类中心的距离，简单理解就是该样本点到聚类中心的欧式距离。该算法由这两部分相乘得到目标函数，算法思想是先规定聚类中心的个数，如果某个样本点到聚类中心的距离小，则隶属度就大，通过求解目标函数，使结果达到最小值，算法在迭代求解目标函数最小值的过程逐渐趋于稳定，分割效果达到一个相对较好的聚类分割效果。FCM的目标函数E为：

其中，C是聚类个数，N是图像像素总数，u_ij表示的是第j个像素点相对于第i个聚类中心的模糊隶属度，m是隶属度的一个因子，即模糊加权系数，m>1，按照聚类有效性问题研究的结果，m的取值范围可以限制在1.5≤m≤2.5，所以取m＝2是最简洁的值，v_i是第i个聚类中心，||x_j-v_i||为第j个像素点与第i个聚类中心之间的欧氏距离。该算法中隶属度u_ij的取值范围是[0,1]，表示各像素点隶属于各聚类中心的程度，并且满足各像素属于不同聚类的模糊隶属度的和为1，即该目标函数能够满足最小化约束条件定义如下：

将最小化约束条件代入重新构建的目标函数中，得到目标函数取得最小值的必要条件：

其中，λ_j(1≤j≤n)为n个有约束的拉格朗日乘子。将上式同时对u_ij和v_i求导，获得目标函数最小化必要条件为：

和

根据上述两式迭代求解模糊隶属度u_ij以及聚类中心v_i，最终实现函数优化并取得最小值。

然而，标准的FCM方法在图像分割中存在以下几个明显的缺点：(1)初始聚类中心的选择对图像分割结果有很大影响，如果初始聚类中心选的好，那可以通过较少的迭代使能量函数达到最优，同时会得到比较好的分割结果，相反则会增加算法的迭代次数，同时也可能使算法陷入局部最优，得不到理想的分割结果；(2)FCM算法只考虑像素点到聚类中心的欧氏距离，没有将像素的空间信息考虑进去，从而导致方法对于噪声比较敏感，当处理有噪声的图像时，不能取得令人满意的分割结果。(3)由于磁共振成像设备的不完整性以及图像本身的特异性，造成图像中具有相同组织的灰度值随着位置的不同而发生改变，使图像在分割时受平滑变化的偏场影响比较严重，导致图像分割时难以取得令人满意的效果。

发明内容

本发明的目的就是为了解决上述问题，提供一种具有抗噪性和偏场校正的医学数字图像分割方法，它可用于处理带有噪声和偏场的医学图像的分割。

为了实现上述目的，本发明采用如下技术方案：

一种具有抗噪性和偏场校正的医学数字图像分割方法，为了能够更好的对带有偏场的图像进行偏场校正和分割，方法以构造的多项式拟合曲面作为聚类中心，通过迭代计算拟合曲面系数，隶属度和偏场完成对医学图像中偏场的纠正和图像的分割；另外，将邻域灰度差影响系数以及局部窗口的方差系数有效结合，设计一种新的影响因子对中心像素点进行合理的衡量，从而增强像素邻域信息的影响，起到削减抑制噪声的作用。

本发明的具体步骤为：

1、初始化隶属度矩阵和多项式拟合曲面的系数，并设置聚类个数、隶属度权重因子以及终止阈值；

2、通过求解目标函数最小化，计算更新偏场、拟合曲面系数以及隶属度；

2-1)计算更新偏场

2-2)计算更新拟合曲面系数a_k0、a_k1、a_k2；

2-3)计算更新隶属度

3、如果(U⁽ⁿ⁾-U⁽ⁿ⁺¹⁾)≥ε，则回到步骤2，否则算法迭代结束，跳到步骤4；

4、根据隶属度最大的准则分类标记每个像素点，得到图像的初步分割结果，记为I₀(x_ij)；

5、将图像损失能量函数E(x_ij)加入到图像的分割结果I₀(x_ij)中，完成最终图像分割。

本方法是一个迭代的过程，其中的三个变量是通过迭代运算相互交替不断更新计算的，每个变量的计算都是在另外两个变量的基础上求得的，在方法中，引入了残余偏场模型和抗噪模糊因子，可以使能量函数在迭代的同时去除偏场和噪声带来的影响，从而得到比较满意的图像分割结果。

本发明的有益效果：

(1)重新修改FCM算法的目标函数，将偏场估计模型加入该目标函数中，能够对MRI图像同时进行偏场校正和分割，能很大程度地消除偏场的影响，取得令人满意的分割效果。

(2)通过定义邻域像素与拟合曲面的空间距离，并将局部窗口的方差系数以及邻域窗口灰度差系数结合起来共同构造出新的抗噪模糊因子，由于加入了抗噪模糊因子并采用多项式拟合曲面代替聚类中心，能够更好的抵抗噪声影响和保持更精确的图像细节。

附图说明

图1是本发明的流程示意图；

具体实施方式

基于模糊聚类的分割算法主要思想是通过对聚类目标函数进行迭代，通过不断更新隶属度矩阵和聚类中心，使目标函数取得最小值，进而获取全局最优解。本发明方法以多项式拟合曲面作为聚类中心，通过迭代计算拟合曲面系数，隶属度以及偏场使目标函数达到最优，然后根据隶属度最大的准则分类标记像素点，并将能量损失函数加入初步分割后的图像中，完成带偏场的医学图像的偏场校正和最终分割。

下面结合附图与实施例对本发明作进一步说明。

如图1所示，包括以下步骤：

1、初始化隶属度矩阵和多项式拟合曲面的系数，并设置聚类个数、隶属度权重因子以及终止阈值；

用[0,1]间的随机数初始化隶属度矩阵U，设置聚类个数C，k＝1,2,...C、隶属度权重因子m以及终止阈值ε，初始化多项式拟合曲面f_k(x_i,y_i)＝a_k0x_i+a_k1y_i+a_k2的系数a_k0、a_k1、a_k2。

2、通过求解目标函数最小化，计算更新偏场、拟合曲面系数以及隶属度：

采用拟合曲面作为聚类中心的图像分割的目标函数定义为：

其中，b_i是当前像素点i的偏差场强度，a_k0、a_k1、a_k2是拟合曲面的系数，x_i、y_i是当前像素点i的横纵坐标，u_ki是当前像素点i相对于第k个曲面函数的隶属度，m是隶属度权重因子，P_i是当前像素点i的灰度值，N和C是像素总数和聚类数，G'_ki是影响因子，且：

其中，u_kj是当前像素点j相对于第k个曲面函数的隶属度，δ_ij是模糊因子，x_j、y_j是当前邻域窗口像素点的横纵坐标，P_j是落入以P_i像素为中心的邻域窗口中的像素。

2-1)计算更新偏场

通常，偏场是一组给定的基函数的线性组合g₁....g_M，其中基函数的数量M足够大，公式为：

目的是找到这组线性组合的一组最优系数w₁....w_M，并将该组最优系数表示为一个列向量w＝(w₁...w_M)^T，基函数g₁(i)....g_M(i)表示为一个列向量函数G(i)＝(g₁(i)...g_M(i))^T，通过公式(6)对偏场求偏导，得到偏场最优系数公式为：

其中，a_k0、a_k1、a_k2是拟合曲面的系数，x_i、y_i是当前像素点i的横纵坐标，u_ki是当前像素点i相对于第k个曲面函数的隶属度，P_i是当前像素点i的灰度值，G(i)是一组基函数。

通过上式的定义来估计偏场的值，公式为：

2-2)计算更新拟合曲面系数a_k0、a_k1、a_k2

根据公式(6)中目标函数对拟合曲面系数a_k0、a_k1、a_k2分别求偏导，得到更新拟合曲面系数的公式：

其中，

b_i是当前像素点i的偏差场强度，x_i、y_i是像素i的横纵坐标，u_ki是像素i相对于第k个曲面函数的隶属度，δ_ij是模糊因子，P_i是像素i的灰度值，x_j、y_j是像素i邻域像素的横纵坐标，P_j是以P_i像素为中心的邻域像素。

模糊因子δ_ij由以下公式给出：

一方面，通过邻域像素与拟合曲面的空间距离来定义邻域像素点的影响程度，空间距离因子δ_sd定义如下：

其中，d_ij是邻域像素与拟合曲面的空间距离。

另一方面，将局部窗口的方差系数以及邻域窗口灰度差系数结合定义一个模糊因子δ_sc。局部窗口的方差系数定义如下：

其中，Var(x)是局部窗口中灰度值方差，是局部窗口中灰度的平均值。

邻域灰度差系数反映邻域像素点和中心点的灰度值相似程度，定义如下：

S_ij＝||x_j-x_i||²,j∈N_i (14)

其中，j是邻域像素点，i是中心像素点，S_ij是第j个像素点与中心像素点之间的灰度差。

为了统一数据量纲，将局部方差系数C_u和邻域灰度差系数S_ij归一化到(0,1]之间，得：

可以得到新定义的模糊因子δ_sc：

其中，f_min是邻域窗口中局部方差数的最小值，f_max是最大值，同理，S_ij(min)是邻域窗口中局部灰度差系数的最小值，S_ij(max)是最大值，参数ζ的取值范围是(0,1]，其目的是为了保证κ_i和值不为零。

结合空间距离因子δ_sd和模糊因子δ_sc，最终放入目标函数的模糊因子δ_ij为：

δ_ij＝δ_sd×δ_sc (18)

2-3)计算更新隶属度

根据公式(6)中的目标函数对隶属度u_ki求偏导，得到更新隶属度的公式：

其中，a_k0、a_k1、a_k2是拟合曲面的系数，x_i、y_i是像素i的横纵坐标，P_i是像素i的灰度值，G'_ki是影响因子，由式子(7)定义。

3、如果(U⁽ⁿ⁾-U⁽ⁿ⁺¹⁾)≥ε，则返回到步骤2，否则算法迭代结束，跳到步骤4；

当连续两次迭代更新的隶属度值变化小于终止条件ε或方法超过最大迭代次数限制时迭代结束，跳到步骤4；否则跳到步骤2继续迭代更新；

4、根据隶属度最大的准则分类标记每个像素点，得到图像的初步分割结果，记为I₀(x_ij)；

5、将图像损失能量函数E(x_ij)加入到图像的初始分割结果I₀(x_ij)中，完成图像最终分割；

为了避免在去偏场过程中产生的图像有用信息丢失，采用奇异值去噪的方法对残余偏场进行处理。具体方法：设原始图像为I(x_ij)，偏场校正后的图像为I₀(x_ij)，图像中残余的偏场可近似表示为：

Bias(x_ij)＝I(x_ij)-I₀(x_ij) (20)

但这部分残余偏场中除了偏场通常还包含部分图像信息，因此对其采用奇异值去噪方法处理，得到偏场Bias⁰(x_ij)，则图像去偏场后的能量损失函数E(x_ij)可定义为：

E(x_ij)＝Bias(x_ij)-Bias⁰(x_ij) (21)

将图像损失能量函数E(x_ij)重新加入到图像的分割结果中，得到图像最终的分割公式为：

I₀'(x_ij)＝I₀(x_ij)+E(x_ij) (22)

综上，本发明是基于FCM算法的基础上，以多项式曲面作为聚类中心进行图像像素的聚类。一方面通过引入偏场模型校正图像偏场带来的影响，并对图像中存留的偏场残差进一步处理，起到了很好的偏场校正效果；另一方面，通过构造新的影响因子，对图像中存在的噪声起到了抑制作用，更加准确的保留了图像的原有信息，本发明的方法很大程度地消除了偏场带来的影响，能够取得令人满意的分割效果，同时加入了抗噪模糊因子并采用多项式拟合曲面代替聚类中心，能够更好的抵抗噪声影响和保持更精确的图像细节。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (1)

1.一种具有抗噪性和偏场校正的医学数字图像分割方法，其特征是，包含以下几个步骤：

步骤1，用[0,1]间的随机数初始化隶属度矩阵U，设置聚类个数C，k＝1,2,...C、隶属度权重因子m以及终止阈值ε，初始化多项式拟合曲面f_k(x_i,y_i)＝a_k0x_i+a_k1y_i+a_k2的系数a_k0，a_k1，a_k2，其中x_i,y_i是像素点i的横纵坐标；

步骤2，通过求解目标函数最小化，计算更新偏场、拟合曲面系数以及隶属度；

步骤3，当连续两次迭代更新的隶属度值变化小于终止条件ε或超过最大迭代次数时迭代结束，跳到步骤4；否则跳到步骤2-1)继续迭代更新；

步骤4，根据隶属度最大的准则分类标记每个像素点，得到图像的初步分割结果，记为I₀(x_ij)；

步骤5，将图像损失能量函数加入到图像的初步分割结果中，完成图像最终分割；

其中，所述步骤2包括：

定义图像分割的目标函数：

其中，b_i是像素点i的偏差场强度，a_k0,a_k1,a_k2是拟合曲面的系数，x_i,y_i是像素点i的横纵坐标，u_ki是像素点i相对于曲面k的隶属度，m是隶属度权重因子，P_i是像素点i的灰度值，N和C是像素总数和聚类数，G′_ki是影响因子，且：

其中，j是以像素点i为中心的邻域窗口中的像素点，u_kj是像素点j相对于曲面k的隶属度，δ_ij是模糊因子，x_i,y_i是像素点i的横纵坐标，P_j是像素点j的灰度值；

对目标函数求偏导，得到偏场、多项式系数和隶属度的更新公式分别为：

2-1)计算偏场公式为：

其中，G(i)是一组基函数，偏场最优系数的计算公式为：

其中，a_k0，a_k1，a_k2是拟合曲面的系数，x_i,y_i是像素点i的横纵坐标，u_ki是像素点i相对于曲面k的隶属度，P_i是像素点i的灰度值；

2-2)计算拟合曲面系数a_k0、a_k1、a_k2，公式为：

其中，

b_i是像素点i的偏差场强度，x_i,y_i是像素点i的横纵坐标，u_ki是像素点i相对于曲面k的隶属度，δ_ij是模糊因子，P_i是像素点i的灰度值，x_j,y_j是像素点j的横纵坐标，j是以像素点i为中心的邻域像素点；

模糊因子δ_ij由以下公式给出：

δ_ij＝δ_sd×δ_sc (18)

其中，空间距离因子δ_sd定义如下：

模糊因子δ_sc定义为：

其中，d_ij是邻域像素与拟合曲面的空间距离，k_i和φ_ij是归一化后的局部方差系数和邻域灰度差系数；

2-3)计算隶属度公式为：

其中，a_k0,a_k1,a_k2和a_l0,a_l1,a_l2分别是拟合曲面k和l的系数，x_i,y_i是像素点i的横纵坐标，P_i是像素点i的灰度值，G′_ki和G′_li分别是像素点i与曲面k和l的影响因子，由式子(7)定义；所述步骤5包括：

图像去偏场后的能量损失函数E(x_ij)定义为：

E(x_ij)＝Bias(x_ij)-Bias⁰(x_ij) (21)

其中，图像中残余的偏场Bias(x_ij)可表示为原始图像I(x_ij)和步骤4得到的图像I₀(x_ij)的差，偏场Bias⁰(x_ij)由对Bias(x_ij)采用奇异值去噪方法处理得到；

将图像损失能量函数E(x_ij)重新加入到步骤4图像的初步分割结果中，得到最终的分割图像I′₀(x_ij)，公式为：

I′₀(x_ij)＝I₀(x_ij)+E(x_ij) (22)。