CN101968851B - 基于字典学习上采样的医学影像处理方法 - Google Patents

Abstract

本发明公开了一种基于字典学习上采样的医学影像处理方法，它属于图像处理技术领域。其实现过程是：输入原始医学影像；对原始医学影像进行切割与增强处理；对切割与增强处理后的医学影像提取特征；根据提取的特征得到训练样本集和测试样本集；从训练样本中找出弱势样本的边界点；根据边界点邻域的情况得到需要生成新样本的个数；利用稀疏连线取点方法生成所需的新样本；将新样本增加进训练样本集中组成新的训练样本集；对新的训练样本集进行分类诊断，得到分类器；采用分类器，对测试样本集进行诊断，得到最终诊断结果。本发明具有对医学影像诊断识别率高泛化能力强的优点，可用于医务工作者评价疾病预后和治疗效果。

Description

基于字典学习上采样的医学影像处理方法

技术领域

本发明属于图像处理技术领域，特别是涉及医学影像处理，可用于监测疾病分布、研究发病机理及疾病辅助诊断。

背景技术

计算机科学技术的迅猛发展对医学领域产生了巨大的影响，人们试图让计算机逐渐代替人类实现自动诊断疾病等极具挑战性的工作，医学图像在临床诊断中发挥着重要的作用。自1895年伦琴发现X射线以来特别是1979年出现了计算机断层扫描CT技术，极大地推动了影像医学的发展。近三十年以来，新的影像学技术更是层出不穷。

然而，由于来源所需价值高以及个人隐私等问题，使得有病变的医学图像比正常图像较少，引起了数据不平衡问题，最终导致了医学图像识别困难。怎样有效地解决这种困难是医学图像领域中亟待解决的问题。目前，通用的方法就是利用上采样的方法，增加有病变医学图像的数据个数，改变分布，以降低数据的不平衡程度。

最原始的上采样方法是复制稀有类的样本，但是这样做对相同的数据重复学习了，耗费了时间，而且对提高稀有类识别率没有太大帮助。

较高级的上采样方法则采用一些启发式技巧，有选择地复制弱势样本，或者生成新的弱势样本。Chawla等人提出的SMOTE算法是一种简单有效的上采样方法，该方法首先为每个弱势样本随机选出几个邻近样本，并且在该样本与这些邻近的样本的连线上随机取点，生成无重复的新的弱势样本。2008年，Haibo He等人提出的ADASYN方法是对边界上的弱势样本按照自身的权重值利用SMOTE的连线取点的方法随机生成新的弱势样本。然而这些方法过多的学习了弱势样本，容易引起对弱势样本的过学习，使本身泛化能力弱，造成测试样本的总体识别率低。

发明内容

本发明的目的在于克服上述已有技术的不足，针对医学影像数据不平衡、诊断识别率低和泛化能力弱等问题，提出了一种基于字典学习上采样的医学影像处理方法，以提高医学影像诊断识别率和增强泛化能力。

为实现上述目的，本发明提供基于字典学习上采样的医学影像处理方法，包括如下步骤：

(1)输入原始医学影像；

(2)采用直方图均衡化和均方差标准化方法，对原始医学影像进行切割与增强处理；

(3)对切割与增强处理后的医学影像提取其灰度共生矩阵或Hu矩特征；定义有标识的训练样本集合T₁和未标识的测试样本集合T₂：

T₁＝{(x₁，u₁)..(x_i，u_i).，(x_l，u_l)}

T₂＝{v₁，..b_i，..v_m}

其中，x_i表示第i个有标识训练样本的特征，u_i表示第i个有标识训练样本的标识，V_i表示第i个未标识的测试样本的特征，l为有标识训练样本的个数；m为未标识的测试样本的个数；

(4)对有标识的训练样本T₁进行基于字典学习上采样的分类诊断，得到分类器C：

4a)从有标识的训练样本T₁中选取出样本个数最少的一类作为弱势样本T₃；

4b)找出弱势样本T₃中的边界点集合B＝{b₁，..b_i，..b_s}，计算每个边界点b_i需要生成的新样本的个数n_i，其中b_i表示第i个边界点，s为边界点的数量；

4c)对选取出的弱势样本T₃采用KSVD的方法训练生成一个字典D；

