CN102663703B - 基于Treelet的Bayer型CFA图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于Treelet的Bayer型CFA图像去噪方法，主要解决现有去噪方法直接应用于CFA图像去噪易出现的噪声大量残余的问题。其实现步骤是：对输入的含噪CFA图像提取高频信息；对高频图像分块处理，逐个对每个图像块提取训练数据，并将训练数据投影到Treelet基矩阵中；通过对投影系数进行收缩，并通过逆投影获得去噪后的高频图像块；拼接所有去噪后的高频图像块，并加上输入CFA图像的低频信息得到去噪后的CFA图像。本发明能直接对CFA图像进行自适应去噪，同时能在很好地保持图像细节的情况下，减少去噪后图像中的噪声残余，可用于对单CCD或CMOS传感器相机捕获的CFA图像的降噪处理。

Description

基于Treelet的Bayer型CFA图像去噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及基于Treelet的Bayer型CFA图像去噪方法，可作为用于对单CCD或CMOS传感器相机捕获的CFA图像的降噪处理。

背景技术

数字彩色图像通常用红、绿、蓝三基色来表示颜色值，因而需要三个不同的二维矩阵来表示。出于成本考虑，目前大部分相机采用的是单个CCD或CMOS传感器，通过在传感器前加一个颜色滤波阵列CFA，只用一个矩阵来表示彩色图像。每个像素点上只有一个颜色值，而另外两个颜色值则根据其邻域信息来插值，这一插值技术被称为“去马赛克”技术。由于电磁效应与热效应对传感器的影响，数码相机捕获的图像通常会引入一定的噪声。噪声对去马赛克过程造成很大的影响，甚至会产生严重的虚假颜色效应。因此，对颜色滤波阵列CFA图像进行去噪处理很有必要，这对提高单传感相机捕获彩色图像的质量具有重要意义。

针对去噪与去马赛克过程的处理顺序的不同，有三种可能的策略：先去马赛克后去噪，先去噪后去马赛克，以及同时进行去马赛克去噪。第一种策略最大的优点在于现有的去噪技术一般都可以直接采用，而不需要附加任何处理。这种策略的缺点也是显而易见的，直接去马赛克使得噪声复杂化从而产生难以去除的虚假颜色效应，这使得后续去噪算法的设计极为困难。第二种策略是先对CFA图像进行去噪处理，然后在对降噪后的CFA图像进行去马赛克。此策略可以很好的解决第一种策略中出现虚假颜色效应的问题。鉴于CFA图像的特点，不能直接采用现有的灰度图像去噪技术，必须进行相应处理。针对Bayer型CFA图像，一般先将其拆分为四个子图像，然后分别对这四个子图像进行去噪，再组合成去噪后的CFA图像，最后进行颜色插值，得到最终的彩色图像。对子图像的处理可以直接采用现有的灰度图像的去噪算法。将CFA图像分割处理不但繁琐，而且并没有考虑到通道间颜色信息的相关性，往往会出现过平滑现象。第三种策略是去马赛克的同时进行降噪处理，即联合去马赛克去噪技术。这类技术在设计算法时需要考虑噪声的影响，虽然将去噪与去马赛克两个步骤合并在一起，但算法设计上更为复杂，如果设计不够好，不一定能提高处理效果和速度。

Zhang等人在文章“PCA-based spatially adaptive denoising of CFA images forsingle-sensor digital cameras”(IEEE Trans.Image Process.，vol.18，No.4，pp.797-812，2009)提出一种基于PCA的空域自适应CFA图像去噪方法，可直接对CFA图像数据进行去噪。该方法改进了单色图像中自适应PCA去噪方法，并将其推广到CFA图像去噪中，有效地降低了噪声对CFA图像去马赛克过程的影响。该方法针对Bayer型CFA图像的特点，首先对CFA图像进行高斯滤波，再通过最小均方误估计对图像高频信息的PCA变换系数进行处理，从而达到抑制噪声的目的。该方法虽然可直接应用于Bayer型CFA图像中，不需要进行分割处理，但不足之处是：对PCA系数处理之后仍然有较多噪声残余，尤其是在高噪声情况下，降噪效果不佳。

