CN103218791B - 基于稀疏自适应字典的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于稀疏自适应字典的图像去噪方法，主要解决现有方法训练字典时存在过拟合和自适应性不足的缺点。其实现过程是：(1)从含噪图像中获取图像块并列化为向量组成训练数据集；(2)利用训练数据集迭代地训练字典，在迭代过程中，将上次迭代所得的字典作为本次迭代的基字典，迭代完成后得到最终的字典和训练数据集在该字典上的编码系数矩阵；(3)利用训练所得的字典和编码系数矩阵得到去噪后的数据集；（4）利用去噪后的数据集重构出去噪后的图像。本发明方法训练出来的字典具有稀疏性和较好的自适应性，提高了图像去噪的效果，可用于对自然图像和医学CT图像的去噪。

Description

基于稀疏自适应字典的图像去噪方法

技术领域

本发明属于图像处理技术领域，特别涉及含加性噪声的图像去噪方法，可用于自然图像和医学CT图像的去噪。

背景技术

图像去噪是图像处理领域中的非常重要的一个问题。在采集获取图像过程，由于设备本身的问题以及在传输过程中的干扰，使得获取到的图像难免含有噪声，从而降低了图像的质量，对后续处理都带来了影响。因此，图像去噪是必要的，也是非常有意义的。在实际当中，大部分噪声可以近似的认为是加性的高斯白噪声，因此去除含噪图像中的高斯白噪声成为图像去噪领域中一个非常重要的方向。

传统的去噪方法一般可以分为两类，一类是基于空域的方法，一类是基于变换域的方法。基于空域的去噪方法中，比较经典的方法有，中值滤波，高斯滤波，双边滤波，以及近几年提出的非局部均值等方法。传统经典的方法的共性是利用图像灰度值在局部窗口内的连续性这个特征来修正像素值，从而进行去噪。这些方法的缺点是在去除噪声的同时模糊了图像的细节信息，例如图像的边缘，纹理等。非局部均值的方法的提出在一定程度上弥补了局部的缺点，它不仅利用图像局部的连续性，也利用图像内容的冗余性，从而取得了一定的效果。基于变换域的去噪方法中，具有代表意义的是小波域的各种去噪方法，小波对点的奇异性表示比较好，但由于它缺少方向性，对图像边缘、轮廓等线性奇异性的结构特征表示不是太好，为了解决这个问题，一些新的具有多尺度多方向特性的变换应运而生，如：Brushlet变换、Curvelet变换、Contourlet变换等。这些变换在一定程度上解决了小波的缺少方向性问题，但也不足以对图像进行通用的表示。

基于冗余字典稀疏表示的去噪方法可以归为基于变换域的图像去噪方法，但是稀疏表示与小波以及多尺度不同的是，它是利用图像信号的稀疏性这个特征，来区分噪声和信号，从而进行图像去噪。而在进行稀疏表示时，能否得到好的字典是十分关键的。在训练字典的方法中，KSVD算法是典型代表。该算法利用离散余弦变换DCT初始化字典，通过含噪图像中的数据来学习训练，得到更新后的字典。因为是在含噪声图像块组成的数据集训练字典，这使得训练出来的字典具有噪声，从而影响了去噪效果。同时该方法中字典学习速度较慢，在训练数据集较少时存在过学习的情况。在KSVD算法基础上提出来的S-KSVD算法虽然提升了速度，同时减轻了KSVD过学习的缺点，但却存在去噪效果相对较差的缺点。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于稀疏自适应字典的图像去噪方法，以提高训练字典的速度和对图像的去噪效果。

实现本发明目的的技术思路是：根据数据和字典都具有稀疏性的特点，在字典学习过程中，更新原子时将新原子限制在已得字典的列空间内。

1.一种基于稀疏自适应字典的图像去噪方法，包括如下步骤：

(1)对含噪图像Ω中的任一点i，以其为中心选取大小为的邻域块，将邻域块列化为向量记为y_i，得到数据集Y＝{y_i|i∈1,2,...,N}，其中N为邻域块的个数；

(2)利用数据集Y训练字典：

2a)设迭代次数J为15，字典中原子的个数M＝4n，初始字典D⁽⁰⁾为大小n×M的离散余弦字典,循环变量k＝1；

2b)对数据集Y中的每一列信号y_i，求解它在字典D^(k-1)上的编码系数α_i,即利用正交匹配追踪OMP方法求解如下优化问题得出α_i：

