CN103400358A - 一种基于小波信息冗余的红外图像去噪算法 - Google Patents

Abstract

本发明公开了一种基于小波信息冗余的红外图像去噪算法，它通过在离散小波变换（DWT）过程中分别以不同的下采样方式获取多组含有相似冗余信息的小波系数，再利用噪声估计参数作为非线性变换函数阈值的基础，对高频小波系数进行非线性变换，抑制高频噪声，然后利用变换后小波系数重构（IDWT）多副图像，利用含相似冗余信息的多副图像加权进一步去除高频噪声，获取高质量红外图像。通过本发明进行红外图像去噪处理，既能保持原有图像的细节信息，又可以显著滤除图像高频噪声，提升图像清晰度，同时该方法可采用单片FPGA实现，有利于实现小型化、低成本、低功耗。

Description

一种基于小波信息冗余的红外图像去噪算法

技术领域

本发明属于红外图像信息处理技术领域，具体涉及一种基于小波信息冗余的红外图像去噪算法。

背景技术

红外成像技术在原理上具有良好的应用优势，在军事上可用于红外热制导、夜间作战等，在工业上可用于电力监测、夜间监控等，在医学上可用于非接触式体温测量等，因而红外成像技术成为现今研究的热门领域。

红外成像技术原理上是将被测对象的红外热辐射信号转换为二维可视化热像的红外图像处理仪器或设备，其输出信号具有灰度集中，对比度低等特点，从而红外图像增强成为必不缺少的手段，比如传统的直方图均衡法（HE）、对比度受限直方图均衡法、结合全局增强法与局部增强法的方式等。图像增强会导致图像噪声的放大，在对比度提升、细节分辨力增强的同时，噪声影响图像平滑度。为进一步提升图像质量，众多学者提出各种去噪算法，其中包括传统的均值滤波、中值滤波及高斯滤波等。Tomasi,C和Manduchi,R提出一种保留边沿的双边滤波方法，在去噪同时一定程度上保留图像细节信息，A,Buades等在双边滤波基础上增加一维信息，提出局部块相似信息滤波去噪方法，图像去噪效果更加明显。

小波变换在空域和频域同时具有良好的局部化特性，不仅可将图像的结构和纹理分别表现在不同分辨率层次上，而且具有检测边缘的能力，因而将小波变换应用于去噪时，能够在保留图像边缘细节特征的基础之上，又有效地去除噪声，提高图像的质量，使去噪后的图像更有利于人或机器的分析。小波去噪方式包括基于硬阈值、软阈值的离散小波变换（DWT）去噪。刘兴淼等提出一种基于小波变换的非线性红外图像增强算法，在其中采用软阈值方式进行小波去噪，张长江等利用离散平稳小波变换对红外图像噪声进行去噪增强，对原始图像进行小波变换获得低频和高频系数，接着根据低频系数的特点设计了非线性函数对低频系数进行增强，并对高频系数进行小波去噪，最后通过小波重构得到增强的图像，获得较为平滑清晰的红外图像。

上述方法，局部块相似信息滤波去噪方法去噪效果最佳，但对每个像素需进行块距离、点距离及灰度距离上的三维运算，计算量相当庞大，且随图像大小、块大小及块数目呈指数倍增长，影响其工程应用，双边滤波有学者对其进行简化加速进行实现后的形式去噪效果一般，基于小波变换的去噪算法大多针对噪声抑制曲线进行改进，且少有提到算法的可实现性。

发明内容

本发明的目的在于提供一种基于小波信息冗余的红外图像去噪算法，该方法可以有效解决相关红外图像预处理算法引起的噪声放大问题，在对背景噪声进行大幅度的抑制的同时，保留图像细节，增加图像清晰度，提升红外图像质量。

为了实现上述目的，本发明采用如下技术方案：

一种基于小波信息冗余的红外图像去噪算法，其特征在于，包括如下步骤：

