CN102063708B - 基于Treelet和非局部均值的图像去噪 - Google Patents

Abstract

本发明公开了一种基于Treelet和非局部均值的图像去噪方法，主要解决现有的非局部均值方法对被高噪声腐蚀的自然图像去噪效果不佳的问题。其实现步骤是：(1)对输入的图像计算其协方差矩阵；(2)对协方差矩阵计算图像Treelet变换的尺度向量；(3)对图像逐像素取滑窗与尺度向量相乘，得到特征向量函数；(4)由特征向量函数对图像逐像素滤波，得到去噪后的图像。本发明具有对高噪声腐蚀的自然图像去噪效果好的优点，能恢复出图像原有的特征，可用于变化检测、目标识别时对图像的预处理。

Description

基于Treelet和非局部均值的图像去噪

技术领域

本发明属于图像处理技术领域，涉及高噪声情况下的自然图像的去噪，可用于开展森林资源调查、土地利用与覆盖的变化检测、环境变化评估、灾害评估、城市规划、国防军情监控、医学影像、天文学影像等领域的数字图像预处理。

背景技术

图像去噪的主要目的是解决实际图像由于噪声干扰导致的图像质量下降问题。通过去噪可以提高图像质量，增大信噪比，更好地体现图像所携带的信息，因此图像去噪技术在很多领域中都占据着很重要的位置。

根据图像的特点和噪声的统计特征，多年来已经提出了很多的图像去噪方法，现有消除噪声即滤波的已有方法主要分为空间域滤波和频率域滤波。图像的空间域滤波方法是直接对图像的灰度做运算，图像的频率域滤波方法则是在某种变换域内对图像的变换系数进行运算，然后反变换回图像空间域的。其中空间域的滤波因不需要对原图像进行变换和反变换，方法直接简洁，在一定领域内较为常用。常用的空间滤波方法有拟合图像和平滑图像的方法。这些方法以平滑数据的方式去除噪声，通常也会模糊数据本身，使得纹理等细节结构无法得到较好的保持。非局部均值去噪方法利用图像的信息冗余特性，通过对待去噪像素点的邻域窗口内的所有像素点进行加权平均取得新的像素灰度值，并且在对每个像素点进行加权平滑的过程中考虑了局部结构的相似性，充分利用了图像的冗余性以及局部结构信息，因此取得到了较好的滤波效果。

在非局部均值方法中其每一个相似块实际就是一种无序含噪的高维数据，若相似窗大小为5×5像素，则图像中每个像素的邻域特征向量为25维，因此计算特征向量间的相似性需要大量的存储空间和计算时间，算法的复杂度较高。Huang等人和Lee等人发现图像中的强度数据在低维空间中更为集中，这意味着图像在低维空间中更能体现其相关性。Azzabou等人最早提出将非局部均值方法中的相似窗的邻域特征向量投影至其低维的子空间中计算，通过主成分分析PCA得到邻域观察值向量的低维子空间，并将邻域特征向量间的相似性权重在低维空间中计算，降低了非局部均值方法的时间消耗，且精度有所提高，但主成分分析获得的低维特征是针对高维数据的特定方向进行投影，当高维数据的分布不明显时，获得的低维特征将会产生误差。Tolga等人同样提出了一种基于字典的主成分分析方法，在提取的前d个低维主成份空间中计算每个像素点的权重取得了更优的去噪效果，但需要对获得的字典耗费大量时间进行训练，算法复杂度高。

在本发明中我们引入了一种新的自适应多尺度分析和表示方法-Treelet，该方法通过对高维数据的逐层降维分解更能反映数据潜在的结构以及变量之间的相关性，对于高维无序、含噪数据的分析要优于PCA等降维方法，因此本发明将Treelet与非局部均值方法相结合，提出了一种基于Treelet和非局部均值的图像去噪方法，该方法相比前两种方法在高噪声下的去噪效果更优。

发明内容

