CN103854262B - 基于结构聚类与稀疏字典学习的医学图像降噪方法 - Google Patents

Abstract

本发明公开了一种基于结构聚类与稀疏字典学习的医学图像降噪方法，克服了现有技术中字典学习方法忽略了图像块的结构相似性，导致损失图像的细节信息，图像的均匀区域不光滑的问题。本发明实现步骤为：(1)输入噪声图像；(2)预滤波；(3)结构聚类；(4)提取各类训练样本集合；(5)分类学习稀疏字典；(6)得到最终的训练字典；(7)得到各类去噪估计值；(8)得到噪声抑制后的图像；(9)输出噪声抑制后的图像。本发明具有有效去除噪声，提高图像的视觉效果，边缘、纹理得到有效的保持，均匀区域光滑的优点，可用于对医学计算机断层扫描、核磁共振图像的降噪处理。

Description

基于结构聚类与稀疏字典学习的医学图像降噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及含有加性噪声的图像降噪技术领域中的一种基于结构聚类与稀疏字典学习的医学图像降噪方法。本发明可用于对医学计算机断层扫描（computedtomographyCT）、核磁共振（magneticresonanceMR）图像进行降噪处理。

背景技术

图像降噪的目的是在去除噪声的同时能够很好地保护图像的特征信息，例如图像的边缘和纹理。在图像的采集、处理和传输的过程中，由于技术的限制和设备本身的固有特性，使得图像不可避免地含有各种噪声，导致图像质量严重下降，极大地影响后续的图像处理。因此，图像降噪具有非常重要的意义。

王昊，康晓东，刘玲玲，耿佳佳发表的论文“基于非下采样Contourlet变换的医学CT图像去噪”（《计算机工程与应用》，2012，48（27），150-154）中提出了一种基于非下采样Contourlet变换的医学CT图像去噪方法。该方法的具体步骤是：首先，对含噪的CT图像进行非下采样Contourlet变换，分解成低频系数矩阵和高频系数矩阵，再对高频系数矩阵利用Context模型处理系数子矩阵并进行阈值去噪处理，最后对去噪后的变换系数进行非下采样Contourlet逆变换得到去噪后的CT图像。该方法虽然弥补了Contourlet变换在平移不变性方面的不足，降低了去噪图像中的伪吉布斯现象，但仍然存在的不足是，该方法恢复得到的图像不够清晰，损失了部分细节信息，增加了后续图像解译的难度。

西安电子科技大学在其专利申请“基于稀疏自适应字典的图像去噪方法”（专利申请号：201310161279.7，授权公告号：CN103218791A）中提出了一种基于稀疏自适应字典的图像去噪方法。该专利技术利用从含噪图像中抽取得到的训练数据集迭代地训练字典，在迭代的过程中，将上次迭代所得的字典作为本次迭代的基字典，利用迭代完成后获得的最终字典和编码系数矩阵对噪声图像进行去噪。该方法虽然能够增强字典的自适应性，但仍然存在的不足是，字典学习方法忽略了图像块的结构相似性，不能有效的保持图像的边缘，对图像的均匀区域去噪效果不够理想。

发明内容

本发明的目的在于克服上述现有技术的不足，提出了一种基于结构聚类与稀疏字典学习的医学图像降噪方法，增强了字典的自适应性和结构性，在有效去除噪声的同时，能够很好地保持图像的边缘和纹理等细节信息。

为实现上述目的，本发明的具体步骤如下：

(1)输入噪声图像：

输入任选的一幅含有加性高斯白噪声的噪声图像。

(2)预滤波：

采用广义K均值K-SVD方法，对噪声图像进行滤波处理，得到滤波后图像。

(3)结构聚类：

(3a)设定滤波后图像中以任一像素点为中心的邻域窗口大小为9×9，利用控制核回归权值公式，求出任一像素点与该像素点的邻域窗口内其他像素点的控制核回归权值向量，取控制核回归权值向量的均值向量作为当前邻域窗口的特征向量；

(3b)采用K均值聚类的方法，对滤波后图像中所有的邻域窗口的特征向量进行聚类，得到滤波后图像中各个邻域窗口的类别。

(4)提取各类训练样本集合：

(4a)从步骤（1）输入的噪声图像中提取大小为9×9的噪声图像块；

(4b)在滤波后图像中，找出与所提取的噪声图像块索引位置相同的邻域窗口，将该邻域窗口的类别作为该噪声图像块的类别；