$\∀i\∈\{\mathrm{1,2},\·\·\·,N\},\underset{{\mathrm{\α}}_i}{\min }{\left|\right|{\mathrm{\α}}_i\left|\right|}_{0}s.t.{\left|\right|y_i-D^{(k-1)}{\mathrm{\α}}_i\left|\right|}_2^2\≤n{\left(\mathrm{C\σ}\right)}^2$

其中，C为一个已知常数系数，σ为噪声方差，||α_i||₀表示向量α_i的零范数，s.t.表示对该问题的约束条件，记编码系数矩阵为

A^(k-1)＝(α₁,α₂,...,α_N)；

2c)利用编码系数矩阵A^(k-1)依次更新字典D^(k-1)中的所有原子，设循环变量l＝1：

2c1)找出在稀疏表示中使用了原子d_l的所有信号，这些信号的编号合w_l，即w_l＝{i|α_i(l)≠0,i∈{1,2,...,N}}，其中，α_i(l)表示向量α_i的第l个元素；

2c2)对于i∈w_l，计算信号y_i在d_l上的表示残差,：

$e_i^l=y_i-\underset{m\≠l}{\mathrm{\Σ}}{d}_m{\mathrm{\α}}_i\left(m\right),$

用所有的i∈w_l构成残差矩阵E_l；

2c3)计算残差矩阵E_l在相应的系数向量g上的投影d′：

d′＝E_l·g，

其中，j_p表示集合w_l中的第p个元素，p∈{1,2,...,|w_l|}，|w_l|表示集合w_l中元素的个数；

2c4)利用字典D^(k-1)对投影d′进行稀疏编码得到编码系数β，即利用正交匹配追踪OMP方法求解优化问题：

$\underset{\mathrm{\β}}{\min }{\left|\right|d^{\′}-D^{(k-1)}\mathrm{\β}\left|\right|}_2^{2}s.t.{\left|\right|\mathrm{\β}\left|\right|}_0\≤L,$

其中L＝10为新原子的稀疏度；

2c5)利用字典D^(k-1)和d′的编码系数β得到新原子：d_new＝D^(k-1)β，以及信号{y_i|i∈w_l}在新原子上的系数g_new＝(E_l)^Td_new；

2c6)更新原子d_l＝d_new，同时更新{y_i|i∈w_l}在d_l上的表示系数为g_new，即α_i(l)＝g_new(i)，i∈w_l；

2c7)判断本次字典更新是否已经完成：若l＝M，则本次字典更新已完成，得到了字典D^(k)和系数矩阵A^(k)，执行步骤2d）,否则将l增加1，返回步骤2c1）；

2d)判断字典训练是否完成：若k＝J，则字典训练已经完成，得到字典D^(J)和系数矩阵A^(J)，执行步骤（3），否则将k增加1，返回步骤2b；

(3)利用第（2）步所得的字典D^(J)和系数矩阵A^(J)，计算去噪后的数据集：

$\hat{Y}=D^{\left(J\right)}{A}^{\left(J\right)};$

(4)利用去噪后的数据集得到去噪后的图像。

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

1.本发明方法训练所得的字典与KSVD方法训练所得的字典相比，得到的原子是稀疏的，在能够防止训练数据较少时KSVD方法过学习而拟合噪声的缺点，进而减少了字典中所含的噪声。

2.与S-KSVD算法相比，本发明方法将上一次迭代中所得字典作为基字典，克服了S-KSVD方法所得字典自适应性不足的缺点，相比S-KSVD，本发明在字典稀疏性和自适应性间取得了很好的折中，因而去噪效果更好。

3．从去噪后图像的视觉效果上看，本发明比原始的KSVD方法和S-KSVD方法去噪效果都要好，尤其是相比S-KSVD方法本发明去噪后的图像峰值信噪比PSNR提高最为明显。

综上，本发明不同于KSVD的只使用数据的稀疏性的方法，也不同于S-KSVD的认为字典本身在一个固定的基字典下具有稀疏性的方法，本发明方法中的基字典是在学习中不断变化的，因而在保证字典稀疏性的同时有更好的自适应性。

附图说明

图1是本发明的流程图；

图2是本发明实验中使用的第一组三维CT测试图像中序号为70的一幅图像；

图3是对图2加入σ＝30的高斯白噪声后的含噪图像；