4d)将每个边界点b_i，通过稀疏连线取点方法得到n_i个新样本，用这些新样本组合成新样本集合T₄；

4e)将新样本集合T₄加入弱势样本T₃中组成当前新的弱势样本T₃′；

4f)将新的弱势样本T₃′与训练集中其他样本共同组成当前新的有标识的训练样本T₁′；

4g)用支持向量机方法对新的有标识的训练样本T₁′进行分类诊断，得到诊断后的分类器C；

(5)采用诊断后得到的分类器C，对未标识测试样本T₂进行诊断，得到对未标识测试样本T₂的最终诊断结果。

本发明与现有的技术相比具有以下优点：

1、本发明采用ADASYN算法中的边界点选取和SMOTE算法中的连线取点策略，结合二者的长处，提高了弱势样本的识别率；

2、本发明利用稀疏表示方法产生边界样本的虚拟点作为新增样本，克服了ADASYN算法和SMOTE算法对弱势样本的过学习问题，提高了总体识别率；

3、本发明采用所有的弱势样本进行字典学习，充分考虑了弱势样本的整体性能，增强了对医学影像处理的泛化能力。

附图说明

图1是本发明的流程图；

图2是本发明仿真采用的原始乳腺X影像图；

图3是本发明对图2进行去冗余后的乳腺X影像图；

图4是本发明对图3采用直方图均衡化和均方差标准化方法后的乳腺X影像图。

具体实施方式

参照图1，本发明基于字典学习上采样的医学影像处理方法，包括如下步骤：

步骤1：采用直方图均衡化和均方差标准化方法，对原始医学影像集中的医学影像进行切割与增强处理，得到视觉效果较好的医学影像集。

1a)输入原始医学影像，其大小为M×N，本实例选择如图2所示的原始乳腺X影像集中的一幅影像，其大小1024×1024；

1b)对输入的原始医学影像，采用计算机自动切割方法，切除其图像的背景和图像中存在的人为印记，得到切割后的乳腺X影像，如图3所示；

1c)对切割后的医学影像采用直方图均衡化和均方差标准化方法去除噪声，得到具有较好视觉效果的乳腺X影像，如图4所示。

步骤2：对切割与增强处理后的乳腺X影像进行Hu矩特征提取。

2a)在得到的视觉效果较好的乳腺X影像上，计算点(x，y)处的(p+q)阶矩m_pq及(p+q)阶中心矩μ_pq：

$m_{\mathrm{pq}}=\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{x}^p{y}^{p}f(x,y)$

${\mathrm{\μ}}_{\mathrm{pq}}=\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{(x-x_c)}^p{\left({y-y}_c\right)}^{q}f(x,y)$

式中，x表示点(x，y)的横坐标值，y表示点(x，y)的横坐标值；

M表示乳腺X影像的行高，N表示乳腺X影像的列高；

p∈{0，1，2，3}，q∈{0，1，2，3}，x^p表示x的p次方，y^q表示y的q次方；

f(x，y)表示点(x，y)处的像素值；

(x_c，y_c)表示视觉效果较好的乳腺X影像的重心坐标，x_c为横坐标值，y_c为纵坐标值；

(x-x_c)^p为(x-x_c)的p次方，(y-y_c)^q为(y-y_c)的q次方；

2b)根据得到的点(x，y)处的(p+q)阶矩m_pq及(p+q)阶中心矩μ_pq，得到点(x，y)的归一化中心矩的计算公式：

${\mathrm{\η}}_{\mathrm{pq}}={\mathrm{\μ}}_{\mathrm{pq}}/m_{\mathrm{pq}}^{\mathrm{\γ}}$

式中，

表示m_pq的γ次方，γ＝(p+q)/2+1；

2c)将2a)中p与q的值代入点(x，y)的归一化中心矩的计算公式中，得到点(x，y)的归一化中心矩η₀₂，η₀₃，η₁₁，η₁₂，η₂₀，η₂₁，η₃₀；

2d)利用点(x，y)处的归一化中心矩，提取视觉效果较好的乳腺X影像的七个Hu矩特征，分别定义为φ₁，φ₂，…，φ₇，即：

φ₁＝η₂₀+η₀₂

${\mathrm{\φ}}_2={({\mathrm{\η}}_{20}-{\mathrm{\η}}_{02})}^2+{4\mathrm{\η}}_{11}^2$