（1）对红外图像进行离散小波分解，分解过程中分别以不同的下采样方式获取多组含有相似冗余信息的小波系数；

（2）对小波系数进行噪声估计；

（3）对高频小波系数进行非线性变换，以步骤（2）中噪声估计得到的噪声估计参数作为非线性变换函数阈值的基础；

（4）利用变换后的小波系数进行反变换重构多副图像，利用含相似冗余信息的多副图像加权进一步去除高频噪声。

进一步地，所述步骤（1）中的下采样方式有四种，分别为：取偶数行及偶数列方式、取奇数行及偶数列方式、取偶数行及奇数列方式，取奇数行及奇数列方式。

进一步地，所述步骤（2）中的噪声估计方法采用噪声方差估计方

进一步地，所述步骤（3）中，利用不完全贝塔函数对高频小波系数进行非线性变换，其中函数阈值选用噪声估计参数的倍数，倍数选择范围为3-6。

进一步地，所述步骤（4）中，利用含相似冗余信息的多副图像对齐通过均值进一步去除高频噪声。

进一步地，所述步骤（1）-（4）中均采用多层处理。

进一步地，采用三层处理。

进一步地，具体步骤包括：

（1）利用低频滤波器Lo_D及高频滤波器Hi_D对图像I(x,y)进行小波分解，通过不同的下采样方式获取四组含有相似信息的小波系数：分别以下标ee、oe、eo及oo表示；

A(x,y)=Lo_D*_v(Lo_D*_hI(x,y))  H(x,y)=Lo_D*_v(Hi_D*_hI(x,y))

V(x,y)=Hi_D*_v(Lo_D*_hI(x,y))  D(x,y)=Hi_D*_v(Hi_D*_hI(x,y))  （1）

其中*h表示行卷积，*v表示列卷积；

则一组小波系数为取偶数行及偶数列方式的下采样，即：

LL_ee=A([x/2]×2,[y/2]×2)  LH_ee=H([x/2]×2,[y/2]×2)

HL_ee=V([x/2]×2,[y/2]×2)  HH_ee=D([x/2]×2,[y/2]×2)    （2）

另取奇数行及偶数列、偶数行及奇数列，奇数行及奇数列可得到另三组小波系数：

LL_oe=A([x/2]×2+1,[y/2]×2)    LH_oe=H([x/2]×2+1,[y/2]×2)

HL_oe=V([x/2]×2+1,[y/2]×2)    HH_oe=D([x/2]×2+1,[y/2]×2)    （3）

LL_eo=A([x/2]×2,[y/2]×2+1)    LH_eo=H([x/2]×2,[y/2]×2+1)

HL_eo=V([x/2]×2,[y/2]×2+1)    HH_eo=D([x/2]×2,[y/2]×2+1)    （4）

LL_oo=A([x/2]×2+1,[y/2]×2+1)    LH_oo=H([x/2]×2+1,[y/2]×2+1)

HL_oo=V([x/2]×2+1,[y/2]×2+1)    HH_oo=D([x/2]×2+1,[y/2]×2+1)    （5）

（2）输入图像I中含有理想图像及噪声：

I=I_ideal+N    （6）

其中N为噪声，I_ideal为理想图像；

利用均值滤波f_average对I进行平滑近似I_ideal，然后从输入图像I中减去近似I_ideal，得到噪声估计；

$\tilde{N}=I-f_{\mathrm{average}}*I---\left(7\right)$

对

进行l层离散小波变换，获取多层小波系数，利用各层的HL、LH及HH小波系数方差衡量噪声等级；

${\mathrm{\δ}}_{\mathrm{HL}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HL}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HL}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{LH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{LH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{LH}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{HH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HH}}^l})}^2}---\left(8\right)$

其中为不同层的小波系数二维分布的高度及宽度；

（3）不完全贝塔函数如下式：