图4是本发明实验中使用的第二组三维CT测试图像中序号为98的一幅图像；

图5是对图4加入σ＝30的高斯白噪声后的含噪图像；

图6是使用现有KSVD方法对图3去噪所得的结果图像；

图7是使用现有S-KSVD方法对图3去噪所得的结果图像；

图8是使用本发明对图3去噪所得的结果图像；

图9是使用现有KSVD方法对图3训练所得字典的各个原子组成图像块后的灰度图像；

图10是使用本发明对图3训练所得字典的各个原子组成图像块后的灰度图像；

图11是使用现有KSVD方法对图5训练所得字典的各个原子组成图像块后的灰度图像；

图12是使用现有S-KSVD方法对图5去噪所得的结果图像；

图13是使用本发明对图5去噪所得的结果图像；

图14是使用现有KSVD方法对图5训练所得字典的各个原子组成图像块后的灰度图像；

图15是使用本发明从图5训练所得字典的各个原子组成图像块后的灰度图像。

具体实施方式

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

步骤1在含噪声图像上获取训练数据集Y：

记输入图像为从噪声图像的左上角第一个点开始，自上而下，自左到右，逐点扫描，依次取以当前点i为中心的大小为的图像块；然后将图像块列化为向量y_i，并组成训练数据集Y＝{y_i|i∈1,2,...,N}，其中N为图像块的个数。

步骤2利用数据集Y训练字典：

2a)设迭代次数J为15，字典中原子的个数M＝4n，初始字典D⁽⁰⁾的大小为n×M的离散余弦字典，循环变量k＝1；

2b)对数据矩阵Y中的每一列信号y_i，i∈{1,2,...,N}，求解它在字典D^(k-1)上的编码系数α_i，即利用正交匹配追踪OMP方法求解如下优化问题得出编码系数α_i：

$\underset{{\mathrm{\α}}_i}{\min }{\left|\right|{\mathrm{\α}}_i\left|\right|}_{0}s.t.{\left|\right|y_i-D^{(k-1)}{\mathrm{\α}}_i\left|\right|}_2^2\≤n{\left(\mathrm{C\σ}\right)}^2,i\∈\{\mathrm{1,2},\·\·\·,N\}$

其中，||α_i||₀表示α_i的零范数，C＝1.15，σ为噪声方差，记编码系数矩阵(α₁,α₂,...,α_N)为A^(k-1)；

2c)利用编码系数矩阵A^(k-1)更新字典D^(k-1)：

2c1)找出在稀疏表示中使用了原子d_l的所有信号，用这些信号的编号组成集合w_l，即w_l＝{i|α_i(l)≠0,i∈{1,2,...,N}}，其中，α_i(l)表示向量α_i的第l个元素；

2c2)对于i∈w_l，计算信号y_i在d_l上的表示残差：用所有的i∈w_l构成残差矩阵E_l；

2c3)计算残差矩阵E_l在相应的系数向量g上的投影d′：

d′＝E_l·g，

其中，j_p表示集合w_l中的第p个元素，,，|w_l|表示集合w_l中元素的个数；

2c4)利用字典D^(k-1)对投影d′进行稀疏编码得到编码系数β，即利用正交匹配追踪OMP方法求解优化问题：

$\underset{\mathrm{\β}}{\min }{\left|\right|d^{\′}-D^{(k-1)}\mathrm{\β}\left|\right|}_2^{2}s.t.{\left|\right|\mathrm{\β}\left|\right|}_0\≤L,$

其中L＝10为新原子的稀疏度；

2c5)利用字典D^(k-1)和投影d′的编码系数β得到新的原子：d_new＝D^(k-1)β，

以及信号{y_i|i∈w_l}在新原子上的系数g_new＝(E_l)^Td_new；

2c6)更新原子d_l，即将字典D^(k-1)中的原子d_l替换为d_new，同时将{y_i|i∈w_l}在d_l上的表示系数g替换为g_new，即α_i(l)＝g_new(i)，i∈w_l；

2c7)判断本次字典更新是否已经完成：若l＝M，则字典D^(k-1)中的所有原子都已更新，得到字典D^(k)和系数矩阵A^(k)，执行步骤2d）,否则将l增加1，返回步骤2c1）；

2d)判断字典训练是否已全部完成：若k＝J，则字典训练已经完成，得到字典D^(J)和系数矩阵A^(J)，执行步骤（3），否则将k增加1，返回步骤2b。