φ₃＝(η₃₀-3η₁₂)²+(3η₂₁-η₀₃)²

φ₄＝(η₃₀+η₁₂)²+(η₂₁+η₀₃)²

φ₅＝(η₃₀-3η₁₂)(η₃₀+η₁₂)φ_x+(η₀₃-3η₂₁)(η₂₁+η₀₃)φ_y

φ₆＝(η₂₀-η₀₂)[(η₃₀+η₁₂)²-(η₀₃+η₂₁)²]+4η₁₁(η₃₀+η₁₂)(η₀₃+η₂₁)

φ₇＝(3η₂₁-η₀₃)(η₃₀+η₁₂)φ_x+(η₃₀-3η₁₂)(η₀₃+η₂₁)φ_y，

其中，

φ_x＝(η₃₀+η₁₂)²-3(η₀₃-3η₂₁)²

φ_y＝(η₀₃+η₂₁)²-3(η₃₀-η₁₂)²

步骤3：根据对视觉效果较好的乳腺X影像提取的特征，得到有标识的训练样本集T₁和未标识的测试样本集T₂：

T₁＝{(x_i，u_i)|x_i∈Rⁿ，u_i∈{-1，1}，i＝1，…，l}

T₂＝{v₁，..v_i，..v_m}

式中，x_i表示第i个有标识训练样本的特征，u_i表示第i个有标识训练样本的标识，v_i表示第i个未标识的测试样本的特征，l为有标识训练样本的个数，m为未标识的测试样本的个数，本发明采用的乳腺X影像的诊断结果包括正常和癌症患者两类，正常的标记为-1，癌症患者标记为1。

步骤4：从有标识的训练样本T₁中选取出样本个数相对较少的一类作为弱势样本集T₁₁，另外一类做为强势样本集T₁₂：

T₁₁＝{(x_i，u_i)|x_i∈Rⁿ，u_i∈{1，2}，i＝1，…l₁}

T₁₂＝{(x_i，u_i)|x_i∈Rⁿ，u_i∈{1，2}，i＝1，…l₂}

式中，l₁为弱势样本集T₁₁的样本个数，l₂为强势样本集T₁₂的样本个数，Rⁿ为n维的实数集合。

步骤5：找出弱势样本T₁₁中的边界点集合B，计算每个边界点B(i)需要生成的新样本的个数n_i，其中i为整数，且1≤i≤s，s为边界点的数量，B和n_i按如下步骤计算：

5.1)计算弱势样本集T₁₁中所有样本的权重集合w₁₁ w _L1 ：

w_L1＝{w_i|i＝1，…l₁}  w _L1 ＝{w _i |i＝1，…l ₁ }：

式中，w_i＝num_i/K，表示样本(x_i，u_i)对应的权重，其中，K为样本x_i的邻域个数，K＝5；num_i代表样本x_i邻域内强势样本的个数；

5.2)根据权重集合w₁₁从弱势样本集T₁₁中选择出边界点集合B以及它的权重集合w_B：

B＝{(x_i，u_i)|0＜w_i＜1，w_i∈w₁₁}

w_B＝{w_i|0＜w_i＜1，w_i∈w₁₁}

其中，(x_i，u_i)为有标识的训练样本集T₁中的第i个样本，w_i为权重集合w₁₁中的第i个权重；

5.3)将权重集合w_B归一化，得到w_B′：

${w_B}^{\′}=\{w_i/\underset{i=1}{\overset{l_3}{\mathrm{\Σ}}}{w}_i|i=1,\·\·\·l_3,w_i\∈w_B\text{}}$

式中，l₃为边界点集合B的个数，w_B为边界点集合B归一化后的权重集合；

5.4)计算每个边界点B(i)需要生成的新样本的个数n_i：

n_i＝w_B′(i)*(l₂-l₁)

其中，w_B′(i)为w_B′中的第i个元素，l₁为弱势样本集T₁₁的样本个数，l₂为强势样本集T₁₂的样本个数。

步骤6：对选取出的弱势样本T₃采用KSVD的方法训练生成一个字典D。

步骤7：将每个边界点B(i)，通过稀疏连线取点方法得到n_i个新样本，每个新样本的生成步骤如下：

7.1)利用步骤6得到的字典D对B(i)稀疏表示得到稀疏系数α₁；