本发明的目的在于针对非局部去噪中高维数据的计算问题，提出一种基于Treelet和非局部均值的图像去噪方法，以减少了计算复杂度，克服现有方法对自然图像去噪产生的视觉模糊现象，提高去噪后图像的可视性。

为实现上述目的，本发明的去噪方法包括如下步骤：

(1)对输入加噪图像逐像素取5×5的滑窗N_i，计算其协方差矩阵

(2)根据滑窗N_i的协方差矩阵

计算图像的Treelet变换的尺度向量Φ：

(2a)依据滑窗N_i的协方差矩阵

计算图像的相似度矩阵

${\hat{M}}_{\mathrm{mn}}=\left|\frac{{\widehat{\mathrm{\Σ}}}_{\mathrm{mn}}}{\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{mm}}{\widehat{\mathrm{\Σ}}}_{\mathrm{nn}}}}\right|+\left|{\widehat{\mathrm{\Σ}}}_{\mathrm{mn}}\right|$

其中，m和n分别为协方差阵

的位置索引；

为

的第m行第n列的矩阵元素，

和分别为的第m行第m列的矩阵元素以及第n行第n列的矩阵元素；

(2b)定义Treelet变换分解的层数l＝0，1，…，20，初始化和变量的下标集：δ⁽⁰⁾＝{1，2，L，25}，当l＝0时，初始化Treelet分解的协方差矩阵为：

相似度矩阵为：

初始化基矩阵B₀为25×25的单位矩阵；

(2c)由当前l层Treelet分解的相似度矩阵

中找出最相似的两个变量：

$(\mathrm{\α},\mathrm{\β})=\arg \max {\hat{M}}^{\left(l\right)}$

其中，α和β分别表示当前l层Treelet分解的协方差矩阵

中的位置索引，分别对应第一主成分和第二主成分；

(2d)对当前l层的Treelet分解的协方差矩阵进行局部主成分分析变换，即通过雅克比旋转矩阵J，计算当前聚类层级l的基矩阵B_l＝B_l-1J，并更新相似度矩阵

和协方差矩阵

雅克比旋转矩阵

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& \cos \left({\mathrm{\θ}}_l\right)& L& -\sin \left({\mathrm{\θ}}_l\right)& L& 0\\ M& & M& O& M& & M\\ 0& L& \sin \left({\mathrm{\θ}}_l\right)& L& \cos \left({\mathrm{\θ}}_l\right)& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

这里的旋转角度θ_l由

且|θ_l|≤π/4计算得到；

(2e)旋转角度后α和β分别代表

第一和第二主成分索引位置，尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，定义当前l层级的尺度向量Φ是尺度函数φ_l和上一层的尺度向量{φ_l-1，jt}_jt≠α，β的合集，将差变量β的下标从和变量下标集δ中去除；

(2f)重复步骤(2c)至步骤(2e)直至分解l＝20层，得到最终Treelet分解的基矩阵：

B＝[Φ，ψ₁...，ψ_l]^T

其中，Φ和ψ分别是正交基矩阵的尺度函数向量和细节函数向量；

(3)对图像每个像素点i取5×5像素的滑动窗N_i与尺度向量Φ相乘，得到尺度空间的特征向量f(N_i)：

$f\left(N_i\right)=\underset{p=1}{\overset{5}{\mathrm{\Σ}}}\⟨y\left(N_i\right),{\mathrm{\φ}}_p\⟩{\mathrm{\φ}}_p$

其中<y(N_i)，φ_p>代表两个向量的内积；

(4)由尺度空间的特征向量f(N_i)对图像逐像素进行非局部均值滤波，得到去噪后的图像。

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

a、本发明由于利用Treelet变换提取了相似窗内高维邻域向量的主要特征向量，有效地减少了计算复杂度；

b、本发明由于利用Treelet分解的主要特征向量取代高维无序数据，即原始数据与Treelet分解的基矩阵中的尺度函数相乘，得到该原始数据的主要特征向量表示，减少含噪数据对图像去噪的影响，使得在高噪声情况下的去噪效果优于传统非局部均值去噪方法。