$\mathrm{IB}(x;a,b)=\frac{1}{B(a,b)}{\∫}_0^x{t}^{a-1}{(1-t)}^{b-1}\mathrm{dt}---\left(9\right)$

其中为归一化系数，a与b分别影响前段及后段曲线斜率；

利用不完全贝塔函数对小波系数进行非线性变换：

${\mathrm{HL}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HL}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HL}}^{l,r}\right)/T_{\mathrm{HL}}^l)*T_{\mathrm{HL}}^l& {\mathrm{HL}}^{l,r}\≤T_{\mathrm{HL}}^l\\ {\mathrm{HL}}^{l,r}& {\mathrm{HL}}^{l,r}>T_{\mathrm{HL}}^l\end{array}\right.$

${\mathrm{LH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{LH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{LH}}^{l,r}\right)/T_{\mathrm{LH}}^l)*T_{\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}\≤T_{\mathrm{LH}}^l\\ {\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}>T_{\mathrm{LH}}^l\end{array}\right.$

${\mathrm{HH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HH}}^{l,r}\right)/T_{\mathrm{HH}}^l)*T_{\mathrm{HH}}^l& {\mathrm{HH}}^{l,r}\≤T_{\mathrm{HH}}^l\\ {\mathrm{HH}}^{l,r}& {\mathrm{HH}}^{l,r}>T_{\mathrm{HH}}^l\end{array}\right.---\left(10\right)$

其中

分别为l层LH,HL,HH的压缩阈值，取

为各层LH,HL,HH估计噪声的倍数，即：

$T_{\mathrm{HL}}^l=g\×{\mathrm{\δ}}_{\mathrm{HL}}^l,T_{\mathrm{LH}}^l=g\×{\mathrm{\δ}}_{\mathrm{LH}}^l,T_{\mathrm{HH}}^l=g\×{\mathrm{\δ}}_{\mathrm{HH}}^l---\left(11\right)$

sign(·)表示取·的符号位，abs表示取绝对值，l=1,2,3表示小波分解层数，r=1,2,3,4分别表示ee,eo,oe,oo四组采样形式；

（4）对一层小波变换来说，四组经过压缩后的小波系数，通过反变换可得到四组具有像素位差的LLR图像，通过相应位移对齐后，对四副具有相似冗余信息的图像进行平均进一步去噪，通过三层计算最后得到即为最终的去噪结果图像；

${\mathrm{LL}}_{\mathrm{new}}^{l-1}=({\mathrm{LLR}}_{\mathrm{ee}}^{l-1}+{\mathrm{Shift}}_{\mathrm{eo}}\left({\mathrm{LLR}}_{\mathrm{eo}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oe}}\left({\mathrm{LLR}}_{\mathrm{oe}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oo}}\left({\mathrm{LLR}}_{\mathrm{oo}}^{l-1}\right))/4---\left(12\right)$

其中 ${\mathrm{LLR}}_{\mathrm{ee}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,1}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,1})$ 为偶数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{eo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,2}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,2})$ 为偶数行及奇数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oe}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,3}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,3})$ 为奇数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,4}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,4})$ 为奇数行及奇数列小波反变换结果；l=1,2,3代表三层小波，l=3时， ${\mathrm{LL}}_{\mathrm{new}}^{l,1}={\mathrm{LL}}^{l,1};{\mathrm{LL}}_{\mathrm{new}}^{l,2}={\mathrm{LL}}^{l,2};{\mathrm{LL}}_{\mathrm{new}}^{l,3}={\mathrm{LL}}^{l,3};{\mathrm{LL}}_{\mathrm{new}}^{l,4}={\mathrm{LL}}^{l,4};$ Shift_eo(·),Shift_oe(·),Shift_oo(·)分别表示对·在水平方向循环左移一个像素、在垂直方向循环上移一个像素及在对角线上向左上循环移一个像素。

由于采用上述技术方案，本发明具有以下有益效果：

（1）对图像进行小波变换，对图像进行多分辨率分解，对不同空间分辨率下的噪声进行滤除。