发明内容

本发明针对上述现有CFA图像去噪技术的缺陷，提出了一种基于Treelet的CFA图像去噪方法，以减少降噪后图像中的噪声残余，在高噪声下提高降噪效果。

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

(1)输入一幅待去噪的Bayer型CFA图像；

(2)提取CFA图像的高频信息：

2a)对输入的CFA图像进行二维高斯低通滤波，得到其低频图像；

2b)将输入CFA图像与其低频图像作差，得到输入CFA图像的高频图像；

(3)对输入CFA图像的高频图像逐像素取6×6像素大小的图像块；

(4)对提取的一个6×6像素大小的图像块进行Treelet去噪：

4a)以当前待去噪的6×6像素大小的图像块的中心为中心，在高频图像中取一个大小为34×34像素大小的窗口作为训练窗；

4b)在训练窗中，以第1行第1列的像素作为将要选取图像块的左上角顶点像素，取一个大小为6×6像素的图像块，然后分别沿水平方向和垂直方向每隔一个像素取一个大小为6×6像素的图像块，得到15×15总共225个图像块，构成了一个训练图像块集合C；

4c)从训练图像块集合C中选取与当前待去噪图像块均方误差最小的40个图像块，并拉成36×1的列向量，构建成大小为36×40训练数据

4d)利用雅克比旋转方法获取训练数据

的Treelet基矩阵B_l；

4e)将训练数据

投影到Treelet基矩阵B_l中，得到训练数据

的Treelet投影系数

4f)利用线性最小均方误差估计方法对投影系数进行收缩得到收缩后的系数

4g)将收缩后的系数

投影到Treelet基矩阵B_l的逆矩阵

中，得到去噪后的数据

完成当前待去噪图像块的去噪；

(5)重复执行步骤(4)完成对每个图像块的去噪，将所有去噪后的图像块按原来在高频图像中位置拼接起来，组成CFA图像去噪后的高频图像；

(6)将CFA图像去噪后的高频图像与步骤2a)中得到的低频图像相加，得到去噪后的CFA图像，并输出。

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

第一，本发明是针对Bayer型的CFA图像设计的一种自适应的去噪方法，能够直接对CFA图像进行去噪处理，而绝大部分现有的去噪技术则需要先将CFA图像分割成四个子图像，然后分别进行去噪处理，因此本发明实现起来更为简便；

第二，本发明通过在邻域中提取训练数据，在Treelet域中进行系数收缩，具有自适应性，当通道间噪声水平不一致时，现有技术的去噪效果差，但本发明仍然能获得很好的降噪效果，尤其在高噪声下，本发明优势更明显。

附图说明

图1为本发明的流程图；

图2为本发明使用的测试图像。

具体实施方式

下面结合附图，对本发明具体实现及效果作进一步的详细描述。

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

步骤1，输入一幅待去噪的Bayer型颜色滤波阵列CFA图像I_v。

步骤2，提取输入CFA图像I_v的高频信息：

2a)对输入的CFA图像I_v进行二维高斯低通滤波，得到其低频图像

首先，根据下面公式设计二维高斯低通滤波器G：

$G(x,y)=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\δ}}{e}^{-\frac{x^2+y^2}{2{\mathrm{\δ}}^2}},$

其中，G(x，y)为所需设计的滤波器中第x行第y列的元素，δ为平滑因子，且δ＝3；

然后，根据以下公式将滤波器与输入图像进行二维卷积，得到输入图像的低频图像

$I_v^l=I_v*G,$

其中，I_v为输入图像，G为高斯低通滤波器，*为二维卷积运算；

2b)将输入图像I_v与其低频图像

作差，得到输入图像的高频图像

即 $I_v^h=I_v-I_v^l.$

步骤3，对高频图像

从第15行第15列的像素开始作为正方形的左上角顶点取第一个图像块，然后分别沿水平方向和垂直方向每隔一个像素取一个6×6像素大小的图像块，直到高频图像

的倒数第20行倒数第20列的像素为止，将高频图像分割成了一系列6×6像素大小的图像块。

步骤4，对提取的一个6×6像素大小的图像块进行Treelet去噪：