实验结果和数据表明在高噪声情况下，本发明基于Treelet和非局部均值的图像去噪效果要优于现有的非局部均值去噪的两种方法。

附图说明

图1是本发明的实现流程图；

图2是本发明的实验图像；

图3是本发明的去噪结果图像。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，输入加噪图像，计算其协方差矩阵

对输入图像中的每个像素i将以其为中心的5×5像素大小窗口的邻域N_i内的像素灰度值向量标记为y(N_i)，计算此邻域的协方差矩阵

$\widehat{\mathrm{\&Sigma;}}=\frac{1}{\left|\mathrm{\&Omega;}\right|}\underset{i\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}(y\left(N_i\right)-\stackrel{\&OverBar;}{y}){(y\left(N_i\right)-\stackrel{\&OverBar;}{y})}^T---\left(1\right)$

其中，Ω为输入图像中所有像素点的集合，|Ω|为输入图像的所有像素的总个数，

为输入图像的像素灰度均值向量；

步骤2，根据滑窗Ni的协方差矩阵

计算图像Treelet变换的尺度向量Φ。

(2a)依据滑窗N_i的协方差矩阵计算图像的相似度矩阵

${\hat{M}}_{\mathrm{mn}}=\left|\frac{{\widehat{\mathrm{\&Sigma;}}}_{\mathrm{mn}}}{\sqrt{{\widehat{\mathrm{\&Sigma;}}}_{\mathrm{mm}}{\widehat{\mathrm{\&Sigma;}}}_{\mathrm{nn}}}}\right|+\left|{\widehat{\mathrm{\&Sigma;}}}_{\mathrm{mn}}\right|---\left(2\right)$

其中，m和n分别为协方差阵的位置索引；为的第m行第n列的矩阵元素，

和

分别为

的第m行第m列的矩阵元素以及第n行第n列的矩阵元素；

(2b)定义Treelet变换分解的层数l＝0，1，…，20，初始化和变量的下标集：δ⁽⁰⁾＝{1，2，L，25}，当l＝0时，初始化协方差矩阵为

相似度矩阵为

初始化基矩阵B₀为25×25的单位矩阵：

$B_0=\left[\begin{array}{cccc}1& 0& L& 0\\ 0& 1& L& 0\\ M& M& O& M\\ 0& 0& L& 1\end{array}\right]---\left(3\right)$

(2c)由当前l层Treelet分解的相似度矩阵

中找出最相似的两个变量：

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \max {\hat{M}}^{\left(l\right)}---\left(4\right)$

其中，α和β分别表示当前l层Treelet分解的协方差矩阵

中的位置索引，分别对应第一主成分和第二主成分；

(2d)对当前l层的Treelet分解的协方差矩阵

进行局部主成分变换，即通过雅克比旋转矩阵J，计算当前聚类层级l的基矩阵B_l＝B_l-1J，并更新相似度矩阵

和协方差矩阵

雅克比旋转矩阵

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& \cos \left({\mathrm{\&theta;}}_l\right)& L& -\sin \left({\mathrm{\&theta;}}_l\right)& L& 0\\ M& & M& O& M& & M\\ 0& L& \sin \left({\mathrm{\&theta;}}_l\right)& L& \cos \left({\mathrm{\&theta;}}_l\right)& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]---\left(5\right)$

这里的旋转角度θ_l由

且|θ_l|≤π/4计算得到；

(2e)旋转角度后

α和β分别代表

第一和第二主成分索引位置，尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，定义当前l层级的尺度向量Φ是尺度函数φ_l和上一层的尺度向量{φ_l-1，jt}_jt≠α，β的合集，将差变量β的下标从和变量下标集δ中去除；

(2f)重复步骤(2c)至步骤(2e)直至分解l＝20层，得到最终Treelet分解的基矩阵：

B＝[Φ，ψ₁...，ψ_l]^T    (6)