7.2)随机选取B(i)近邻中的一个弱势样本，利用步骤6得到的字典D对该样本稀疏表示得到稀疏系数α₂；

7.3)通过α₁和α₂计算所要生成新样本的稀疏系数α：

α＝α₁+(α₂-α₁)*rand

其中，rand为0到1的一个随机数；

7.4)将字典D与稀疏系数α相乘，得到所要生成的一个新样本T_new。

步骤8：将步骤7生成的所有新样本加入弱势样本T₁₁中组成当前新的弱势样本T₁₁′。

步骤9：将新的弱势样本T₁₁′与强势样本集T₁₂共同组成当前新的有标识的训练样本T₁′。

步骤10：用支持向量机方法对训练样本T₁′进行分类诊断，得到诊断后的分类器。

步骤11：采用诊断后得到的分类器，对未标识测试样本T₂进行诊断，得到类别标签h₁，h₂，...，h_l(h_i∈{-1，1})，若类别标签为1，诊断结果为癌变；若类别标签为-1，诊断结果为正常。

步骤12：输出未标识测试样本T₂的最终诊断结果。

本发明的效果可以通过以下对乳腺X影像仿真数据进一步说明：

1.实验条件和内容

实验仿真环境为：windows XP，SPI，CPU Pentium(R)4，基本频率2.4GHZ，软件平台为Matlab7.0.4运行。仿真选用的原始乳腺X影像来源于公共数据集MIAS，共获取322幅原始乳腺X影像。

实验内容包括：分别采用复制样本的上采样方法COPY、SMOTE、ADASYN以及本发明四种方法对原始乳腺X影像进行试实验，检测患有癌症的图像。实验结果的验证中，我们列出了四种评价指标来证明本发明的优越性，分别为样本总识别率T_Ratio，小类样本识别率Recall，总的评价指标F_measure和G_mean。

为了定义这些评价指标，需要用到机器学习方法的基础评价指标——混淆矩阵，如表1所示。其中，用正类表示弱势样本，负类表示强势样本。

表1两类混淆矩阵

	预测正类	预测反类
			实际正类	TP(true positive)	FN(false negative)；
实际反类	FP(false positive)；	TN(true negative)

在一个两类混淆矩阵中，实际为正类，预测也为正类的样本数量称为正确正类TP；实际为正类，预测为反类的称为错误反类FN；实际为反类，预测为正类的称为错误正类FP；实际为反类，预测为反类的称为正确反类TN。

利用混淆矩阵定义的四种评价指标如表2所示：

表2四种评价指标的定义

2、仿真结果

在仿真试验中，支持向量机里面的核映射方法选用的是高斯核，核参数取的是0.005。本发明的字典优化所用的是KSVD的方法。KSVD的原子个数取10，迭代次数也是10，稀疏表示误差是0.01。实验保证各种方法在同样的数据下运行。

表3为四种方法对322幅图像整体的诊断结果。其中，有癌变的图像是115副，没有癌变的正常样本是207副，不平衡度为1∶1.8。实验结果为20次实验的平均结果，每次实验都是随机采取70％的样本作为有标记的训练样本，其余的30％作为无标记的测试样本。

表3四种方法对322幅图像整体的诊断结果

	Ratio	Recall	F_measure	G_mean
					COPY	0.48434	0.48727	0.46739	0.52909
SMOTE	0.49219	0.59545	0.52316	0.55133
					ADASYN	0.47656	0.51000	0.47478	0.52506
D_SMOTE	0.52344	0.66636	0.56835	0.59570

从表3可以看到，本发明在所有的评价指标中都是最优的，这充分说明本发明在处理医学图像不平衡问题上的优越性。

将处理好的医学图像平行切成4份，将得到1288副图像。根据癌变的位置，对这1288副图像标记类别，1代表有癌变的，-1代表正常的。随机去掉420副正常图像，最终得到有癌变的图像120副，正常图像748副，共968副，不平衡度为1∶6.23。四种方法对968幅图像整体的实验结果如表4所示。实验结果为20次实验的平均结果，每次实验都是随机采取50％的样本作为有标记的训练样本，其余的50％作为无标记的测试样本。

表4四种方法对968幅图像整体的实验结果

	Ratio	Recall	F_measure	G_mean
					COPY	0.82855	0.49764	0.49688	0.65083