4a)以当前待去噪的6×6像素大小的图像块的中心为中心，在高频图像

中取一个大小为34×34像素的窗口作为训练窗；

4b)在训练窗中，以第1行第1列的像素作为将要选取图像块的左上角顶点像素，取一个大小为6×6像素的图像块，然后分别沿水平方向和垂直方向每隔一个像素取一个大小为6×6像素的图像块，得到15×15总共225个图像块，构成了一个训练图像块集合C；

4c)计算训练图像块集合C中所有训练图像块与当前待去噪图像块的均方误差，将225个训练图像块按均方误差由小到大的顺序排列并选出前40个，将这40个训练图像块均拉成36×1的列向量，按上述由小到大的顺序排列组成一个36×40的二维矩阵X，该矩阵的第一列为当前待去噪的图像块，其余各列均为训练图像块，则所需提取的训练数据

4d)利用雅克比旋转方法获取训练数据

的Treelet基矩阵B_l，通过如下步骤进行：

4d1)计算训练数据矩阵

的协方差矩阵

和相关系数矩阵

$\widehat{\mathrm{\Σ}}=\left\{{\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}\right\},$

$\hat{M}=\left\{{\hat{M}}_{\mathrm{ij}}\right\},$

其中，

表示协方差矩阵中的元素，i和j分别表示矩阵的行和列，E表示数学期望，x_i表示训练数据矩阵的第i列，x_j表示训练数据矩阵的第j列，T表示矩阵转置操作，

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

为x_i和x_j的协方差，

和

分别为x_i和x_j各自的方差。

4d2)初始化Treelet分解的层数l＝1，初始化相似度矩阵

为协方差矩阵

初始化基矩阵为B₀，B₀为36×36大小的单位矩阵；

4d3)利用下式，找出相似度矩阵中最相似的两个变量：

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

其中，α和β分别表示矩阵中相似度最大的两个变量的位置索引，argmax表示在矩阵中寻找最大值元素的位置，是l-1层相似度矩阵；

4d4)对步骤4d3)中得到的两个变量进行雅克比旋转，旋转角度θ_l由以下方程组计算得到：

$\left\{\begin{array}{c}{J}^T{\widehat{\mathrm{\Σ}}}_{\mathrm{\α\β}}^{(l-1)}J=0\\ J^T{\widehat{\mathrm{\Σ}}}_{\mathrm{\β\α}}^{(l-1)}J=0\\ {\mathrm{\θ}}_l\≤\mathrm{\π}/4\end{array}\right.$

其中，J为雅克比旋转矩阵，矩阵形式如下：

其中，θ_l为旋转角度，

为第l-1层索引α对应的变量与索引β对应的变量的协方差，

为第l-1层索引β对应的变量与索引α对应的变量的协方差；

4d6)将基矩阵B_l和相似度矩阵

分别更新为：

B_l＝B_l-1J，

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

为l层相似度矩阵；

4d7)重复步骤4d3)至步骤4d6)直至l＝34层，得到第34层Treelet分解的基矩阵。

4e)将训练数据投影到Treelet基矩阵B_l中，得到训练数据

的Treelet投影系数： $\tilde{Y}=\tilde{X}\·B_l;$

4f)利用线性最小均方误差估计方法对投影系数

进行收缩，得到收缩后的系数

4f1)计算投影系数的协方差

${\mathrm{\Ω}}_{\tilde{Y}}=\frac{1}{40}\underset{i=1}{\overset{40}{\mathrm{\Σ}}}(y_i-{\stackrel{\‾}{y}}_i){(y_i-{\stackrel{\‾}{y}}_i)}^T,$

其中，y_i为

的第i列，

为

的第i列中所有元素的均值，∑为累加运算，T为矩阵转置运算；

4f2)根据协方差

对投影系数进行收缩，得到收缩后的系数

$\hat{Y}=\frac{{\mathrm{\Ω}}_{\tilde{Y}}-{\mathrm{\Ω}}_{V_y}}{{\mathrm{\Ω}}_{\tilde{Y}}}\tilde{Y},$

其中，

为噪声投影后的协方差，