其中，Φ和ψ分别是正交基矩阵的尺度函数向量和细节函数向量，尺度函数向量为原始数据主要成分的基表示，细节函数向量为原始数据细节成分的基表示，原始数据与对应的向量相乘分别得到主要成分的特征向量和细节成分的特征向量。

步骤3，计算5×5像素的滑动窗内像素的尺度空间的特征向量f(N_i)。

对图像每个像素点i取5×5像素的滑动窗N_i与尺度向量Φ相乘，得到尺度空间的特征向量f(N_i)：

$f\left(N_i\right)=\underset{p=1}{\overset{5}{\mathrm{\&Sigma;}}}\&lang;y\left(N_i\right),{\mathrm{\&phi;}}_p\&rang;{\mathrm{\&phi;}}_p---\left(7\right)$

其中，<y(N_i)，φ_p>代表两个向量的内积，通过上式提取了原始数据的5个主要的特征向量来取代原始数据，达到降低数据维数，减少含噪数据对去噪效果的影响。

步骤4，由尺度空间中的特征向量f(N_i)对图像逐像素进行非局部均值滤波。

(4a)在以像素点i为中心的21×21像素大小的搜索区域S_i内，逐个计算邻域N_j的特征向量集f(N_j)与(N_i)的欧氏距离：||f(N_i)-f(N_j)||²；

(4b)将上述欧式距离代入如下公式，得到像素点i去噪后的灰度值v’(i)：

$v^{\&prime;}\left(i\right)=\underset{y\&Element;S_i}{\mathrm{\&Sigma;}}\frac{1}{Z\left(i\right)}{e}^{-{\left|\right|f\left(N_i\right)-f\left(N_j\right)\left|\right|}^2/h^2}v\left(i\right)---\left(8\right)$

其中，归一化系数Z(i)由如下公式计算

$Z\left(i\right)=\underset{j\&Element;S_i}{\mathrm{\&Sigma;}}{e}^{-\frac{{\left|\right|f\left(N_i\right)-f\left(N_j\right)\left|\right|}^2}{h^2}}---\left(9\right)$

本发明的效果可以进一步通过以下内容进行说明：

1.实验数据

实验所使用的输入图像如图2所示，图2(a)为Lena图像，图2(b)为Peppers图像，图2(c)为Barara图像，三幅图像均为自然图像，大小均为512×512像素，灰度级为256。

2.实验评价指标

图像去噪效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，图像质量非常好，感觉很清晰则去噪效果好，反之效果较差。在客观上评价一幅图像的去噪效果本发明采用峰值信噪比PSNR来衡量。令真实图像为v(i)，去噪后的图像为v’(i)，令v_max＝max{v(i)，i∈Ω}。则峰值信噪比按如下公式计算：