(4c)合并噪声图像中类别相同的噪声图像块，得到6类噪声图像块集合，将6类噪声图像块集合中的每一个图像块向量化，得到6类图像块向量集合，将图像块向量集合作为训练样本集合。

(5)分类学习稀疏字典：

采用稀疏字典学习方法，利用6类训练样本集合进行稀疏字典学习，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数。

(6)得到最终的训练字典：

按照下式，得到每一类训练样本集合最终的训练字典：

$\mathrm{Dic}=\mathrm{\Φ}\·\tilde{S}$

其中，Dic表示每一类训练样本集合最终的训练字典，Φ表示基础字典，表示每一类训练样本集合最终更新后的稀疏字典。

(7)得到各类去噪估计值：

按照下式，得到每一类训练样本集合的去噪估计值：

Y＝Dic·Γ

其中，Y表示每一类训练样本集合的去噪估计值，Dic表示每一类训练样本集合最终的训练字典，Γ表示每一类训练样本集合最终更新后的稀疏编码系数。

(8)得到噪声抑制后的图像：

(8a)对步骤（1）输入的噪声图像中的任意一个像素点位置，从6类训练样本集合的去噪估计值中提取出与该像素点位置对应的所有像素值，将所有像素值的均值作为该像素点的去噪结果；

(8b)对步骤（1）输入的噪声图像中的所有位置的像素点，按照步骤(8a)的步骤，得到所有像素点的去噪结果，获得最终噪声抑制后的图像。

(9)输出噪声抑制后的图像。

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

第一，由于本发明利用6类训练样本集合进行稀疏字典学习，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数，克服了现有技术中恢复得到的图像不够清晰，损失部分细节信息的缺点，使得本发明具有有效去除图像的噪声，恢复后的图像视觉效果好的优点。

第二，由于本发明利用控制核回归权值向量的均值向量，作为邻域窗口的特征向量，采用K均值聚类的方法，对邻域窗口的特征向量进行聚类，克服了现有技术中字典学习方法忽略了图像块的结构相似性的缺点，使得本发明具有有效捕捉图像的结构特征，恢复后的图像边缘保持清晰，均匀区域光滑的优点。

附图说明

图1是本发明的流程图；

图2是本发明与现有技术对计算机断层扫描图像降噪效果对比图；

图3是本发明与现有技术对核磁共振图像降噪效果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1，本发明的具体步骤如下：

步骤1：输入噪声图像。

输入任选的一幅含有加性高斯白噪声的噪声图像。

步骤2：预滤波。

采用广义K均值K-SVD方法，对噪声图像进行滤波处理，得到滤波后图像。

所述的广义K均值K-SVD方法的具体步骤如下：

第一步，设置图像块的大小为8×8，超完备字典的大小为64×256，初始化超完备字典为离散余弦变换DCT字典。

第二步，按照下式，求出稀疏编码系数：

${\mathrm{\α}}_{\mathrm{mn}}=\arg \min \mathrm{\μ}{\left|\right|{\mathrm{\α}}_{\mathrm{mn}}\left|\right|}_0+{\left|\right|{\mathrm{D\α}}_{\mathrm{mn}}-R_{\mathrm{mn}}X\left|\right|}_2^2$

其中，α_mn表示稀疏编码系数，m，n表示任意图像块在噪声图像中的坐标位置，argmin表示求最小值操作，μ表示残差控制因子，||·||₀表示取零范数操作，表示取二范数的平方操作，D表示超完备字典，R_mn表示图像块提取操作符，X表示去噪后的图像。

第三步，按照下式，计算误差矩阵：

$E_k=\{R_{\mathrm{mn}}X-\underset{t\≠k}{\mathrm{\Σ}}{d}_t{\mathrm{\α}}_{\mathrm{mn}}\left(t\right)|(m,n)\∈w_k\}$

其中，E_k表示误差矩阵，k表示超完备字典的第k个原子的编号，R_mn是图像块提取操作符，m，n表示任意图像块在噪声图像中的坐标位置，X表示去噪后的图像，d_t表示超完备字典的第t个原子，α_mn表示稀疏编码系数，α_mn(t)表示稀疏编码系数α_mn的第t个元素，w_k表示使用了超完备字典的第k个原子d_k的所有图像块的坐标位置构成的集合。

第四步，按照下式，对误差矩阵进行奇异值分解：

E_k＝UΔV^T

其中，E_k表示误差矩阵，U表示左奇异矩阵，Δ表示奇异值矩阵，V表示右奇异矩阵，V^T表示右奇异矩阵V的转置。