（2）在算法中加入噪声估计策略，可实现对系统自适应，根据系统动态调整小波系统非线性变换曲线，有效抑制噪声，同时增强细节，从而达到最优的去噪效果。

（3）利用离散小波变换（DWT）采样所丢弃的冗余信息对图像进一步去噪，提升图像质量。

（4）采用一片FPGA即可实现所有的计算步骤，达到实时处理的效果，减小了系统功耗，降低了实现成本，提高了系统的集成度。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1是本发明整体算法的框图；

图2是单层算法框图；

图3是非线性变换曲线；

图4是不同场景下含噪声红外图像（其中，4a为雨天场景，4b为室内场景，4c为晴天场景）；

图5是处理后的红外图像（其中，5a为雨天场景，5b为室内场景，5c为晴天场景,5d为5a中的I部细节示意图，5e为5b中的I部细节示意图，5f为5c中的I部细节示意图）。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

（1）利用低频滤波器Lo_D及高频滤波器Hi_D对图像I(x,y)进行小波分解，通过不同的下采样方式获取四组含有相似信息的小波系数：分别以下标ee、oe、eo及oo表示；

A(x,y)=Lo_D*_v(Lo_D*_hI(x,y))    H(x,y)=Lo_D*_v(Hi_D*_hI(x,y))

V(x,y)=Hi_D*_v(Lo_D*_hI(x,y))    D(x,y)=Hi_D*_v(Hi_D*_hI(x,y))    （1）

其中*h表示行卷积，*v表示列卷积；

则一组小波系数为取偶数行及偶数列方式的下采样，即：

LL_ee=A([x/2]×2,[y/2]×2)    LH_ee=H([x/2]×2,[y/2]×2)

HL_ee=V([x/2]×2,[y/2]×2)    HH_ee=D([x/2]×2,[y/2]×2)    （2）

另取奇数行及偶数列、偶数行及奇数列，奇数行及奇数列可得到另三组小波系数：

LL_oe=A([x/2]×2+1,[y/2]×2)    LH_oe=H([x/2]×2+1,[y/2]×2)

HL_oe=V([x/2]×2+1,[y/2]×2)    HH_oe=D([x/2]×2+1,[y/2]×2)    （3）

LL_eo=A([x/2]×2,[y/2]×2+1)    LH_eo=H([x/2]×2,[y/2]×2+1)

HL_eo=V([x/2]×2,[y/2]×2+1)    HH_eo=D([x/2]×2,[y/2]×2+1)    （4）

LL_oo=A([x/2]×2+1,[y/2]×2+1)    LH_oo=H([x/2]×2+1,[y/2]×2+1)

HL_oo=V([x/2]×2+1,[y/2]×2+1)    HH_oo=D([x/2]×2+1,[y/2]×2+1)    （5）

（2）对小波系数进行噪声估计。

图像噪声随成像环境而变化，例如不同的成像系统，其噪声会有差异，成像系统不同的工作环境同样影响噪声等级。去噪算法为达到自适应性，根据实时图像评估当前噪声等级，以此为依据提供高频小波系数非线性变换函数阈值，根据场景动态调整阈值，从而保证在不同的成像环境下，达到最优的去噪效果。

输入图像I中含有理想图像及噪声：

I=I_ideal+N    （6）

其中N为噪声，I_ideal为理想图像；

利用均值滤波f_average对I进行平滑近似I_ideal，然后从输入图像I中减去近似I_ideal，得到噪声估计；

$\tilde{N}=I-f_{\mathrm{average}}*I---\left(7\right)$

对进行l层离散小波变换，获取多层小波系数，利用各层的HL、LH及HH小波系数方差衡量噪声等级；

${\mathrm{\δ}}_{\mathrm{HL}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HL}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HL}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{LH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{LH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{LH}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{HH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HH}}^l})}^2}---\left(8\right)$