$\mathrm{PSNR}=10\log \left[\frac{v_{\max }^2}{\frac{1}{\left|\mathrm{\&Omega;}\right|}\underset{i\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{[v^{\&prime;}\left(i\right)-v\left(i\right)]}^2}\right]---\left(10\right)$

评价去噪效果的方法是在一幅清晰的图像上加入随机高斯白噪声，然后在加噪的图像上进行去噪的实验。

2.实验结果与分析

采用以上参数的选取方法，表1列出了在不同噪声水平下本发明与文献“Buades A.，Coll B.，Morel J.M..A non-local algorithm for imagedenoising.Conf Computer Vision andPattern Recognition，CVPR2005，60-65.”中非局部均值去噪方法(NLM)以及文献“Tasdizen T..Principal components for non-local means image denoising.Intern ConfImage Processing，ICIP2008，73-77.”中基于PCA的非局部均值去噪方法(PNLM)进行对比实验，其中NLM方法的平滑参数h根据文献取为10σ～15σ中的最优值，PNLM方法采用文献中参数选取的方法，本发明中相似窗大小均为5×5，搜索窗为21×21，尺度空间维数为5。

图3为图2(a)的Lena图像在噪声水平为50时三种方法的实验结果，其中图3(a)是现有NLM算法的去噪结果，图3(b)是现有PNLM算法的去噪结果，图3(c)是本发明的去噪结果，从图3可以看出，本发明的去噪效果明显优于现有的两种非局部均值去噪算法。

表1为噪声水平为10～50时三种方法的实验评价指标。

表1NLM、PNLM和本发明方法在不同噪声水平下的去噪结果(PSNR)

从表1中可以看到本发明方法在去噪效果上略优于PNLM方法，但更优于NLM方法，尤其在高噪声情况下优势更加明显，如Lena图像在噪声标准差为50时NLM和PNLM方法的PSNR值分别为26.24和27.47，本发明方法比NLM方法高出了1db以上。

Claims (2)

1.一种基于Treelet和非局部均值的图像去噪方法，包括如下步骤：

(1)对输入加噪图像逐像素取5×5的滑窗N_i，计算其协方差矩阵 

(2)根据滑窗N_i的协方差矩阵 

计算图像的Treelet变换的尺度向量Φ：

(2a)依据滑窗N_i的协方差矩阵 

计算图像的相似度矩阵 

其中，m和n分别为协方差阵 

的位置索引； 

为 

的第m行第n列的矩阵元素， 

和 分别为 

的第m行第m列的矩阵元素以及第n行第n列的矩阵元素；

(2b)定义Treelet变换分解的层数l＝0，1，…，20，初始化和变量的下标集：δ⁽⁰⁾＝{1，2，…，25}，当l＝0时，初始化Treelet分解的协方差矩阵为： 

相似度矩阵为： 

初始化基矩阵B₀为25×25的单位矩阵；

(2c)由当前l层Treelet分解的相似度矩阵 

中找出最相似的两个变量：

其中，α和β分别表示当前l层Treelet分解的协方差矩阵 

中的位置索引，分别对应第一主成分和第二主成分；

(2d)对当前l层的Treelet分解的协方差矩阵 

进行局部主成分分析变换，即通过雅克比旋转矩阵J，计算当前聚类层级l的基矩阵B_l＝B_l-1J，并更新相似度矩阵

和协方差矩阵

雅克比旋转矩阵 

这里的旋转角度θ_l由 

计算得到；

(2e)旋转角度后 

α和β分别代表 

第一和第二主成分索引位置，尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，定义当前l层级的尺度向量Φ是尺度函数φ_l和上一层的尺度向量 

的合集，将差变量β的下标从和变量下标集δ中去除；

(2f)重复步骤(2c)至步骤(2e)直至分解l＝20层，得到最终Treelet分解的基矩阵：

B＝[Φ，ψ₁...，ψ_l]^T

其中，Φ和ψ分别是正交基矩阵的尺度函数向量和细节函数向量；

(3)对图像每个像素点i取5×5像素的滑动窗N_i与尺度向量Φ相乘，得到尺度空间的特征向量f(N_i)：

其中<y(N_i)，φ_p>代表两个向量的内积；

(4)由尺度空间的特征向量f(N_i)对图像逐像素进行非局部均值滤波，得到去噪后的图像：

(4a)在以像素点i为中心的21×21像素大小的搜索区域S_i内，逐个计算邻域N_j的特征向量集f(N_j)与f(N_i)的欧氏距离：||f(N_i)-f(N_j)||²；

(4b)将上述欧式距离代入如下公式，得到像素i的去噪后的灰度值v’(i)：

其中，归一化系数Z(i)由如下公式计算： 

2.根据权利要求1所述的图像去噪方法，其中步骤(1)所述的对输入图像逐像素取5×5的滑窗，计算其协方差矩阵 按如下公式计算：

其中，Ω为原始去噪图像中所有像素点的向量，|Ω|为去噪图像中的像素个数，y(N_i)为5×5像素大小的相似滑窗内的像素灰度值向量， 

为输入图像的像素灰度均值向量。

Images