第五步，用左奇异矩阵U的第一列更新超完备字典的第k列原子d_k，用右奇异矩阵V的第一列乘以残差矩阵E_k的最大特征值Δ(1,1)更新集合w_k对应的图像块的稀疏编码系数α_mn的第k个元素α_mn(k)。

第六步，对超完备字典中所有原子分别重复执行步骤2的第三步至第五步的操作，得到更新后的超完备字典和稀疏编码系数。

第七步，重复执行步骤2的第二步至第六步5次，得到最终更新后的超完备字典和最终更新后的稀疏编码系数。

第八步，按照下式，得到滤波后图像：

$Z^{\′}={(\mathrm{\λE}+\underset{\mathrm{mn}}{\mathrm{\Σ}}{R}_{\mathrm{mn}}^T{R}_{\mathrm{mn}})}^{-1}(\mathrm{\λZ}+\underset{\mathrm{mn}}{\mathrm{\Σ}}{R}_{\mathrm{mn}}^T\hat{D}{\widehat{\mathrm{\α}}}_{\mathrm{mn}})$

其中，Z′表示滤波后图像，λ表示拉格朗日系数，E表示单位矩阵，R_mn表示重叠块提取操作符，m，n表示任意图像块在噪声图像中的坐标位置，表示重叠块提取操作符R_mn的转置，表示对进行求逆运算，Z表示噪声图像，表示最终更新后的超完备字典，表示最终更新后的稀疏编码系数。

步骤3：结构聚类。

设定滤波后图像中以任一像素点为中心的邻域窗口大小为9×9，利用控制核回归权值公式，求出任一像素点与该像素点的邻域窗口内其他像素点的控制核回归权值向量，取控制核回归权值向量的均值向量作为当前邻域窗口的特征向量。

控制核回归权值公式如下：

${\mathrm{\ω}}_{a,b}=\frac{\sqrt{\det \left(M_b\right)}}{2\mathrm{\π}{h}^2}\exp \{-\frac{{(p_a-p_b)}^T{M}_b(p_a-p_b)}{2h^2}\}$

其中，ω_a,b表示控制核回归权值，a表示滤波后图像中任意的像素点，b表示以像素点a为中心的邻域窗口内其他的像素点，表示开方操作，det(·)表示求方阵的行列式操作，M_b表示由像素点b在水平和垂直方向上的估计梯度构成的对称梯度协方差矩阵，h表示平滑参数，exp(·)表示取指数操作，p_a和p_b分别表示像素点a和像素点b在滤波后图像中的位置，(·)^T表示转置操作。

采用K均值聚类的方法，对滤波后图像中所有的邻域窗口的特征向量进行聚类，得到滤波后图像中各个邻域窗口的类别。

K均值聚类的方法的具体步骤如下：

第一步，从滤波后图像中任意选取6个邻域窗口的特征向量作为6类邻域窗口的聚类中心。

第二步，按照下式，计算滤波后图像的邻域窗口与6类聚类中心的特征向量的欧式距离：

${\mathrm{dis}}_{j,q}=\left|\right|{\stackrel{\‾}{\mathrm{\ω}}}_j-v_q\left|\right|$

其中，dis_j,q表示滤波后图像的邻域窗口与6类聚类中心的特征向量的欧式距离，j表示滤波后图像中邻域窗口的序号，q表示滤波后图像中邻域窗口的类别，q＝1,2,...,6，表示滤波后图像中第j个邻域窗口的特征向量，v_q表示滤波后图像中的第q类邻域窗口的聚类中心，||·||表示求欧式距离操作。

第三步，将滤波后图像中各个邻域窗口分配到与滤波后图像中6类聚类中心的特征向量的欧式距离最小的对应类别中。

第四步，对滤波后图像中类别相同的邻域窗口的特征向量求平均向量，作为滤波后图像中该类邻域窗口新的聚类中心。

第五步，将步骤3的第二步至第四步重复进行100次，得到滤波后图像中各个邻域窗口的类别。

步骤4：提取各类训练样本集合。

从步骤1输入的噪声图像中提取大小为9×9的噪声图像块。

在滤波后图像中，找出与所提取的噪声图像块索引位置相同的邻域窗口，将该邻域窗口的类别作为该噪声图像块的类别。

合并噪声图像中类别相同的噪声图像块，得到6类噪声图像块集合，将6类噪声图像块集合中的每一个图像块向量化，得到6类图像块向量集合，将图像块向量集合作为训练样本集合。

步骤5：分类学习稀疏字典。