其中H^l,W^l为不同层的小波系数二维分布的高度及宽度；

（3）对各层小波系数进行非线性变换。

根据小波分析的原理，在高频小波系数中，噪声主要集中于幅值较小的小波系数，图像的细节信息主要集中于幅值较大的小波系数，并且绝对值较小的小波系数噪声成分较大，绝对值大的小波系数含噪声成分较少。针对这一特性，对绝对值小于阈值的小波系数进行压缩处理，保留绝对值大于阈值的小波系数，从而达到去噪目的。本发明采用不完全贝塔函数处理高频小波系数。不完全贝塔函数具有两个特性：1）经过原点及（1,1）点，且单调递增；2）原点及（1,1）点附近斜率小，中间斜率大。

$\mathrm{IB}(x;a,b)=\frac{1}{B(a,b)}{\∫}_0^x{t}^{a-1}{(1-t)}^{b-1}\mathrm{dt}---\left(9\right)$

其中

为归一化系数，a与b分别影响前段及后段曲线斜率；图3中为a=8，b=8时的曲线。

利用不完全贝塔函数对小波系数进行非线性变换：

${\mathrm{HL}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HL}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HL}}^{l,r}\right)/T_{\mathrm{HL}}^l)*T_{\mathrm{HL}}^l& {\mathrm{HL}}^{l,r}\≤T_{\mathrm{HL}}^l\\ {\mathrm{HL}}^{l,r}& {\mathrm{HL}}^{l,r}>T_{\mathrm{HL}}^l\end{array}\right.$

${\mathrm{LH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{LH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{LH}}^{l,r}\right)/T_{\mathrm{LH}}^l)*T_{\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}\≤T_{\mathrm{LH}}^l\\ {\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}>T_{\mathrm{LH}}^l\end{array}\right.$

${\mathrm{HH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HH}}^{l,r}\right)/T_{\mathrm{HH}}^l)*T_{\mathrm{HH}}^l& {\mathrm{HH}}^{l,r}\≤T_{\mathrm{HH}}^l\\ {\mathrm{HH}}^{l,r}& {\mathrm{HH}}^{l,r}>T_{\mathrm{HH}}^l\end{array}\right.---\left(10\right)$

其中

分别为l层LH,HL,HH的压缩阈值，取

为各层LH,HL,HH估计噪声的倍数，即：

$T_{\mathrm{HL}}^l=g\×{\mathrm{\δ}}_{\mathrm{HL}}^l,T_{\mathrm{LH}}^l=g\×{\mathrm{\δ}}_{\mathrm{LH}}^l,T_{\mathrm{HH}}^l=g\×{\mathrm{\δ}}_{\mathrm{HH}}^l---\left(11\right)$

sign(·)表示取·的符号位，abs表示取绝对值，l=1,2,3表示小波分解层数，r=1,2,3,4分别表示ee,eo,oe,oo四组采样形式；倍数选择范围为3-6。

（4）基于冗余信息的去噪。

对一层小波变换来说，四组经过压缩后的小波系数，通过反变换可得到四组具有像素位差的LLR图像，通过相应位移对齐后，对四副具有相似冗余信息的图像进行平均进一步去噪，通过三层计算最后得到

即为最终的去噪结果图像；

${\mathrm{LL}}_{\mathrm{new}}^{l-1}=({\mathrm{LLR}}_{\mathrm{ee}}^{l-1}+{\mathrm{Shift}}_{\mathrm{eo}}\left({\mathrm{LLR}}_{\mathrm{eo}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oe}}\left({\mathrm{LLR}}_{\mathrm{oe}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oo}}\left({\mathrm{LLR}}_{\mathrm{oo}}^{l-1}\right))/4---\left(12\right)$