的值等于噪声的方差，噪声方差为已知参数，实际中利用小波系数的绝对值中值法来估计，仿真时取值为5²，10²，15²，20²，25²，30²五组噪声方差之一。

4g)将收缩后的系数

投影到Treelet基矩阵B_l的逆矩阵

中，得到去噪后的数据

$\hat{X}=\hat{Y}\·B_l^T,$

其中，为第l层的基矩阵的转置，也是B_l的逆矩阵，

的第一列即为当前待去噪图像块降噪后的高频部分，将其拉回6×6图像块，即完成当前图像块的去噪。

步骤5，重复执行步骤(4)完成对每个图像块的去噪，将所有去噪后的图像块按原来在高频图像中位置拼接起来，组成CFA图像去噪后的高频图像

步骤6，将CFA图像去噪后的高频图像

与步骤2a)中得到的低频图像

相加，得到去噪后的CFA图像

即

并输出结果。

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

1.实验图像

实验所使用的测试图像为Kodak标准彩色图像库中的两幅图像，如图2所示，其中，图2(a)为Kodak标准彩色图像库中第6幅图像，图2(b)为Kodak标准彩色图像库中第18幅图像，两幅图像大小均为512×768像素，均具有红、绿、蓝三个颜色通道，每个颜色通道灰度级为256。

2.实验评价指标

评价颜色滤波阵列CFA图像去噪效果的方法是对一幅清晰干净的图像加入噪声后按Bayer型CFA格式进行下采样，从而获得含噪的Bayer型CFA图像，然后再进行去噪和去马赛克实验，得到去噪后的彩色图像，最后通过与原图对比进行评价。

CFA图像去噪克效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，如果结果图像噪声去除的很干净，且图像细节保持好，则说明该去噪方法效果好，反之效果较差。在客观上评价一幅图像的去马赛克效果好坏，本发明采用常用的彩色峰值信噪比CPSNR来衡量。

令真实彩色图像为I，图像大小为m×n×3，去噪后的彩色图像为

则彩色峰值信噪比按如下公式计算：

$\mathrm{CPSNR}=10{\log }_{10}\frac{{255}^2}{\frac{1}{3\mathrm{mn}}\underset{k=1}{\overset{3}{\mathrm{\Σ}}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(\hat{I}(i,j,k)-I(i,j,k))}^2},$

其中，k为通道序号索引，k＝1，2，3，i和j为图像中的像素位置索引，i＝1，2，...，m，j＝1，2，...，n，m和n分别为图像的行数和列数。

3.实验结果与分析

为了说明本发明的性能与优势，本发明与另外一种CFA图像去噪方法进行了对比实验。对比方法为Zhang等人在文章“PCA-based spatially adaptive denoising of CFAimages for single-sensor digital cameras”(IEEE Trans.Image Process.，vol.18，No.4，pp.797-812，2009)中提出一种基于PCA的空域自适应CFA图像去噪方法(记为PCA)。为了客观地进行实验对比，对去噪后的CFA图像都采用Zhang等人在文章“Colordemosaicking via directional linear minimum mean square-error estimation”(IEEE Trans.Image Process.，vol.14，no.12，pp.2167-2178，Dec.2005.)提出的去马赛克方法，得到最终的结果图，然后再进行比较。

试验中对Kodak彩色图像库中两幅图像进行了测试，实验评价指标采用彩色峰值信噪比CPSNR，实验结果如表1所示。

表1对两幅测试图像试验结果

表1中的粗体部分为最大值，从表1可见，当加入的噪声标准差为5时，基于PCA的方法对图2(a)结果的结果好，而对于图2(b)本发明的结果好；当加入的噪声标准差大于5时，本发明对两幅图像的测试结果都比基于PCA方法的结果好，并且噪声方差越大，本发明的优势越明显，尤其当噪声水平为30时，对于两幅图像本发明结果的CPSNR值要比基于PCA方法分别高出0.79dB和0.54dB。

Claims (3)

1.一种基于Treelet的Bayer型CFA图像去噪方法，包括如下步骤： 

（1）输入一幅待去噪的Bayer型CFA图像； 

（2）提取CFA图像的高频信息： 

