CN102222327A

CN102222327A - 基于Treelet变换和最小均方误差估计的图像去噪方法

Info

Publication number: CN102222327A
Application number: CN 201110199552
Authority: CN
Inventors: 王桂婷; 焦李成; 张林刚; 公茂果; 钟桦; 张小华; 田小林; 侯彪; 王爽
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2011-07-15
Filing date: 2011-07-15
Publication date: 2011-10-19

Abstract

本发明公开了基于Treelet变换和最小均方误差估计的图像去噪方法，主要解决现有去噪方法对被高斯白噪声腐蚀的自然图像去噪效果不佳的问题。其实现步骤如下，(1)输入一幅待去噪图像；(2)搜索相似图像块；(3)Treelet变换获得基矩阵；(4)投影；(5)计算去噪后系数；(6)计算去噪后相似图像块矩阵；(7)计算相似度；(8)矩阵赋值；(9)判断是否处理完图像内所有图像块，若处理完，则进行步骤(10)，否则转入步骤(2)；(10)结果输出。本发明具有对含有高斯白噪声的自然图像去噪效果好的优点，能恢复出图像原有的特征，可用于图像分割、目标识别，变换检测等对图像的预处理。

Description

基于Treelet变换和最小均方误差估计的图像去噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种自然图像处理滤波技术领域的Treelet变换和最小均方误差估计的去噪方法。该方法可用于开展森林资源调查、灾害评估、城市规划、医学影像、天文学影像等领域的数字图像预处理。

背景技术

由于受各种条件的限制，图像在获取、编码、传输过程中会受到各种各样的噪声影响，这对图像分割，目标识别等其他一些后续处理工作带来了不利，因此对图像进行去噪很有必要。图像去噪解决了图像受到噪声干扰图像质量下降的问题，抑制噪声影响，提高了图像质量，是图像后处理的基础。

目前，人们更多采用频域滤波的方法来进行去噪图像。空域滤波是直接对图像的灰度进行处理，主要有全变差，中值滤波，双边滤波，非局部均值。空域滤波去噪方法主要是对图像进行模糊或者拟合操作来达到去噪的目的，导致的结果是图像中的高频信息(如边缘)会被滤除。

西安电子科技大学在其申请的专利“基于Treelet变换和非局部均值的图像去噪方法”(专利申请号201110001952.1，公开号CN102063708A)中公开了一种利用相似图像块灰度值Treelet变换的系数来计算计算图像块之间的相似度的非局部均值去噪方法。该方法可在高噪声情况下能提高去噪效果，但是该方法仍存在的不足是，对图像中细节信息保持不好，去噪后图像过平滑，导致处理后的图像结构信息(边缘、纹理、点)会被模糊或滤除。

Thaipanich，T.等人在“Improved image denoising with adaptive nonlocal means algorithm.IEEE Transactions on Consumer Electronics，vol.56，no.4，Nov.2010，2623-2630.”中提出了一种自适应的非局部均值去噪方法。该方法首先利用图像的梯度信息对图像进行分类，在属于同类的像素中搜索相似图像块，最后用加权平均方法进行去噪。该方法存在的不足是，直接使用含噪图像中的灰度计算图像块之间的相似度，获得的相似性权值不准确，使得最后去噪结果与原始图像存在偏差，另外在寻找相似快的时候对图像块进行了旋转操作，计算复杂度比较高。

发明内容

本发明的目的在于克服现有技术中的不足，提出一种基于Treelet变换和最小均方误差估计的图像去噪方法，在图像块周围搜索相似图像块，采用Treelet变换对相似图像块进行变换，再用最小均方误差估计对去噪后系数进行估计。

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

(1)输入一幅待去噪图像；

(2)选择相似图像块

2a)在待去噪图像中，任取一个像素点为中心，以5～11个像素点为边长，确定一个正方形的中心图像块，以21～41个像素为边长，确定一个正方形的搜索窗；

2b)在搜索窗中逐行扫描选取与中心图像块等大小的所有图像块；

2c)选择相似图像块

按照下式计算步骤2a)中的中心图像块和步骤2c)中选取的图像块的相似度：

d (y_{i}, y_{j}) = E {{| | y_{i} - y_{j} | |}_{2}^{2}}

其中，d(y_i，y_j)为中心图像块和选取的图像块的相似度，y_i表示2a)中中心图像块的灰度特征向量，y_j表示步骤2b)中选取的图像块的灰度特征向量，E表示数学期望，‖‖₂表示二范数；

按照下式计算与中心图像块相似的图像块中心像素位置集合：