其中 ${\mathrm{LLR}}_{\mathrm{ee}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,1}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,1})$ 为偶数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{eo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,2}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,2})$ 为偶数行及奇数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oe}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,3}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,3})$ 为奇数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,4}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,4})$ 为奇数行及奇数列小波反变换结果；l=1,2,3代表三层小波，l=3时， ${\mathrm{LL}}_{\mathrm{new}}^{l,1}={\mathrm{LL}}^{l,1};{\mathrm{LL}}_{\mathrm{new}}^{l,2}={\mathrm{LL}}^{l,2};{\mathrm{LL}}_{\mathrm{new}}^{l,3}={\mathrm{LL}}^{l,3};{\mathrm{LL}}_{\mathrm{new}}^{l,4}={\mathrm{LL}}^{l,4};$ Shift_eo(·),Shift_oe(·),Shift_oo(·)分别表示对·在水平方向循环左移一个像素、在垂直方向循环上移一个像素及在对角线上向左上循环移一个像素。

综上所述，本发明提出一种基于小波信息冗余的红外图像去噪算法，它通过在离散小波变换（DWT）过程中分别以不同的下采样方式获取多组含有相似冗余信息的小波系数，再利用噪声估计参数及特定函数对小波系数进行非线性变换，抑制高频噪声，然后利用变换后小波系数重构（IDWT）多副图像，利用含相似冗余信息的多副图像加权进一步去除高频噪声，获取高质量红外图像。此算法可在单片FPGA中进行实现，利用ALTERACYCLONIII芯片实现后的处理帧频达到50fps，满足实时性要求。

为说明算法环境适应性，分别选取雨天、室内及晴天下的不同的场景（依次为图4a、图4b及图4c）进行处理。图4是三张384×288原始噪声红外图像，从图中可以看出，图像经过预处理后，信息显示相对完整清晰，但噪声放大现象明显，尤其在信息量低的情况下，见图4a，背景噪声严重，而晴天下场景由于场景复杂信息量大，虽然背景噪声相对被弱化，在信息量小的区域（树枝丛中），图像噪声依然明显。通过图1中三层小波处理，每一层都采用图2的流程，可得到如图5（图5a、图5b及图5c）的处理后图像，另将局部细节图截取以作对比（如图5d、图5e及图5f）。从图5整体效果可以看出，图像平滑效果明显，背景噪声被大幅度抑制，而同时从所截取的细节图可以看出图像细节更加清晰，具有更好的视觉效果。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于小波信息冗余的红外图像去噪算法，其特征在于，包括如下步骤：

（1）对红外图像进行离散小波分解，分解过程中分别以不同的下采样方式获取多组含有相似冗余信息的小波系数；

（2）对小波系数进行噪声估计；

（3）对高频小波系数进行非线性变换，以步骤（2）中噪声估计得到的噪声估计参数作为非线性变换函数阈值的基础；

（4）利用变换后的小波系数进行反变换重构多副图像，利用含相似冗余信息的多副图像加权进一步去除高频噪声。

2.根据权利要求1所述的红外图像去噪算法，其特征在于，所述步骤（1）中的下采样方式有四种，分别为：取偶数行及偶数列方式、取奇数行及偶数列方式、取偶数行及奇数列方式，取奇数行及奇数列方式。

3.根据权利要求1所述的红外图像去噪算法，其特征在于，所述步骤（2）中的噪声估计方法采用噪声方差估计方法。

4.根据权利要求1所述的红外图像去噪算法，其特征在于，所述步骤（3）中，利用不完全贝塔函数对高频小波系数进行非线性变换，其中函数阈值选用噪声估计参数的倍数，倍数选择范围为3-6。

5.根据权利要求1所述的红外图像去噪算法，其特征在于，所述步骤（4）中，利用含相似冗余信息的多副图像对齐通过均值进一步去除高频噪声。

6.根据权利要求1所述的红外图像去噪算法，其特征在于，所述步骤（1）-（4）中均采用多层处理。

7.根据权利要求6所述的红外图像去噪算法，其特征在于，采用三层处理。