SMOTE	0.81235	0.59516	0.50001	0.70453
					ADASYN	0.80210	0.60386	0.49782	0.70336
D_SMOTE	0.90012	0.56518	0.62115	0.73446

从表4可以看到，本发明在不平衡度比较高的情况下，对弱势样本的识别率相对较差。这是因为本发明防止对弱势样本足够多的学***衡度较大时，只有对弱势样本足够大的学习与训练才能把弱势样本识别率提高上去，但是这样影响了总的识别率与分类器的泛化能力。本发明在总识别率与泛化能力上就有明显的优势。

Claims (1)

1.一种基于字典学习上采样的医学影像处理方法，包括如下步骤：

(1)输入原始医学影像；

(2)采用直方图均衡化和均方差标准化方法，对原始医学影像进行切割与增强处理；

(3)对切割与增强处理后的医学影像提取其灰度共生矩阵或Hu矩特征；定义有标识的训练样本集合T₁和未标识的测试样本集合T₂：

T₁＝{(x₁，u₁)..(x_i，u_i).，(x_l，u_l)}

T₂＝{v₁，..v_i，..v_m}

其中，x_i表示第i个有标识训练样本的特征，u_i表示第i个有标识训练样本的标识，v_i表示第i个未标识的测试样本的特征，l为有标识训练样本的个数，m为未标识的测试样本的个数；

(4)对有标识的训练样本T₁进行基于字典学习上采样的分类诊断，得到分类器C：

4a)从有标识的训练样本T₁中选取出样本个数最少的一类作为弱势样本T₃；

4b)找出弱势样本T₃中的边界点集合B＝{b₁，..b_i，..b_s}，计算每个边界点b_i需要生成的新样本的个数n_i，其中b_i表示第i个边界点，s为边界点的数量：

4b1)计算弱势样本集T₃中所有样本的权重集合w_L1：

w_L1＝{w_i|i＝1，…l₁}

式中，l₁为弱势样本集T₃的个数，w_i＝num_i/K，表示样本(x_i，u_i)对应的权重，

其中，K为样本(x_i，u_i)的邻域个数，K＝5；num_i代表样本(x_i，u_i)邻域内强势样本的个数；

4b2)根据权重集合w_L1从弱势样本集T₃中选择出边界点集合B以及边界点的权重集合w_B：

$w_B=\{w_b^1,..w_b^i,..w_b^s\}$

B＝{b₁，..b_i，..b_s}

其中，

表示b_i的权重；

4b3)将权重集合w_B归一化，得到w_B′：

${w_B}^{\′}=\{w_b^i/\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{w}_b^i|i=1,\·\·\·s\}$

其中，w_B′为边界点集合B归一化后的权重集合；

4b4)计算每个边界点b_i需要生成的新样本的个数：

n_i＝w_B′(i)*(l₂-l₁)

其中，w_B′(i)为w_B′中的第i个元素，l₂为训练样本集合T₁中除弱势样本集T₃之外的其他样本的个数；

4c)对选取出的弱势样本T₃采用KSVD的方法训练生成一个字典D；

4d)将每个边界点b_i，通过稀疏连线取点方法得到n_i个新样本，用这些新样本组合成新样本集合T₄：

4d1)利用步骤4c)得到的字典D对b_i稀疏表示得到稀疏系数α₁；

4d2)随机选取b_i近邻中的一个弱势样本，利用步骤4c)得到的字典D对该样本稀疏表示得到稀疏系数α₂；

4d3)通过α₁和α₂计算所要生成新样本的稀疏系数α：

α＝α₁+(α₂-α₁)*rand

其中，rand为0到1的一个随机数；

4d4)将字典D与稀疏系数α相乘，得到所要生成的一个新样本T_new；

4e)将新样本集合T₄加入弱势样本T₃中组成当前新的弱势样本T₃′；

4f)将新的弱势样本T₃′与训练样本集合T₁中除弱势样本集T₃之外的其他样本共同组成当前新的有标识的训练样本T₁′；

4g)用支持向量机方法对新的有标识的训练样本T₁′进行分类诊断，得到诊断后的分类器C；

(5)采用诊断后得到的分类器C，对未标识测试样本T₂进行诊断，得到对未标识测试样本T₂的最终诊断结果。

Images