采用稀疏字典学习方法，利用6类训练样本集合进行稀疏字典学习，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数。

稀疏字典学习方法的具体步骤如下：

第一步，分别设置一个大小为81×100的基础字典和一个大小为100×100的稀疏字典，令基础字典为离散余弦变换DCT字典，初始化稀疏字典为单位矩阵。

第二步，按照下式，计算稀疏编码系数：

$\begin{array}{ccc}\∀i& \min {\left|\right|{\mathrm{\γ}}_i\left|\right|}_0& s.t.{\left|\right|y_i-\mathrm{\ΦS}{\mathrm{\γ}}_i\left|\right|}_2^2\≤{\mathrm{C\σ}}^2\end{array}$

其中，表示取训练样本序号中的任意值，i表示训练样本的序号，min表示求最小值操作，||·||₀表示取零范数操作，γ_i表示稀疏编码系数，s.t.表示约束条件，表示取二范数的平方操作，y_i表示由训练样本的序号对应的训练样本，Φ表示基础字典，S表示稀疏字典，C表示高斯白噪声的增益，σ²表示高斯白噪声的方差。

第三步，按照下式，得到使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合：

I_l＝{i|γ_i(l)≠0,i∈{1,2,...,N}}

其中，I_l表示使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合，i表示训练样本的序号，γ_i表示稀疏编码系数，γ_i(l)表示稀疏编码系数γ_i的第l个元素，N表示每一类训练样本的总个数。

第四步，按照下式，计算残差矩阵：

$E_l=\{y_i-\underset{f\≠l}{\mathrm{\Σ}}{\mathrm{\Φs}}_f{\mathrm{\γ}}_i\left(f\right)|i\∈I_l\}$

其中，E_l表示残差矩阵，l表示稀疏字典的第l个原子的编号，i表示训练样本的序号，y_i表示由训练样本的序号对应的训练样本，∑表示求和操作，f表示稀疏字典的第f个原子的编号，Φ表示基础字典，s_f表示稀疏字典的第f个原子，γ_i表示稀疏编码系数，γ_i(f)表示稀疏编码系数γ_i的第f个元素，I_l表示使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合。

第五步，按照下式，计算残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影：

z＝E_l·g

其中，z表示残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影，E_l表示残差矩阵，g表示集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量。

第六步，按照下式，计算投影z的稀疏编码系数：

其中，表示投影z的稀疏编码系数，argmin表示求最小值操作，表示取二范数的平方操作，z表示残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影，Φ表示基础字典，s.t.表示约束条件，||·||₀表示取零范数操作，L表示稀疏字典原子的稀疏度，L＝24。

第七步，按照下式，对投影z的稀疏编码系数进行规则化：

其中，表示对投影z的稀疏编码系数进行规则化后的稀疏编码系数，||·||₂表示取二范数操作，Φ表示基础字典。

第八步，按照下式，对稀疏字典的第l个原子进行更新：

其中，表示对稀疏字典的第l个原子更新后的新原子，表示对投影z的稀疏编码系数进行规则化后的稀疏编码系数。

第九步，按照下式，对集合I_l对应的训练样本的稀疏编码系数进行更新：

${\widetilde{\mathrm{\γ}}}_i\left(l\right)={\left({\left(E_l\right)}^T\mathrm{\Φ}{\tilde{s}}_l\right)}^T$

其中，表示对集合I_l对应的训练样本的稀疏编码系数进行更新后的稀疏编码系数，i表示训练样本的序号，l表示稀疏字典的第l个原子的编号，E_l表示残差矩阵，(·)^T表示转置操作，Φ表示基础字典，表示对稀疏字典第l个原子更新后的新原子。

第十步，对稀疏字典中所有原子分别重复执行步骤5的第三步至第九步的操作，得到更新后的稀疏字典和稀疏编码系数。

第十一步，将步骤5的第二步至第十步重复执行15次，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数。

步骤6：得到最终的训练字典。

按照下式，得到每一类训练样本集合最终的训练字典：

$\mathrm{Dic}=\mathrm{\Φ}\·\tilde{S}$

其中，Dic表示每一类训练样本集合最终的训练字典，Φ表示基础字典，表示每一类训练样本集合最终更新后的稀疏字典。

步骤7：得到各类去噪估计值。

按照下式，得到每一类训练样本集合的去噪估计值：

Y＝Dic·Γ

其中，Y表示每一类训练样本集合的去噪估计值，Dic表示每一类训练样本集合最终的训练字典，Γ表示每一类训练样本集合最终更新后的稀疏编码系数。