步骤3利用步骤2所得的字典D^(J)和系数矩阵A^(J)，计算去噪后的数据集：

$\hat{Y}=D^{\left(J\right)}{A}^{\left(J\right)};$

步骤4利用去噪后的数据集得到去噪后的图像：

4a)对含噪图像Ω中的任意一点j，从去噪后的数据集中取出与该点对应的所有像素值，取它们的均值作为该点的去噪结果；

4b)对含噪图像Ω中的各个点都按照步骤4a)得到含噪图像Ω中所有点的去噪结果，进而得到去噪后的整幅图像。

本发明的效果可通过以下实验进一步证明：

一、实验条件

1.实验条件：实验所使用的图像为三维CT图像，共两组：第一组是某女性踝关节CT扫描图Vis.F.Ankle，大小为260×190×150，第二组是某男性头部CT扫描图Vis.F.Head，大小为260×190×206。图2和图3是分别从第一组和第二组CT图像中抽取的序号分别为70和98的两幅图像。实验中的关键参数设置为：n＝64，J＝15，C＝1.15，L＝10。实验中，各种方法都是使用MATLAB语言编程实现。

2.实验内容与结果：

实验1，在图2所示图像上加入标准差σ＝30的加性高斯白噪声，得到含噪图像如图3；在上述实验条件下，分别使用KSVD方法和本发明方法对图3训练字典，所得的字典分别如图9和图10；在上述实验条件下，分别使用KSVD方法、S-KSVD方法和本发明方法对图3去噪，得到去噪结果图像分别为图6、图7和图8。

比较图7和图8可以看出，S-KSVD方法所得的图像细节信息损失严重，图像边缘明显变模糊，本发明方法在去除噪声的同时能够保持图像细节，从视觉效果上看明显优于S-KSVD方法。

比较图6和图8可以看出，本发明方法在图3上的去噪效果与KSVD方法相当，但是比较图9和图10可以看出，KSVD方法训练所得的字典含有较多的噪声，本发明方法训练所得的字典与之相比噪声明显较少，克服了KSVD方法易过拟合噪声的缺点。

实验2，在图3所示图像上加入标准差σ＝30的加性高斯白噪声，得到含噪图像如图5所示；在上述实验条件下，分别使用KSVD方法和本发明方法对图5训练字典，所得的字典分别如图14和图15；在上述实验条件下，分别使用KSVD方法、S-KSVD方法和本发明方法对图5去噪，得到去噪结果图像分别为图11、图12和图13。

比较图12和图13可以看出，S-KSVD方法所得的图像细节信息损失严重，图像边缘明显变模糊，本发明方法在去除噪声的同时能够保持图像细节，从视觉效果上看明显优于S-KSVD方法。

比较图11和图13可以看出，本发明方法在图5上的去噪效果略优于KSVD；比较图14和图15可以看出，KSVD方法训练所得字典含有较多的噪声，本发明方法训练所得字典与之相比噪声明显较少，克服了KSVD方法易过拟合噪声的缺点。

实验3，在上述第一组三维CT图像Vis.F.Ankle中，分别加入标准差σ＝10,20,30的加性高斯白噪声，分别使用KSVD方法、S-KSVD方法和本发明方法对该组图像中的每一幅分别处理，计算去噪后的三维CT图像整体的峰值信噪比PSNR，结果如表1。

实验4，在上述第二组三维CT图像Vis.F.Head中，分别加入标准差σ＝10,20,30的加性高斯白噪声，分别使用KSVD方法、S-KSVD方法和本发明方法对该组图像中的每一幅分别处理，计算去噪后的三维CT图像整体的峰值信噪比PSNR，结果如表1。

表格1三种方法所得结果图像的峰值信噪比PSNR

从表1可以看出，本发明方法在三维CT图像上的去噪效果明显优于KSVD和S-KSVD，峰值信噪比PSNR有显著的提高。

综上，本发明方法与其余两种方法相比具有明显的优势：相比KSVD方法和S-KSVD方法，本发明方法去噪性能更好，所得结果图像的PSNR有显著提高，相比S-KSVD方法本发明所得图像的PSNR提高最为显著。本发明方法训练出来的字典即是稀疏的又具有较好的自适应性，克服了KSVD方法训练字典时易过拟合噪声和S-KSVD方法训练字典时自适应性不足的缺点。本发明方法在字典的稀疏性和自适应性之间找到了很好的折中，无论在客观指标还是主观效果上，都表现出了很好的性能。