8.根据权利要求1所述的红外图像去噪算法，其特征在于，具体步骤包括：

（1）利用低频滤波器Lo_D及高频滤波器Hi_D对图像I(x,y)进行小波分解，通过不同的下采样方式获取四组含有相似信息的小波系数：分别以下标ee、oe、eo及oo表示；

A(x,y)=Lo_D*_v(Lo_D*_hI(x,y))  H(x,y)=Lo_D*_v(Hi_D*_hI(x,y))

V(x,y)=Hi_D*_v(Lo_D*_hI(x,y))  D(x,y)=Hi_D*_v(Hi_D*_hI(x,y))  （1）

其中*h表示行卷积，*v表示列卷积；

则一组小波系数为取偶数行及偶数列方式的下采样，即：

LL_ee=A([x/2]×2,[y/2]×2)  LH_ee=H([x/2]×2,[y/2]×2)

HL_ee=V([x/2]×2,[y/2]×2)  HH_ee=D([x/2]×2,[y/2]×2)    （2）

另取奇数行及偶数列、偶数行及奇数列，奇数行及奇数列可得到另三组小波系数：

LL_oe=A([x/2]×2+1,[y/2]×2)  LH_oe=H([x/2]×2+1,[y/2]×2)

HL_oe=V([x/2]×2+1,[y/2]×2)  HH_oe=D([x/2]×2+1,[y/2]×2)    （3）

LL_eo=A([x/2]×2,[y/2]×2+1)  LH_eo=H([x/2]×2,[y/2]×2+1)

HL_eo=V([x/2]×2,[y/2]×2+1)  HH_eo=D([x/2]×2,[y/2]×2+1)    （4）

LL_oo=A([x/2]×2+1,[y/2]×2+1)  LH_oo=H([x/2]×2+1,[y/2]×2+1)

HL_oo=V([x/2]×2+1,[y/2]×2+1)  HH_oo=D([x/2]×2+1,[y/2]×2+1)    （5）

（2）输入图像I中含有理想图像及噪声：

I=I_ideal+N    （6）

其中N为噪声，I_ideal为理想图像；

利用均值滤波f_average对I进行平滑近似I_ideal，然后从输入图像I中减去近似I_ideal，得到噪声估计；

$\tilde{N}=I-f_{\mathrm{average}}*I---\left(7\right)$

对进行l层离散小波变换，获取多层小波系数，利用各层的HL、LH及HH小波系数方差衡量噪声等级；

${\mathrm{\δ}}_{\mathrm{HL}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HL}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HL}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{LH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{LH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{LH}}^l})}^2}$

${\mathrm{\δ}}_{\mathrm{HH}}^l=\sqrt{\underset{i,j}{\overset{H^l,W^l}{\mathrm{\Σ}}}\frac{1}{H^l*W^l}{({\mathrm{HH}}_{i,j}^l-\stackrel{\‾}{{\mathrm{HH}}^l})}^2}---\left(8\right)$

其中H^l,W^l为不同层的小波系数二维分布的高度及宽度；

（3）不完全贝塔函数如下式：

$\mathrm{IB}(x;a,b)=\frac{1}{B(a,b)}{\∫}_0^x{t}^{a-1}{(1-t)}^{b-1}\mathrm{dt}---\left(9\right)$

其中

为归一化系数，a与b分别影响前段及后段曲线斜率；

利用不完全贝塔函数对小波系数进行非线性变换：

${\mathrm{HL}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HL}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HL}}^{l,r}\right)/T_{\mathrm{HL}}^l)*T_{\mathrm{HL}}^l& {\mathrm{HL}}^{l,r}\≤T_{\mathrm{HL}}^l\\ {\mathrm{HL}}^{l,r}& {\mathrm{HL}}^{l,r}>T_{\mathrm{HL}}^l\end{array}\right.$