步骤8：得到噪声抑制后的图像。

对步骤1输入的噪声图像中的任意一个像素点位置，从6类训练样本集合的去噪估计值中提取出与该像素点位置对应的所有像素值，将所有像素值的均值作为该像素点的去噪结果。

对步骤1输入的噪声图像中的所有位置的像素点，按照上面的步骤，得到所有像素点的去噪结果，获得最终噪声抑制后的图像。

步骤9：输出噪声抑制后的图像。

下面结合仿真图对本发明的效果做进一步说明。

1，仿真实验条件：

本发明的硬件测试平台是：处理器为InterCorei3-2370MCPU，主频为2.40GHz，内存4GB，软件平台为：Windows7旗舰版32位操作系统和MatlabR2012b。本发明输入的噪声图像分别为计算机断层扫描图像和核磁共振图像，其中计算机断层扫描图像大小为260×190，核磁共振图像大小为256×256，格式都为PNG。

2，仿真内容：

在仿真中用到的三个与本发明进行效果对比的现有技术所采用对比方法，分别如下：

MichaelElad等人在文献“ImageDenoisingViaSparseandRedundantRepresentationsOverLearnedDictionaries[J].IEEETransactionsonImageProcessing,2006,15(12):3736-3745.”中提出的利用冗余DCT字典方法和广义K均值K-SVD方法的2种图像去噪方法，分别简称ODCT方法和K-SVD方法。

RonRunbinstein等人在文献“Doublesparsity:learningsparsedictionariesforsparsesignalapproximation[J].IEEETransactionsonSignalProcessing,2010,58(3):1553-1564.”中提出的利用稀疏K-SVD方法的图像去噪方法，简称S-KSVD方法。

图2是本发明与现有技术对计算机断层扫描图像降噪效果对比图。其中，图2(a)是干净的计算机断层扫描图像，大小为260×190。图2(b)是对图2(a)添加噪声方差为25的高斯白噪声得到的噪声计算机断层扫描图像。图2(c)是现有技术中采用ODCT方法对图2(b)进行降噪的结果示意图。图2(d)是为现有技术中采用K-SVD方法对图2(b)进行降噪的结果示意图。图2(e)为现有技术中采用S-KSVD方法对图2(b)进行降噪的结果示意图。图2(f)为本发明方法对图2(b)进行降噪的结果示意图。

图3是本发明与现有技术对核磁共振图像降噪效果对比图，其中，图3(a)是干净的核磁共振图像，大小为256×256，图3(b)是对图3(a)添加噪声方差为10的高斯白噪声得到的噪声核磁共振图像，图3(c)是现有技术中采用ODCT方法对图3(b)进行降噪的结果示意图，图3(d)为现有技术中采用K-SVD方法对图3(b)进行降噪的结果示意图，图3(e)为现有技术中采用S-KSVD方法对图3(b)进行降噪的结果示意图，图3(f)为本发明方法对图3(b)进行降噪的结果示意图。

3.仿真结果分析：

图2和图3分别为本发明与现有技术对计算机断层扫描图像和核磁共振图像进行降噪的效果对比图。从图2(c)和图3(c)可以看出，现有技术中的ODCT方法获得的去噪结果示意图中，在图像的边缘处出现了很明显的振铃效应，导致图像边缘变模糊，细节不清晰，恢复得到的去噪图像质量低。从图2(d)和图3(d)可以看出，现有技术中的K-SVD方法获得的结果示意图中，噪声去除的较干净，图像边缘的振铃效应明显地降低，细节信息得到了较好的恢复，但是图像的均匀区域恢复效果不佳。从图2(e)和图3(e)可以看出，现有技术中的S-KSVD方法获得的结果示意图中，噪声得到了有效地抑制，但是在图像的边缘仍然存在比较明显的振铃效应，细节信息比较模糊。从图2(f)和图3(f)可以看出，本发明方法获得的去噪图中，不仅噪声去除得非常干净，图像的边缘和纹理信息恢复的较好，而且均匀区域最光滑。

对图2(a)中的干净计算机断层扫描图像分别添加噪声方差为10,20,25,30,40,50的高斯白噪声，获得含有不同程度的高斯白噪声的噪声计算机断层扫描图像，分别使用现有技术中的ODCT方法，K-SVD方法，S-KSVD方法和本发明方法对这些噪声计算机断层扫描图像进行去噪，计算去噪后的计算机断层扫描图像整体的峰值信噪比PSNR，结果如表1所示。