S_{i} = {j | j &Element; X, d (y_{i}, y_{j}) < th + σ_{n}^{2}}

其中，S_i为相似图像块中心像素位置集合，j为相似图像块中心像素位置，X为输入的含噪图像，d(y_i，y_j)为中心图像块和选取的图像块的相似度，th为阈值，

为图像中的噪声方差；

如果步骤2b)搜索窗内图像块的中心像素位置在相似图像块中心像素位置集合S_i中，该图像块即为中心图像块的相似图像块；

2d)将步骤2a)中得到的中心图像块和步骤2c)中的相似图像块的像素灰度向量组成相似图像块矩阵；

(3)获得基矩阵；

3a)按照下式计算相似图像块矩阵Y去均值的矩阵

\overset{&OverBar;}{Y} = Y - E [Y]

其中，

为Y去均值的矩阵，Y为相似图像块矩阵，E[Y]为Y的行均值；

3b)对相似图像块去均值的矩阵进行Trelet变换，获得基矩阵；

(4)将步骤(2)得到的相似图像块去均值矩阵投影到步骤(3)中的基矩阵上；

(5)用最小均方误差估计方法计算去噪后系数；

5a)按照下式计算去噪后系数矩阵的方差：

{\hat{σ}}_{k}^{2} = \max [0, \frac{1}{p} Σ_{s = 1}^{p} c_{ks}^{2} - σ_{n}^{2}]

其中，

为去噪后系数矩阵

方差的第k行，max是取最大值操作，p表示系数矩阵

的行数，k表示系数矩阵的行位置，c_ki表示投影系数矩阵c的第k行第s列元素；

5b)按照下式计算收缩因子：

W_{k} = \frac{{\hat{σ}}_{k}^{2}}{{\hat{σ}}_{k}^{2} + σ_{n}^{2}}

其中，W_k为第k行的收缩因子，

为去噪后系数矩阵

方差的第k行，

为图像中的噪声方差，k∈{1，2，…，49}；

5c)计算去噪后系数矩阵：

\hat{c} = W \cdot c

其中，

为去噪后系数矩阵，W为收缩因子，c为投影系数矩阵；

(6)按照下式计算去噪后相似图像块矩阵

其中，B^T是基矩阵的转置，

是去噪后系数矩阵，E[Y]是Y的行均值；

(7)计算相似度

按下式计算去噪后中心图像块和相似图像块之间的相似度：

其中，wei为相似度，

和

分别为y_i和y_j去噪后的结果，‖‖₂表示二范数，h为惩罚因子；

(8)矩阵赋值：将两个矩阵中步骤2c)S_i记录的位置的元素分别赋值为相似度和相似度与去噪后相似图像块灰度的乘积；

(9)判断是否处理完图像内所有图像块，若处理完，则进行步骤(10)，否则转入步骤(2)，选择中心图像块时进行下采样操作，水平和垂直方向每Q个像素选择一个中心图像块；

(10)进行权值归一化处理，得到去噪结果。

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

第一，本发明采用Treelet变换对相似图像块进行处理，由于噪声对Treelet变换结果影响小，克服了现有技术受噪声影响较大的缺点，使得本发明的去噪结果对图像中的细节信息得到较为完整的保持。

第二，本发明计算去噪系数方法是一种无偏估计，克服了现有技术计算存在偏差的缺点，能够准确的计算去噪后系数和图像块之间的相似度，从而能够准确计算待估计像素点的恢复值，因此能够较好平滑噪声，同时保持和恢复图像的边缘和纹理细节，图像的去噪效果明显增强。

第三，本发明选择中心图像块时进行下采样操作，克服了现有技术对像素点处理存在的计算时间长的缺点，减少了运算时间，实现过程简单，且可以并行实现。

附图说明

图1为本发明的流程图；

图2为现有技术中基于KSVD的去噪方法和SA_DCT去噪方法的效果图；

图3为输入的噪声标准差为20的图像和本发明的去噪效果图。

具体实施方式

下面结合附图1对本发明的步骤作进一步的详细描述。

步骤1，输入一幅待去噪图像。

步骤2，搜索相似图像块。

2a)选择中心图像块和搜索窗。

对输入的含噪图像X，任取一个像素点为中心，以5～11个像素点为边长，取正方形区域为中心图像块。本发明实施例中选取大小为7×7的中心像素块。以21～41个像素为边长，确定一个正方形的搜索窗。本发明实施例中选取大小为39×39的搜索窗。