Claims (2)

1.一种基于稀疏自适应字典的图像去噪方法，包括如下步骤：

(1)对含噪图像Ω中的任一点i，以其为中心选取大小为的邻域块，将邻域块列化为向量记为y_i，得到数据集Y＝{y_i|i∈1,2,...,N}，其中N为邻域块的个数，n＝64；

(2)利用数据集Y训练字典：

2a)设迭代次数J为15，字典中原子的个数M＝4n，初始字典D⁽⁰⁾为大小n×M的离散余弦字典,循环变量k＝1；

2b)对数据集Y中的每一列信号y_i，求解它在字典D^(k-1)上的编码系数α_i,即利用正交匹配追踪OMP方法求解如下优化问题得出α_i：

$\begin{array}{ccc}\∀i\∈\{\mathrm{1,2},...,N\},\underset{{\mathrm{\α}}_i}{\min }{\left|\right|{\mathrm{\α}}_i\left|\right|}_0& s.t.& {\left|\right|y_i-D^{(k-1)}{\mathrm{\α}}_i\left|\right|}_2^2\≤{\left(\mathrm{C\σ}\right)}^2\end{array}$

其中，C为一个已知常数系数，σ为噪声方差，||α_i||₀表示向量α_i的零范数，s.t.表示对该问题的约束条件，记编码系数矩阵为

A^(k-1)＝(α₁,α₂,...,α_N)；

2c)利用编码系数矩阵A^(k-1)依次更新字典D^(k-1)中的所有原子，设循环变量l＝1：

2c1)找出在稀疏表示中使用了原子d_l的所有信号，这些信号的编号组成集合w_l，即w_l＝{i|α_i(l)≠0,i∈{1,2,...,N}}，其中，α_i(l)表示向量α_i的第l个元素；

2c2)对于i∈w_l，计算信号y_i在d_l上的表示残差

$e_i^l=y_i-\underset{m\≠l}{\mathrm{\Σ}}{d}_m{\mathrm{\α}}_i\left(m\right),$

用所有的i∈w_l构成残差矩阵E_l；

2c3)计算残差矩阵E_l在相应的系数向量g上的投影d′：

d′＝E_l·g，

其中，j_p表示集合w_l中的第p个元素，p∈{1,2,...,|w_l|}，|w_l|表示集合w_l中元素的个数；

2c4)利用字典D^(k-1)对投影d′进行稀疏编码得到编码系数β，即利用正交匹配追踪OMP方法求解优化问题：

$\begin{array}{ccc}\underset{\mathrm{\β}}{\min }{\left|\right|d^{\′}-D^{(k-1)}\mathrm{\β}\left|\right|}_2^2& s.t.& {\left|\right|\mathrm{\β}\left|\right|}_0\≤L\end{array},$

其中L＝10为新原子的稀疏度；

2c5)利用字典D^(k-1)和d′的编码系数β得到新原子：d_new＝D^(k-1)β，以及信号{y_i|i∈w_l}在新原子上的系数g_new＝(E_l)^Td_new；

2c6)更新原子d_l＝d_new，同时更新{y_i|i∈w_l}在d_l上的表示系数为g_new，即α_i(l)＝g_new(i)，i∈w_l；

2c7)判断本次字典更新是否已经完成：若l＝M，则本次字典更新已完成，得到了字典D^(k)和系数矩阵A^(k)，执行步骤2d),否则将l增加1，返回步骤2c1)；

2d)判断字典训练是否完成：若k＝J，则字典训练已经完成，得到字典D^(J)和系数矩阵A^(J)，执行步骤(3)，否则将k增加1，返回步骤2b；

(3)利用第(2)步所得的字典D^(J)和系数矩阵A^(J)，计算去噪后的数据集：

$\hat{Y}=D^{\left(J\right)}{A}^{\left(J\right)};$

(4)利用去噪后的数据集得到去噪后的图像。

2.根据权利要求1所述的基于稀疏自适应字典的图像去噪方法，其特征在于步骤(4)，按如下步骤进行：

4a)对含噪图像Ω中的任意一点j，从去噪后的数据集中取出与该点对应的所有像素值，取它们的均值作为该点的去噪结果；

4b)对含噪图像Ω中的各个点都按照步骤4a)得到含噪图像Ω中所有点的去噪结果，进而得到去噪后的整幅图像