表1对含有不同噪声方差的计算机断层扫描图像进行降噪的峰值信噪比（PSNR）

从表2中可以看出，与其他三种现有技术的ODCT方法，K-SVD方法和S-KSVD方法相比，本发明方法获得的降噪核磁共振图像的PSNR值最高，表明本发明方法的去噪效果比现有技术的ODCT方法，K-SVD方法和S-KSVD方法好。

对图3(a)中的干净核磁共振图像分别添加噪声方差为10,20,25,30,40,50的高斯白噪声，获得含有不同程度的高斯白噪声的噪声核磁共振图像，分别使用现有技术中的ODCT方法，K-SVD方法，S-KSVD方法和本发明方法对这些噪声核磁共振图像进行去噪，计算去噪后的核磁共振图像整体的峰值信噪比PSNR，结果如表2所示。

表2对含有不同噪声方差的核磁共振图像进行降噪的峰值信噪比（PSNR）

从表2中可以看出，与其他三种现有技术的ODCT方法，K-SVD方法和S-KSVD方法相比，本发明方法获得的降噪核磁共振图像的PSNR值最高，表明本发明方法的去噪效果比现有技术的ODCT方法，K-SVD方法和S-KSVD方法好。

Claims (5)

1.一种基于结构聚类与稀疏字典学习的医学图像降噪方法，包括如下步骤：

(1)输入噪声图像：

输入任选的一幅含有加性高斯白噪声的噪声图像；

(2)预滤波：

采用广义K均值K-SVD方法，对噪声图像进行滤波处理，得到滤波后图像；

(3)结构聚类：

(3a)设定滤波后图像中以任一像素点为中心的邻域窗口大小为9×9，利用控制核回归权值公式，求出任一像素点与该像素点的邻域窗口内其他像素点的控制核回归权值向量，取控制核回归权值向量的均值向量作为当前邻域窗口的特征向量；

(3b)采用K均值聚类的方法，对滤波后图像中所有的邻域窗口的特征向量进行聚类，得到滤波后图像中各个邻域窗口的类别；

(4)提取各类训练样本集合：

(4a)从步骤(1)输入的噪声图像中提取大小为9×9的噪声图像块；

(4b)在滤波后图像中，找出与所提取的噪声图像块索引位置相同的邻域窗口，将该邻域窗口的类别作为该噪声图像块的类别；

(4c)合并噪声图像中类别相同的噪声图像块，得到6类噪声图像块集合，将6类噪声图像块集合中的每一个图像块向量化，得到6类图像块向量集合，将图像块向量集合作为训练样本集合；

(5)分类学习稀疏字典：

采用稀疏字典学习方法，利用6类训练样本集合进行稀疏字典学习，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数；

(6)得到最终的训练字典：

按照下式，得到每一类训练样本集合最终的训练字典：

$Dic=\mathrm{\Φ}\·\tilde{S}$

其中，Dic表示每一类训练样本集合最终的训练字典，Φ表示基础字典，表示每一类训练样本集合最终更新后的稀疏字典；

(7)得到各类去噪估计值：

按照下式，得到每一类训练样本集合的去噪估计值：

Y＝Dic·Γ

其中，Y表示每一类训练样本集合的去噪估计值，Dic表示每一类训练样本集合最终的训练字典，Γ表示每一类训练样本集合最终更新后的稀疏编码系数；

(8)得到噪声抑制后的图像：

(8a)对步骤(1)输入的噪声图像中的任意一个像素点位置，从6类训练样本集合的去噪估计值中提取出与该像素点位置对应的所有像素值，将所有像素值的均值作为该像素点的去噪结果；

(8b)对步骤(1)输入的噪声图像中的所有位置的像素点，按照步骤(8a)的步骤，得到所有像素点的去噪结果，获得最终噪声抑制后的图像；

(9)输出噪声抑制后的图像。

2.根据权利要求1所述基于结构聚类与稀疏字典学习的医学图像降噪方法，其特征在于：步骤(2)所述的广义K均值K-SVD方法，按如下步骤进行：

第一步，设置图像块的大小为8×8，超完备字典的大小为64×256，初始化超完备字典为离散余弦变换DCT字典；

第二步，按照下式，求出稀疏编码系数：

${\mathrm{\α}}_{mn}=\mathrm{argmin}\mathrm{\μ}\left|\right|{\mathrm{\α}}_{mn}|{|}_0+||{\mathrm{D\α}}_{mn}-R_{mn}X|{|}_2^2$