2b)选取图像块。在搜索窗中逐行扫描选取与中心图像块等大小的所有图像块。

2c)选择相似图像块

d (y_{i}, y_{j}) = E {{| | y_{i} - y_{j} | |}_{2}^{2}}

S_{i} = {j | j &Element; X, d (y_{i}, y_{j}) < th + σ_{n}^{2}}

其中，S_i为相似图像块中心像素位置集合，j为相似图像块中心像素位置，X为输入的含噪图像，d(y_i，y_j)为中心图像块和选择的图像块的相似度，th为阈值，本发明实施例中th取值为200，

为图像中的噪声方差；

2d)将步骤2a)中得到的中心图像块和步骤2c)中的相似图像块的像素灰度向量组成相似图像块矩阵。

步骤3，获得基矩阵。

3a)按照下式计算相似图像块矩阵Y去均值的矩阵

\overset{&OverBar;}{Y} = Y - E [Y]

其中，

为Y去均值的矩阵，Y为相似图像块矩阵，E[Y]为Y的行均值。

3b)对相似图像块去均值的矩阵进行Treelet变换，获得基矩阵。

b1)按照下式，计算灰度特征向量组成的矩阵Y的协方差矩阵

和相关系数矩阵

{\hat{Σ}}_{ij} = E [(y_{i} - E y_{i}) {(y_{j} - E y_{j})}^{T}]

{\hat{M}}_{ij} = {\hat{Σ}}_{ij} / \sqrt{{\hat{Σ}}_{ii} {\hat{Σ}}_{jj}}

其中，

表示协方差矩阵中的元素，E表示数学期望，y_i表示中心图像块的灰度特征向量，y_j表示与中心图像块相似的图像块的灰度特征向量，T表示矩阵转置操作，

表示相关系数矩阵中的元素，和

分别为y_i和y_j的协方差。

b2)初始化Treelet分解的层数l＝1，…，L，本发明实施例中选取分解的最高层数L＝48。当l＝0时，初始化相似度矩阵

为初始化基矩阵B₀为与协方差矩阵相同大小的单位矩阵；初始化和变量的下标集：δ＝{1，2，…，p}，其中p为矩阵Y的维数。

(b3)按照下式，找出相似度矩阵

中最相似的两个变量：

(α, β) = \arg \max {\hat{M}}^{(l - 1)}

其中，α和β分别表示矩阵

中最相似的两个变量的位置，arg max表示在矩阵中寻找最大值，是l-1层相似度矩阵。

(b4)对步骤3c)中得到的两个变量进行主成分变换，得到去相关后的第一主成分和第二主成分。这里，主成分变换是通过雅克比旋转实现。旋转角度θ_l由

且θ_l≤π/4计算得到，其中协方差矩阵

雅克比旋转矩阵J由下式计算：

其中，θ_l为选择角度。

按照下式更新基矩阵B_l和

B_l＝B_l-1J

{\hat{M}}^{(l)} = J^{T} {\hat{M}}^{(l - 1)} J

其中，B_l第l层级的基矩阵，B_l-1第l层级的基矩阵，J为雅克比旋转矩阵，J^T为雅克比旋转矩阵的转置，

为l层相似度矩阵。

(b5)旋转角度后使

即α对应的变量为主成分分析后的第一主成分，β对应的变量为第二主成分。尺度向量φ_l和细节向量ψ_l分别为此分解层级下的基矩阵B_l的第α和β列。

(b6)重复步骤(b2)至步骤(b5)直至l＝L层，得到第L层Trelet分解的基矩阵，记为B。

B＝[φ_l，ψ₁，…ψ_l]^T

其中，φ_l表示最高层的尺度，ψ_l表示第l层的细节向量，T表示矩阵转置操作。

步骤4，将步骤(2)得到的相似图像块去均值矩阵投影到步骤3b)中的基矩阵上。

按照下式将矩阵投影到基矩阵B上：

c = B \overset{&OverBar;}{Y}

c是相似图像块在基矩阵B上的投影系数，B是获得的基矩阵，

是相似图像块去均值矩阵。

步骤5，根据最小均方误差估计方法，计算去噪后系数矩阵；

5a)按照下式计算去噪后系数矩阵的方差：

{\hat{σ}}_{k}^{2} = \max [0, \frac{1}{p} Σ_{s = 1}^{p} c_{ks}^{2} - σ_{n}^{2}]

其中，

为去噪后系数矩阵

方差的第k行，max是取最大值操作，p表示系数矩阵

的行数，k表示系数矩阵

的行位置，c_ki表示投影系数矩阵c的第k行第s列元素；