2a）对输入的CFA图像进行二维高斯低通滤波，得到其低频图像； 

2b）将输入CFA图像与其低频图像作差，得到输入CFA图像的高频图像； 

（3）对输入CFA图像的高频图像逐像素取6×6像素大小的图像块； 

（4）对提取的一个6×6像素大小的图像块进行Treelet去噪： 

4a）以当前待去噪的6×6像素大小的图像块的中心为中心，在高频图像中取一个大小为34×34像素大小的窗口作为训练窗； 

4b）在训练窗中，以第1行第1列的像素作为将要选取图像块的左上角顶点像素，取一个大小为6×6像素的图像块，然后分别沿水平方向和垂直方向每隔一个像素取一个大小为6×6像素的图像块，得到15×15总共225个图像块，构成了一个训练图像块集合C； 

4c）从训练图像块集合C中选取与当前待去噪图像块均方误差最小的40个图像块，并拉成36×1的列向量，构建成大小为36×40训练数据

4d）利用雅克比旋转方法获取训练数据

的Treelet基矩阵B_l： 

4d1)计算训练数据矩阵的协方差矩阵

和相关系数矩阵

其中，

表示协方差矩阵中的元素，i和j分别表示矩阵的行和列，E表示数学期望，x_i表示训练数据矩阵

的第i列，x_j表示训练数据矩阵

的第j列，T表示矩阵转置操作，

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

为x_i和x_j的协方差，

和分别为x_i和x_j各自的方差； 

4d2)初始化Treelet分解的层数l=1，初始化相似度矩阵

为相关系数矩阵

 初始化基矩阵为B₀，B₀为36×36大小的单位矩阵； 

4d3)利用下式，找出第l-1层相似度矩阵

中最相似的两个变量： 

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

中相似度最大的两个变量的位置索引，arg max表示在矩阵中寻找最大值元素的位置，

是第l-1层相似度矩阵； 

4d4)对步骤4d3)中得到的两个变量进行雅克比旋转，旋转角度θ_l由以下方程组计算得到： 

其中，J为雅克比旋转矩阵，矩阵形式如下： 

其中，θ_l为旋转角度，为第l-1层索引α对应的变量与索引β对应的变量的协方差，

为第l-1层索引β对应的变量与索引α对应的变量的协方差； 

4d5）将基矩阵B_l和相似度矩阵分别更新为： 

B_l＝B_l-1J， 

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

为第l层相似度矩阵； 

4d6)重复步骤4d3）至步骤4d5）直至l=34层，得到第34层Treelet分解的基矩阵； 

4e）将训练数据

投影到Treelet基矩阵B_l中，得到训练数据

的Treelet投影系数

4f）利用线性最小均方误差估计方法对投影系数

进行收缩得到收缩后的系数

4g）将收缩后的系数

投影到Treelet基矩阵B_l的逆矩阵

中，得到去噪后的数据完成当前待去噪图像块的去噪； 

（5）重复执行步骤（4）完成对每个图像块的去噪，将所有去噪后的图像块按原来在高频图像中位置拼接起来，组成CFA图像去噪后的高频图像； 

（6）将CFA图像去噪后的高频图像与步骤2a）中得到的低频图像相加，得到去噪后的CFA图像，并输出。 

2.根据权利要求1所述的方法，其特征在于：步骤4e）中所述的将训练数据

投影到Treelet基矩阵B_l中，得到训练数据的Treelet投影系数

是由训练数据矩阵

与基矩阵相乘得到，即

3.根据权利要求1所述的方法，其特征在于：步骤4f）中所述的利用线性最小均方误差估计方法对投影系数

进行收缩，通过以下步骤进行： 

4f1）计算投影系数

的协方差

其中，y_i为的第i列，

为

的第i列中所有元素的均值，Σ为累加运算，T为矩阵转置运算； 

4f2）根据协方差

对投影系数

进行收缩，得到收缩后的系数

其中，

为噪声投影后的协方差，

的值等于噪声的方差，噪声方差为已知参数，实际中利用小波系数的绝对值中值法进行估计，仿真时取值为5²，10²，15²，20²，25²，30²六组噪声方差之一。 

Images