其中，α_mn表示稀疏编码系数，m，n表示任意图像块在噪声图像中的坐标位置，argmin表示求最小值操作，μ表示残差控制因子，||·||₀表示取零范数操作，表示取二范数的平方操作，D表示超完备字典，R_mn表示图像块提取操作符，X表示去噪后的图像；

第三步，按照下式，计算误差矩阵：

$E_k=\{R_{mn}X-\underset{t\≠k}{\Σ}{d}_t{\mathrm{\α}}_{mn}\left(t\right)|(m,n)\∈w_k\}$

其中，E_k表示误差矩阵，k表示超完备字典的第k个原子的编号，R_mn是图像块提取操作符，m，n表示任意图像块在噪声图像中的坐标位置，X表示去噪后的图像，d_t表示超完备字典的第t个原子，α_mn表示稀疏编码系数，α_mn(t)表示稀疏编码系数α_mn的第t个元素，w_k表示使用了超完备字典的第k个原子d_k的所有图像块的坐标位置构成的集合；

第四步，按照下式，对误差矩阵进行奇异值分解：

E_k＝UΔV^T

其中，E_k表示误差矩阵，U表示左奇异矩阵，Δ表示奇异值矩阵，V表示右奇异矩阵，V^T表示右奇异矩阵V的转置；

第五步，用左奇异矩阵U的第一列更新超完备字典的第k列原子d_k，用右奇异矩阵V的第一列乘以残差矩阵E_k的最大特征值Δ(1,1)更新集合w_k对应的图像块的稀疏编码系数α_mn的第k个元素α_mn(k)；

第六步，对超完备字典中所有原子分别重复执行第三步至第五步的操作，得到更新后的超完备字典和稀疏编码系数；

第七步，重复执行第二步至第六步5次，得到最终更新后的超完备字典和最终更新后的稀疏编码系数；

第八步，按照下式，得到滤波后图像：

$Z^{\′}={(\mathrm{\λ}E+\underset{mn}{\Σ}{R}_{mn}^T{R}_{mn})}^{-1}(\mathrm{\λ}Z+\underset{mn}{\Σ}{R}_{mn}^T\hat{D}{\widehat{\mathrm{\α}}}_{mn})$

其中，Z′表示滤波后图像，λ表示拉格朗日系数，E表示单位矩阵，R_mn表示重叠块提取操作符，m，n表示任意图像块在噪声图像中的坐标位置，表示重叠块提取操作符R_mn的转置，表示对进行求逆运算，Z表示噪声图像，表示最终更新后的超完备字典，表示最终更新后的稀疏编码系数。

3.根据权利要求1所述基于结构聚类与稀疏字典学习的医学图像降噪方法，其特征在于：步骤(3a)所述的控制核回归权值公式如下：

${\mathrm{\ω}}_{a,b}=\frac{\sqrt{\det \left(M_b\right)}}{2{\mathrm{\πh}}^2}\exp \{-\frac{{(p_a-p_b)}^T{M}_b(p_a-p_b)}{2h^2}\}$

其中，ω_a,b表示控制核回归权值，a表示滤波后图像中任意的像素点，b表示以像素点a为中心的邻域窗口内其他的像素点，表示开方操作，det(·)表示求方阵的行列式操作，M_b表示由像素点b在水平和垂直方向上的估计梯度构成的对称梯度协方差矩阵，h表示平滑参数，exp(·)表示取指数操作，p_a和p_b分别表示像素点a和像素点b在滤波后图像中的位置，(·)^T表示转置操作。

4.根据权利要求1所述基于结构聚类与稀疏字典学习的医学图像降噪方法，其特征在于：步骤(3b)所述的K均值聚类的方法，按如下步骤进行：

第一步，从滤波后图像中任意选取6个邻域窗口的特征向量作为6类邻域窗口的聚类中心；

第二步，按照下式，计算滤波后图像的邻域窗口与6类聚类中心的特征向量的欧式距离：

${\mathrm{dis}}_{j,q}=\left|\right|{\stackrel{\‾}{\mathrm{\ω}}}_j-v_q\left|\right|$

其中，dis_j,q表示滤波后图像的邻域窗口与6类聚类中心的特征向量的欧式距离，j表示滤波后图像中邻域窗口的序号，q表示滤波后图像中邻域窗口的类别，q＝1,2,...,6，表示滤波后图像中第j个邻域窗口的特征向量，v_q表示滤波后图像中的第q类邻域窗口的聚类中心，||·||表示求欧式距离操作；