5b)按照下式计算收缩因子：

W_{k} = \frac{{\hat{σ}}_{k}^{2}}{{\hat{σ}}_{k}^{2} + σ_{n}^{2}}

其中，W_k为第k行的收缩因子，

为去噪后系数矩阵

方差的第k行，

为图像中的噪声方差，k∈{1，2，…，49}；

5c)计算去噪后系数矩阵

\hat{c} = W \cdot c

其中，

为去噪后系数矩阵，W为收缩因子，c为投影系数矩阵，本发明实施例中收缩因子W取值范围为0-1；

步骤6，按照下式计算去噪后相似图像块矩阵

其中，B^T是基矩阵的转置，

是去噪后系数，E[Y]是Y的行均值；

步骤7，计算相似度。

计算经过去噪之后的相似图像块矩阵内相似图像块与中心图像块y_i的相似度wei，计算公式如下所示：

其中，

分别是y_i与y_j去噪后的结果，‖‖₂表示二范数，h是惩罚因子，本发明实施例中h取值为5σ_n，σ_n为图像中噪声的标准差。

步骤8，矩阵赋值：将两个矩阵中步骤2c)S_i记录的位置的元素分别赋值为相似度和相似度与去噪后相似图像块灰度的乘积。

步骤9，判断是否处理完图像内所有图像块，若处理完，则进行步骤(10)，否则转入步骤(2)，选择中心图像块时进行下采样操作，水平和垂直方向每Q个像素选择一个中心图像块，本发明实施例中Q取值为4。

步骤10，结果输出。

将两矩阵中相同位置的值做除法运算，进行权值归一化处理，得到去噪结果。

下面结合附图2、附图3对本发明的效果做进一步的描述。

附图2是现有技术基于KSVD的去噪方法和基于SA_DCT去噪方法取自文献中的效果图。

附图2(a)，附图2(b)和附图2(c)为现有技术基于KSVD的去噪方法和的效果图，去噪效果图来源于文献“Michael Elad，Michal Aharon，Image Denoising Via Sparse and Redundant Representations Over Learned Dictionaries，IEEE Trans.on Image Processing，Vol.15，No.12，Dec 2006”第8页。

附图2(d)，附图2(e)和附图2(f)为现有技术基于SA_DCT去噪方法，去噪效果图是来源于文献“Alessandro Foi，Vladimir Katkovnik，Pointwise Shape-Adaptive DCT for High-Quality Denoising and Deblocking of Grayscale and Color Images，IEEE Trans.on Image Processing，Vol.16，No.5，May 2007”第6页。

附图3是本发明去噪后获得的效果图。

附图3(a)，附图3(b)和附图3(c)是本发明输入的待去噪的图像，图像加了标准差为20的随机高斯白噪声，附图3(a)为Lena图像，附图3(b)为Barbara图像，附图3(c)为Cameraman图像，所有图像都是灰度图像，灰度级为256级，其中Lena和Barbara为512×512像素大小，Cameraman为256×256像素大小。附图3(d)，附图3(e)和附图3(f)是本发明的对附图3(a)，附图3(b)和附图3(c)三幅图像去噪效果图。

将本发明的去噪效果图3和现有技术去噪效果图2进行比较，对本发明的效果作进一步描述。

下表列出了附图3(a)，附图3(b)和附图3(c)三幅图像在本发明去噪下的结果和现有技术基于KSVD的方法的去噪结果、基于SA_DCT方法去噪结果的峰值信噪比PSNR对比。

从表中可以看出，Lena图像所采用的三种去噪方法中本发明的图像去噪结果的PSNR值比基于SA_DCT的方法的PSNR值高0.01dB，比基于KSVD的方法的PSNR值要高0.18dB，说明本发明的去噪结果的偏差最小，显示出了本发明对同质区域噪声抑制能力和保护点，线，边缘高频信息的能力，去噪结果图像中细节信息得到较为完整的保持。对其他两幅图像Barbara和Cameraman，本发明去噪结果的PSNR值比基于SA_DCT方法和基于KSVD方法两种方法都要高，而且三种方法图像去噪结果PSNR的平均值本发明也是最高的，从附图3(d)，附图3(e)和附图3(f)可以看出本发明的去噪结果细节信息保持完整，图像清晰，失真小，说明本发明与基于KSVD的去噪方法和SA_DCT方法相比，在对细节信息，像线和纹理进行恢复时表现出了更好的效果，对同质区域的平滑性能很好，而且对图像的边缘细节信息也保持的比较完整。