${\mathrm{LH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{LH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{LH}}^{l,r}\right)/T_{\mathrm{LH}}^l)*T_{\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}\≤T_{\mathrm{LH}}^l\\ {\mathrm{LH}}^l& {\mathrm{LH}}^{l,r}>T_{\mathrm{LH}}^l\end{array}\right.$

${\mathrm{HH}}_{\mathrm{shrink}}^{l,r}=\left\{\begin{array}{cc}\mathrm{sign}\left({\mathrm{HH}}^{l,r}\right)*\mathrm{IB}(\mathrm{abs}\left({\mathrm{HH}}^{l,r}\right)/T_{\mathrm{HH}}^l)*T_{\mathrm{HH}}^l& {\mathrm{HH}}^{l,r}\≤T_{\mathrm{HH}}^l\\ {\mathrm{HH}}^{l,r}& {\mathrm{HH}}^{l,r}>T_{\mathrm{HH}}^l\end{array}\right.---\left(10\right)$

其中

分别为l层LH,HL,HH的压缩阈值，取

为各层LH,HL,HH估计噪声的倍数，即：

$T_{\mathrm{HL}}^l=g\×{\mathrm{\δ}}_{\mathrm{HL}}^l,T_{\mathrm{LH}}^l=g\×{\mathrm{\δ}}_{\mathrm{LH}}^l,T_{\mathrm{HH}}^l=g\×{\mathrm{\δ}}_{\mathrm{HH}}^l---\left(11\right)$

sign(·)表示取·的符号位，abs表示取绝对值，l=1,2,3表示小波分解层数，r=1,2,3,4分别表示ee,eo,oe,oo四组采样形式；

（4）对一层小波变换来说，四组经过压缩后的小波系数，通过反变换可得到四组具有像素位差的LLR图像，通过相应位移对齐后，对四副具有相似冗余信息的图像进行平均进一步去噪，通过三层计算最后得到

即为最终的去噪结果图像；

${\mathrm{LL}}_{\mathrm{new}}^{l-1}=({\mathrm{LLR}}_{\mathrm{ee}}^{l-1}+{\mathrm{Shift}}_{\mathrm{eo}}\left({\mathrm{LLR}}_{\mathrm{eo}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oe}}\left({\mathrm{LLR}}_{\mathrm{oe}}^{l-1}\right)+{\mathrm{Shift}}_{\mathrm{oo}}\left({\mathrm{LLR}}_{\mathrm{oo}}^{l-1}\right))/4---\left(12\right)$

其中 ${\mathrm{LLR}}_{\mathrm{ee}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,1}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,1}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,1})$ 为偶数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{eo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,2}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,2}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,2})$ 为偶数行及奇数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oe}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,3}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,3}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,3})$ 为奇数行及偶数列小波反变换结果； ${\mathrm{LLR}}_{\mathrm{oo}}^{l-1}=\mathrm{IDWT}({\mathrm{LL}}_{\mathrm{new}}^{l,4}+{\mathrm{LH}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HL}}_{\mathrm{shrink}}^{l,4}+{\mathrm{HH}}_{\mathrm{shrink}}^{l,4})$ 为奇数行及奇数列小波反变换结果；l=1,2,3代表三层小波，l=3时， ${\mathrm{LL}}_{\mathrm{new}}^{l,1}={\mathrm{LL}}^{l,1};{\mathrm{LL}}_{\mathrm{new}}^{l,2}={\mathrm{LL}}^{l,2};{\mathrm{LL}}_{\mathrm{new}}^{l,3}={\mathrm{LL}}^{l,3};{\mathrm{LL}}_{\mathrm{new}}^{l,4}={\mathrm{LL}}^{l,4};$ Shift_eo(·),Shift_oe(·),Shift_oo(·)分别表示对·在水平方向循环左移一个像素、在垂直方向循环上移一个像素及在对角线上向左上循环移一个像素。

9.根据权利要求1-8任一项所述的红外图像去噪算法，其特征在于，所述算法在单片FPGA中进行实现。

Images