第三步，将滤波后图像中各个邻域窗口分配到与滤波后图像中6类聚类中心的特征向量的欧式距离最小的对应类别中；

第四步，对滤波后图像中类别相同的邻域窗口的特征向量求平均向量，作为滤波后图像中该类邻域窗口新的聚类中心；

第五步，将第二步至第四步重复进行100次，得到滤波后图像中各个邻域窗口的类别。

5.根据权利要求1所述基于结构聚类与稀疏字典学习的医学图像降噪方法，其特征在于：步骤(5)所述的稀疏字典学习方法的具体步骤如下：

第一步，分别设置一个大小为81×100的基础字典和一个大小为100×100的稀疏字典，令基础字典为离散余弦变换DCT字典，初始化稀疏字典为单位矩阵；

第二步，按照下式，计算稀疏编码系数：

$\begin{array}{cccc}\∀i& min\left|\right|{\mathrm{\γ}}_i|{|}_0& s.t.& \left|\right|y_i-{\mathrm{\ΦS\γ}}_i|{|}_2^2\≤{\mathrm{C\σ}}^2\end{array}$

其中，表示取训练样本序号中的任意值，i表示训练样本的序号，min表示求最小值操作，||·||₀表示取零范数操作，γ_i表示稀疏编码系数，s.t.表示约束条件，表示取二范数的平方操作，y_i表示由训练样本的序号对应的训练样本，Φ表示基础字典，S表示稀疏字典，C表示高斯白噪声的增益，σ²表示高斯白噪声的方差；

第三步，按照下式，得到使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合：

I_l＝{i|γ_i(l)≠0,i∈{1,2,...,N}}

其中，I_l表示使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合，i表示训练样本的序号，γ_i表示稀疏编码系数，γ_i(l)表示稀疏编码系数γ_i的第l个元素，N表示每一类训练样本的总个数；

第四步，按照下式，计算残差矩阵：

$E_l=\{y_i-\underset{f\≠l}{\Σ}{\mathrm{\Φs}}_f{\mathrm{\γ}}_i\left(f\right)|i\∈I_l\}$

其中，E_l表示残差矩阵，l表示稀疏字典的第l个原子的编号，i表示训练样本的序号，y_i表示由训练样本的序号对应的训练样本，∑表示求和操作，f表示稀疏字典的第f个原子的编号，Φ表示基础字典，s_f表示稀疏字典的第f个原子，γ_i表示稀疏编码系数，γ_i(f)表示稀疏编码系数γ_i的第f个元素，I_l表示使用了稀疏字典的第l个原子的所有训练样本的序号构成的集合；

第五步，按照下式，计算残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影：

z＝E_l·g

其中，z表示残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影，E_l表示残差矩阵，g表示集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量；

第六步，按照下式，计算投影z的稀疏编码系数：

其中，表示投影z的稀疏编码系数，argmin表示求最小值操作，表示取二范数的平方操作，z表示残差矩阵在集合I_l对应的训练样本的稀疏编码系数中第l个元素构成的系数向量上的投影，Φ表示基础字典，s.t.表示约束条件，||·||₀表示取零范数操作，L表示稀疏字典原子的稀疏度，L＝24；

第七步，按照下式，对投影z的稀疏编码系数进行规则化：

其中，表示对投影z的稀疏编码系数进行规则化后的稀疏编码系数，||·||₂表示取二范数操作，Φ表示基础字典；

第八步，按照下式，对稀疏字典的第l个原子进行更新：

其中，表示对稀疏字典的第l个原子更新后的新原子，表示对投影z的稀疏编码系数进行规则化后的稀疏编码系数；

第九步，按照下式，对集合I_l对应的训练样本的稀疏编码系数进行更新：

${\widetilde{\mathrm{\γ}}}_i\left(l\right)={\left({\left(E_l\right)}^T\mathrm{\Φ}{\tilde{s}}_l\right)}^T$

其中，表示对集合I_l对应的训练样本的稀疏编码系数进行更新后的稀疏编码系数，i表示训练样本的序号，l表示稀疏字典的第l个原子的编号，E_l表示残差矩阵，(·)^T表示转置操作，Φ表示基础字典，表示对稀疏字典第l个原子更新后的新原子；

第十步，对稀疏字典中所有原子分别重复执行第三步至第九步的操作，得到更新后的稀疏字典和稀疏编码系数；

第十一步，将第二步至第十步重复执行15次，得到每一类训练样本集合最终更新后的稀疏字典和稀疏编码系数。