CN103606132A - 基于空域和时域联合滤波的多帧数字图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于空域和时域联合滤波的多帧数字图像去噪方法。其方案是：输入低光照环境下采集的同一场景的多帧数字图像；选取多帧图像中最清晰的一张作为基准图像；对多帧数字图像进行全局和局部配准；以基准图像为参考，在空域和时间域中建立与基准图像中局部区域相似的相似组，并采用相似组对基准图像中的局部区域进行协同滤波；利用滤波后图像中的亮度信息以及像素分布信息，对滤波后图像依次进行颜色校正和对比度增强，得到去噪后的图像。本发明能有效抑制数字图像中的噪声和运动模糊等因素对于图像质量的影响，并在抑制噪声的同时有效的保留图像中的细节纹理，可用于提高低光照环境下数字图像采集设备中的图像质量。

Description

基于空域和时域联合滤波的多帧数字图像去噪方法

技术领域

本发明属于图像处理技术领域，特别是涉及图像去噪方法，可用于消费类数字图像采集设备、医学影像、天文影像以及环境变化评估等领域的数字图像预处理。

背景技术

图像去噪是图像处理领域的一个热点问题，也是一个极具挑战的研究方向。图像中的噪声会阻碍人们对图像内容的理解，采用图像去噪可以有效的抑制图像中噪声对图像质量的影响，提高人们对图像内容的认知程度，以便对图像作进一步的处理。

根据数字图像的特点及其统计特性，多年来已经有很多学者提出了多种不同的去噪算法，按照其实现方法大致可以分为空域和频域两类，且其中大多都是基于图像局部信息的平滑处理。基于图像局部信息的平滑处理会使图像丢失很多细节信息，去噪效果十分不理想。2005年A.Buades，B.Coll等人对双边滤波去噪算法进行了改进，提出了一种非局部均值滤波的去噪方法。该方法打破了以往采用的图像空间域中“局部平滑”的思想，转而利用自然图像中广泛存在的“非局部”的空域相似性，通过在整幅图像中搜索相似块进行加权平均，以达到抑制噪声的效果。2007年，K.Dabov，Alessandro Foi等人在非局部均值图像去噪方法的基础上，提出一种块匹配三维协同滤波的去噪方法。该方法结合了图像空域的非局部相似性和频域中的稀疏性，在获取空域相似块的基础上，对多个空域相似块构成的分组在频域中进行三维协同滤波，可以在有效抑制图像噪声的同时尽可能的保留图像中的细节信息，是目前公认的性能较好的图像去噪方法。然而，对于单帧的图像去噪问题而言，在未知噪声统计特性的前提下从噪声图像中恢复出原始的无噪声目标图像，必须在一定的约束条件之下才可以获得近似的最优解。已有的方法大都将对图像中噪声的统计特性进行理想的先验性假设作为求解近似最优解的约束条件，由于与实际数字图像采集设备所获取的自然场景图像中的噪声统计特性存在较大差异，因而去噪效果并不理想。

相比于单帧图像，使用针对同一场景的多帧图像可以提供更为丰富的时域信息，进而为图像去噪研究提供了新的思路。对于低光照环境下基于多帧的图像去噪，为了避免由于长时间曝光所导致的各帧图像中的运动模糊，要求所采集的针对同一场景的多帧图像必须是在较短的曝光时间内获取的，而较短的曝光时间势必会由于进光量减少导致图像中的噪声强度显著增强。这种情况下，如何有效的利用多帧图像中的空域和时域相似性来抑制噪声的影响，是一项非常复杂而又富有挑战性的工作。此外，由于拍摄者手部的抖动所造成的多帧图像之间存在相对位移，以及场景中的局部运动所导致的图像局部模糊，都为使用多帧图像去噪设置了障碍。直接将已有的用于单帧的图像去噪方法扩展到多帧图像去噪，由于没有充分利用多帧图像所提供的丰富的时域信息，同时不能很好的解决上述的多帧图像之间的相对位移和场景中局部运动对图像质量造成的影响，因而很难获得理想的去噪效果。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于空域和时域联合滤波的多帧数字图像去噪方法，在抑制噪声的同时消除低光照环境下多帧图像之间的相对位移和场景中局部运动所导致的图像模糊，提高图像的去噪效果，增强图像的主观质量。

本发明目的技术方案是：通过对低光照环境下采集的多帧图像进行全局和局部配准，分析和判别场景中的局部运动区域，利用局部运动区域的清晰度测度和多帧图像的时域相关性消除运动模糊；结合多帧图像的空域和时域相似性对图像进行空域和时域的融合滤波，抑制图像中的噪声；结合场景的亮度特性对去噪图像进行颜色校正和进行对比度增强，实现多帧图像去噪。其实现步骤包括如下：

（1）输入N帧低光照环境下采集的针对同一场景的自然图像z_i,i=1,2,…,N；

（2）计算各帧图像的平均水平梯度

和平均垂直梯度将与之和作为图像清晰度测度G_i，选取G_i最大的一帧，作为参考图像z_r,1≤r≤N，并将该参考图像的清晰度测度记为G_max；

（3）计算各帧图像相对于参考图像z_r的全局运动矢量V_i；

（4）计算各帧图像清晰度测度G_i与参考图像清晰度测度G_max的比值R_i，将R_i大于阈值0.99所对应的第i帧图像记为候选参考图像z_r′，计算各帧图像相对于该候选参考图像z_r′的候选全局运动矢量

（5）计算全局运动矢量V_i和候选全局运动矢量的测度，并依据该测度从参考图像z_r和候选参考图像z_r′中选取基准图像

并将该基准图像

对应的全局运动矢量作为基准全局运动矢量

（6）将输入的各帧图像都划分成大小相同的J个搜索块，利用上述基准全局运动矢量计算各帧图像中每个搜索块相对于基准图像

中对应搜索块的局部运动矢量V_{i_j},1≤i≤N且i≠s,1≤j≤J；

（7）将各帧图像中的每个搜索块都划分为大小相同的Q个融合块，使用步骤（6）中所得局部运动矢量V_{i_j}，在基准图像

中找出与各帧图像中的每个融合块所对应的基准融合块，选取各帧图像中与同一基准融合块对应的融合块构成时域融合块组；

（8）分别计算每个时域融合块组中每个融合块与其对应基准融合块的像素差值，并将这些像素差值的平均值作为该时域融合块组中每个融合块与其对应基准融合块的时域相似性测度δ_{i_j_k},1≤i≤N且i≠s,1≤j≤J,1≤k≤Q；

（9）依据步骤（8）中所得的时域相似性测度δ_{i_j_k}，计算时域滤波权值ω_{i_j_k}：

${\mathrm{\&omega;}}_{i\_j\_k}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;16\\ (48-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/32,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;48\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>48\end{array}\right.;$

（10）用所述时域滤波权值ω_{i_j_k}对每个时域融合块组所对应的基准融合块进行时域加权融合滤波，得到时域融合滤波图像z_fu；

（11）依据步骤（8）中所得各个时域融合块组中每个融合块与其对应基准融合块的像素差值，统计像素差值大于阈值24的像素个数占融合块总像素个数的百分比，计算这些百分比的标准差σ_j,k，若σ_j,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域；

（12）对步骤（11）中所判定的局部运动区域，依据时域相似性测度δ_{i_j_k}，按照下式计算局部运动区域的时域滤波权值ω′_{i_j_k}:

${\mathrm{\&omega;}}_{i\_j\_k}^{\&prime;}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;6\\ (12-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/6,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;12\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>12\end{array}\right.;$

（13）用局部运动区域的时域滤波权值ω′_{i_j_k}对所判定的局部运动区域所对应的基准融合块进行加权融合滤波，使用该融合结果覆盖时域融合滤波图像z_fu中对应位置的融合结果，得到最终的时域融合图像并对该时域融合图像

进行单帧空域非局部均值滤波，得到空域滤波图像z_f；

（14）计算空域滤波图像z_f的平均亮度

根据该平均亮度对空域滤波图像z_f进行亮度校正和色度校正，得到亮度和色度校正图像z_c，对z_c再进行局部对比度增强，得到最终的去噪结果图像z_out。

本发明与现有方法相比具有以下优点：

1.本发明由于在对多帧图像进行全局和局部配准的过程中，依据图像清晰度测度和全局运动矢量测度进行基准图像的选取，并使用该基准图像作为参考进行全局和局部配准，提高了配准的精度，为后续时域融合滤波提供了准确的时域相关信息，有助于提升图像去噪的效果；

2.本发明由于在对多帧图像进行时域融合的过程中，对场景中的局部运动区域进行了判定，并对局部运动区域采用了与场景中其他区域不同的时域融合方法，可在保证时域去噪效果的同时消除场景中由于局部运动所引起的图像局部模糊，提高了去噪图像的主观质量；

3.本发明由于对空域滤波图像进行了亮度和颜色校正，降低了光照不足所引起的图像亮度和色度偏差，并对亮度和色度校正后图像进行局部对比度增强，消除了空时域滤波所导致的对比度下降，进一步提升了去噪图像的主观质量。

附图说明

图1是本发明的流程图；

图2是本发明中使用的图像下采样示意图；

图3是本发明中进行中空域非局部均值滤波权值搜索示意图；

图4是本发明测试使用的低光照环境下采集的6帧自然图像；

图5是现有非局部均值滤波去噪算法的结果图像；

图6是现有块匹配三维协同滤波去噪算法的结果图像；

图7是本发明仿真得到的去噪结果图像；

图8是输入的第1帧自然图像的局部原始分辨率显示。

图9是现有非局部均值滤波去噪算法的结果图像的局部原始分辨率显示。

图10是现有块匹配三维协同滤波去噪算法的结果图像的局部原始分辨率显示。

图11是本发明仿真所得去噪结果图像的局部原始分辨率显示。

具体实施方式

下面结合附图对本发明的实施例做详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和过程，但本发明的保护范围不限于下述的实施例。

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

步骤1：输入多帧图像。

（1a）输入低光照环境下采集的针对同一场景的N帧自然图像z_i,i=1,2,…,N，本实施例中取N=6，但不局限于N=6的情况；

（1b）输入各帧图像的分辨率为W×H像素，本实施例中W=3264，H=2448，但不局限于输入图像分辨率为3264×2448的情况；

（1c）本实施例中，输入图像的格式为YUV，对于其他图像格式，例如RGB、YCbCr、HSL等均可转化为YUV格式后进行输入。

步骤2：选取参考图像。

（2a）对步骤1中输入的6帧图像z_i,i=1,2,…,6的亮度分量Y进行2级下采样：

图2给出本步骤进行下采样的原理示意图，图2中a,b,c,d表示待采样图像中大小为2×2的像素块，x表示下采样图像中由a,b,c,d下采样所得到的像素，x的计算公式如下：

$x=\frac{(a+b+c+d+2)}{4}---1),$

依照图2，第1级下采样过程如下：

输入图像z_i的亮度分量Y，对该亮度分量Y中每个大小为2×2的像素块进行下采样，

使用式1）计算该像素块的下采样值，所有2×2像素块的下采样值构成该亮度分量Y的第1级下采样图像z_{1_i}；

第2级的下采样过程如下：

输入图像z_i亮度分量Y的第1级下采样图像z_{1_i}，对该下采样图像z_{1_i}中每个大小为2×2的像素块使用式1）计算该像素块的下采样值，所有2×2像素块的下采样值构成输入图像z_i亮度分量Y的第2级下采样图像z_{2_i}；

所述对输入图像z_i的亮度分量Y的下采样不限于本实施例中所使用的方法，也可以采用已有的其他下采样方法，例如从每2×2个像素块中选择1个像素作为下采样图像中对应位置的像素；

（2b）分别计算下采样图像z_{2_i}中每个像素的水平梯度和垂直梯度，以像素坐标为(x,y)的像素为例，其水平梯度G_{i_h}(x,y)和垂直梯度G_{i_v}(x,y)的计算公式如下：

$\left\{\begin{array}{c}{G}_{i\_h}(x,y)=S_h*A(x,y)\\ G_{i\_v}(x,y)=S_v*A(x,y)\end{array}\right.,$

其中，A(x,y)表示下采样图像z_{2_i}中以像素坐标(x,y)为中心大小为3×3的图像块，S_h和S_v分别表示3×3大小的水平Sobel算子S_h和垂直Sobel算子S_v：

$S_h=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right],S_v=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right],$

“*”表示卷积运算；

所述计算水平梯度和垂直梯度的Sobel算子不限于本实施例中所使用的Sobel算子，也可以使用其他形式的Sobel算子；

（2c）分别计算各帧的下采样图像z_{2_i}的平均水平梯度

和垂直梯度

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{G}}_{i\_h}=\frac{1}{W_{2\_i}\&CenterDot;H_{2\_i}}{\mathrm{\&Sigma;}}_{x=1}^{W_{2\_i}}{\mathrm{\&Sigma;}}_{y=1}^{H_{2\_i}}{G}_{i\_h}(x,y)\\ {\stackrel{\&OverBar;}{G}}_{i\_v}=\frac{1}{W_{2\_i}\&CenterDot;H_{2\_i}}{\mathrm{\&Sigma;}}_{x=1}^{W_{2\_i}}{\mathrm{\&Sigma;}}_{y=1}^{H_{2\_i}}{G}_{i\_v}(x,y)\end{array}\right.,$

其中，W_{2_i}和H_{2_i}分别表示下采样图像z_{2_i}的宽度和高度；

（2d）分别计算各帧的下采样图像z_{2_i}的图像清晰度测度G_i：

$G_i={\stackrel{\&OverBar;}{G}}_{i\_h}+{\stackrel{\&OverBar;}{G}}_{i\_v},$

选取下采样图像z_{2_i}中G_i最大的一帧，选取其对应的输入图像z_i作为参考图像z_r,1≤r≤N，并记该参考图像z_r的清晰度测度为G_max；

所述计算图像清晰度测度G_i不限于本实施例中所使用平均水平梯度

和垂直梯度

之和，也可以使用其他形式的图像清晰度测度G_i，例如使用平均水平梯度和垂直梯度

的平方和等。

步骤3：计算全局运动矢量V_i。

（3a）对步骤1中所述的输入图像z_i,i=1,2,…,6的亮度分量Y，利用与步骤（2a）相同的方法进行M=4级下采样，得到输入图像z_i亮度分量Y的下采样图像z_{m_i}，以及参考图像z_r亮度分量Y的下采样图像z_{m_r},1≤m≤4,1≤r≤6；

所述的下采样级数M不限于本实施例中所选取的M=4，M的取值与输入图像的大小有关，可依照第M级下采样图像的分辨率不大于256×256像素确定下采样级数M；

（3b）以第4级参考图像的下采样图像z_{4_r}为参考，使用第4级下采样图像z_{4_i}相对于该参考的中心位置进行水平和垂直位移，并计算位移后下采样图像z_{4_i}与参考图像的下采样图像z_{4_r}的均方误差D_{4_i}：

$D_{4\_i}=\frac{1}{W_{4\_i}\&CenterDot;H_{4\_i}}{\mathrm{\&Sigma;}}_{x=-W_{4\_i}}^{W_{4\_i}}{\mathrm{\&Sigma;}}_{y=-H_{4\_i}}^{H_{4\_i}}(p(x,y)-q(x+u_{4\_i},y+v_{4\_i}){)}^2,$

-W_{4_i}≤u_{4_i}≤W_{4_i}，-H_{4_i}≤v_{4_i}≤H_{4_i}，

其中，W_{4_i}和H_{4_i}分别表示下采样图像z_{4_i}的宽度与高度，p(x,y)表示参考图像的下采样图像z_{4_r}中坐标为(x,y)的像素，q(x,y)表示下采样图像z_{4_i}中坐标为(x,y)的像素，u_{4_i}表示水平位移，v_{4_i}表示垂直位移；

（3c）选取使均方误差D_{4_i}最小的水平位移u_{4_i}和垂直位移v_{4_i}，分别作为全局运动矢量V_{4_i}的水平分量V_{4_i_h}和垂直分量V_{4_i_v}；

所述运动矢量选取准则，不限于本实施例中所采用的均方误差最小准则，也可以采用其他准则，例如绝对差之和（SAD）最小准则等；

（3d）以步骤（3c）中全局运动矢量V_{4_i}所指向的第3级参考图像的下采样图像z_{3_r}中的像素位置为中心，使用第3级下采样图像z_{3_i}在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z_{3_i}相对于参考图像的下采样图像z_{3_r}的均方误差D_{3_i}：

$D_{3\_i}=\frac{1}{W_{3\_i}\&CenterDot;H_{3\_i}}{\mathrm{\&Sigma;}}_{x=-W_{3\_i}}^{W_{3\_i}}{\mathrm{\&Sigma;}}_{y=-H_{3\_i}}^{H_{3\_i}}(p(x,y)-q(x+u_{3\_i},y+v_{3\_i}){)}^2,$

-W_{3_i}≤u_{3_i}≤W_{3_i}，-H_{3_i}≤v_{3_i}≤H_{3_i}，

其中，W_{3_i}和H_{3_i}分别表示下采样图像z_{3_i}的宽度与高度，p(x,y)表示参考图像的下采样图像z_{3_r}中坐标为(x,y)的像素，q(x,y)表示下采样图像z_{3_i}中坐标为(x,y)的像素，u_{3_i}表示水平位移，v_{3_i}表示垂直位移；

选取使均方误差D_{3_i}最小的水平位移u_{3_i}和垂直位移v_{3_i}，对全局运动矢量V_{4_i}的水平分量V_{4_i_h}和垂直分量V_{4_i_v}进行修正，得到全局运动矢量V_{3_i}；

所述的全局运动矢量修正范围不限于本实施例中所采用的5×5像素大小的邻域，也可以采用不同尺度的邻域，例如7×7像素大小的邻域等；

（3e）以步骤（3d）中全局运动矢量V_{3_i}所指向的第2级参考图像的下采样图像z_{2_r}中的像素位置为中心，使用第2级下采样图像z_{2_i}在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z_{2_i}相对于参考图像的下采样图像z_{2_r}的均方误差D_{2_i}：

$D_{2\_i}=\frac{1}{W_{2\_i}\&CenterDot;H_{2\_i}}{\mathrm{\&Sigma;}}_{x=-W_{2\_i}}^{W_{2\_i}}{\mathrm{\&Sigma;}}_{y=-H_{2\_i}}^{H_{2\_i}}(p(x,y)-q(x+u_{2\_i},y+v_{2\_i}){)}^2,$

-W_{2_i}≤u_{2_i}≤W_{2_i}，-H_{2_i}≤v_{2_i}≤H_{2_i}，

其中，W_{2_i}和H_{2_i}分别表示下采样图像z_{2_i}的宽度与高度，p(x,y)表示参考图像的下采样图像z_{2_r}中坐标为(x,y)的像素，q(x,y)表示下采样图像z_{2_i}中坐标为(x,y)的像素，u_{2_i}表示水平位移，v_{2_i}表示垂直位移；

选取使均方误差D_{2_i}最小的水平位移u_{2_i}和垂直位移v_{2_i}，对全局运动矢量V_{3_i}的水平分量V_{3_i_h}和垂直分量V_{3_i_v}进行修正，得到全局运动矢量V_{2_i}；

（3f）以步骤（3d）中全局运动矢量V_{2_i}所指向的第1级参考图像的下采样图像z_{1_r}中的像素位置为中心，使用第1级下采样图像z_{1_i}在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z_{1_i}相对于参考图像的下采样图像z_{1_r}的均方误差D_{1_i}：

$D_{1\_i}=\frac{1}{W_{1\_i}\&CenterDot;H_{1\_i}}{\mathrm{\&Sigma;}}_{x=-W_{1\_i}}^{W_{1\_i}}{\mathrm{\&Sigma;}}_{y=-H_{1\_i}}^{H_{1\_i}}(p(x,y)-q(x+u_{1\_i},y+v_{1\_i}){)}^2,$

-W_{1_i}≤u_{1_i}≤W_{1_i}，-H_{1_i}≤v_{1_i}≤H_{1_i}，

其中，W_{1_i}和H_{1_i}分别表示下采样图像z_{1_i}的宽度与高度，p(x,y)表示参考图像的下采样图像z_{1_r}中坐标为(x,y)的像素，q(x,y)表示下采样图像z_{1_i}中坐标为(x,y)的像素，u_{1_i}表示水平位移，v_{1_i}表示垂直位移；

选取使均方误差D_{1_i}最小的水平位移u_{1_i}和垂直位移v_{1_i}，对全局运动矢量V_{2_i}的水平分量V_{2_i_h}和垂直分量V_{2_i_v}进行修正，得到全局运动矢量V_{1_i}；

（3g）以步骤（3f）中全局运动矢量V_{1_i}所指向的参考图像z_r的亮度分量Y中的像素位置为中心，使用输入图像z_i的亮度分量Y在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后输入图像z_i的亮度分量Y相对于参考图像z_r的亮度分量Y的均方误差D_i：

$D_i=\frac{1}{W\&CenterDot;H}{\mathrm{\&Sigma;}}_{x=-W}^W{\mathrm{\&Sigma;}}_{y=-H}^H(p(x,y)-q(x+u_i,y+v_i){)}^2,$

-W≤u_i≤W，-H≤v_i≤H，

其中，W和H分别表示输入图像z_i的宽度与高度，p(x,y)表示参考图像z_r的亮度分量Y中坐标为(x,y)的像素，q(x,y)表示输入图像z_i的亮度分量Y中坐标为(x,y)的像素，u_i表示水平位移，v_i表示垂直位移；

选取使均方误差D_i最小的水平位移u_i和垂直位移v_i，对全局运动矢量V_{1_i}的水平分量V_{1_i_h}和垂直分量V_{1_i_v}进行修正，得到全局运动矢量V_i。

步骤4：确定候选参考图像。

（4a）计算除参考图像z_r外，输入各帧图像的清晰度测度G_i与参考图像z_r的清晰度测度G_max之比R_i，若满足：R_i<0.875，则判定该R_i对应的第i帧输入图像z_i为模糊图像，并从输入图像中剔除该R_i对应的第i帧图像，保留剩余图像用于后续处理，并记剩余图像的帧数为N′，否则，保留所有输入图像用于后续处理；

以下实施例以保留所有输入图像用于后续处理为例，对于剩余图像帧数为N′的情况，可将后续处理中的帧数6替换为N′即可；

所述的模糊图像判定准则不限于本实施例所使用的准则，也可以采用其他的判定准则，例如直接使用清晰度测度G_i的大小作为模糊图像的判定准则；

（4b）对步骤（4a）中所保留的各帧图像，若存在清晰度测度G_i与参考图像z_r的清晰度测度G_max之比R_i满足：R_i>0.99，则选取该R_i对应的第i帧图像作为候选参考图像z_r′，否则，选择参考图像z_r作为基准图像

并选择全局运动矢量V_i作为基准全局运动矢量

并转入步骤6；

（4c）以候选参考图像z_r′为输入图像中的第1帧，参考图像z_r为输入图像中的第2帧为例，即r′=1,r=2，从步骤（3g）所得的全局运动矢量V_i中，选取候选参考图像z_r′相对参考图像z_r的全局运动矢量V_r′：

$V_{r^{\&prime;}}=V_1=[V_{1\_h},V_{1\_v}{]}^T,$

其中，V_{1_h}和V_{1_v}分别表示全局运动矢量V₁的水平分量和垂直分量；

（4d）根据步骤（4c）所得候选参考图像z_r′相对参考图像z_r的全局运动矢量V_r′，计算参考图像z_r相对于候选参考图像z_r′的候选全局运动矢量

${\hat{V}}_{r^{\&prime;}}={\hat{V}}_2=[-V_{1\_h},-V_{1\_v}{]}^T,$

（4e）计算输入的第i帧相对于候选参考图像z_r′的候选全局运动矢量

${\hat{V}}_i=V_i+{\hat{V}}_r=[V_{i\_h}-V_{1\_h},V_{i\_v}-V_{1\_v}{]}^T,i=\mathrm{3,4,5,6},$

使用

和

构成输入各帧相对于候选参考图像z_r′的候选全局运动矢量i=1,2,…,6且i≠r′。

步骤5：选取基准图像和基准全局运动矢量。

（5a）计算各帧图像相对于参考图像z_r的全局运动矢量测度S_r：

$S_r={\mathrm{\&Sigma;}}_{i=1}^6\left(\right|V_{i\_h}|+|V_{i\_v}\left|\right),$

其中，V_{i_h}表示第i帧输入图像相对于参考图像z_r的全局运动矢量V_i的水平分量，V_{i_v}表示第i帧输入图像相对于参考图像z_r的全局运动矢量V_i的垂直分量；

所述全局运动矢量测度不限于本实施例中所使用的全局运动矢量的水平分量与垂直分量的绝对值之和，也可以使用其他全局运动矢量测度，例如全局运动矢量的水平分量与垂直分量的平方和等；

（5b）计算各帧图像相对于候选参考图像z_r′的全局运动矢量测度S_r′：

$S_{r^{\&prime;}}={\mathrm{\&Sigma;}}_{i=1}^6\left(\right|{\hat{V}}_{i\_h}|+|{\hat{V}}_{i\_v}\left|\right),$

其中，

表示第i帧输入图像相对于候选参考图像z_r′的全局运动矢量

的水平分量，

表示第i帧输入图像相对于候选参考图像z_r′的全局运动矢量的垂直分量；

（5c）比较全局运动矢量测度S_r和S_r′，若满足条件：S_r<S_r′+25，则选取候选参考图像z_r′为基准图像

并选取与z_r′对应的全局运动矢量

作为基准全局运动矢量

否则，选取参考图像z_r为基准图像

并选取与z_r对应的全局运动矢量V_i作为基准全局运动矢量

步骤6：计算局部运动矢量。

（6a）将步骤（2a）中所述的1级下采样图像z_{1_i}划分为J=64个大小相同的搜索块，将基准图像

的1级下采样图像

划分为J=64个大小相同的基准搜索块；

所述的划分搜索块的个J=64数不限于本实施例中所使用的64，可根据输入图像的分辨率进行合理选择；

（6b）以基准全局运动矢量

所指向的第j个基准搜索块中的像素位置为参考，以该参考位置为中心，对第i帧输入图像的下采样图像z_{1_i}的第j个搜索块进行水平位移和垂直位移，计算所述第j个搜索块相对于第j个基准搜索块的均方误差D_{1_i_j}：

$D_{1\_i\_j}=\frac{1}{W_{S1}\&CenterDot;H_{S1}}{\mathrm{\&Sigma;}}_{x=-W_{S1}}^{W_{S1}}{\mathrm{\&Sigma;}}_{y=-H_{S1}}^{H_{S1}}(p(x,y)-q(x+u_{1\_i\_j},y+v_{1\_i\_j}){)}^2,$

-7≤u_{1_i_j}≤7，-7≤v_{1_i_j}≤7，

其中，W_S1和H_S1分别表示所述搜索块的宽度与高度，p(x,y)表示基准搜索块中坐标为(x,y)的像素，q(x,y)表示所述搜索块中坐标为(x,y)的像素，u_{1_i_j}表示水平位移，v_{1_i_j}表示垂直位移；

（6c）选取使均方误差D_{1_i_j}最小的水平位移u_{1_i_j}和垂直位移v_{1_i_j}，对基准全局运动矢量的水平分量

和垂直分量

进行修正，得到所述搜索块相对于基准搜索块的局部运动矢量V_{1_i_j}：

$V_{1\_i\_j}=[V_{1\_i\_j\_h},V_{1\_i\_j\_v}{]}^T={[{\tilde{V}}_{i\_h}+u_{1\_i\_j},{\tilde{V}}_{i\_v}+v_{1\_i\_j}]}^T,$

其中，V_{1_i_j_h}和V_{1_i_j_v}分别表示局部运动矢量V_{1_i_j}的水平分量和垂直分量；

（6d）将输入图像z_i的亮度分量Y划分为J=64个大小相同的搜索块，将基准图像

的亮度分量Y划分为J=64个大小相同的基准搜索块；

（6e）以局部运动矢量V_{1_i_j}所指向的基准图像的亮度分量Y的第j个基准搜索块中的像素位置为参考，以该参考位置为中心，对输入图像z_i的第j个搜索块进行水平位移和垂直位移，计算所述第j个搜索块相对于第j个基准搜索块的均方误差D_{i_j}：

$D_{i\_j}=\frac{1}{W_S\&CenterDot;H_S}{\mathrm{\&Sigma;}}_{x=-W_S}^{W_S}{\mathrm{\&Sigma;}}_{y=-H_S}^{H_S}(p(x,y)-q(x+u_{i\_j},y+v_{i\_j}){)}^2,$

-5≤u_{i_j}≤5，-5≤v_{i_j}≤5，

其中，W_S和H_S分别表示所述搜索块的宽度与高度，p(x,y)表示所述基准搜索块中坐标为(x,y)的像素，q(x,y)表示所述搜索块中坐标为(x,y)的像素，u_{i_j}表示水平位移，v_{i_j}表示垂直位移；

（6f）选取使均方误差D_{i_j}最小的水平位移u_{i_j}和垂直位移v_{i_j}，对局部运动矢量V_{1_i_j}的水平分量V_{1_i_j_h}和垂直分量V_{1_i_j_v}进行修正，得到所述搜索块相对于基准搜索块的局部运动矢量V_{i_j}：

V_{i_j}=[V_{i_j_h},V_{i_j_v}]^T=[V_{1_i_j_h}+u_{i_j},V_{1_i_j_v}+v_{i_j}]^T，

其中，V_{i_j_h}和V_{i_j_v}分别表示局部运动矢量V_{i_j}的水平分量和垂直分量。

步骤7：选取融合块组。

（7a）将步骤（6d）中每个搜索块划分为Q=108个大小相同的融合块，并将步骤（6d）中每个基准搜索块划分为Q=108个大小相同的基准融合块；

所述的划分融合块的个数Q不限于本实施例中所使用的108，可根据输入图像的分辨率进行合理选择；

（7b）以第k个基准融合块为参考，在各帧输入图像z_i中找出使用步骤（6f）中所得局部运动矢量V_{i_j}所指向的融合块，用所找出的各帧图像中的融合块构成一个融合块组；1≤k≤108；

（7c）使用步骤（2d）中所述方法，计算融合块组中每个融合块的清晰度测度G_i,j,k，将清晰度测度G_i,j,k小于阈值T_Blur的融合块判定为模糊融合块，并将该模糊融合块从所述融合块组中剔除，阈值T_Blur的计算公式如下：

T_Blur=μ_g-3σ_g，

其中，μ_g表示融合块组中所有融合块清晰度测度的均值，σ_g融合块组中所有融合块清晰度测度的方差；

以下实施例按照融合块组中不存在模糊融合块的情况为例，对于存在模糊融合块的情况，可在针对融合块组进行处理时，忽略对模糊融合块的处理，并对融合块组中融合块的个数进行相应调整；

（7d）经步骤（7c）剔除模糊融合块后的融合块组，即为所选取的融合块组。

步骤8：时域融合滤波去噪。

（8a）计算步骤（7d）中所得融合块组中，每个融合块与对应的基准融合块的时域相似性测度δ_{i_j_k}：

${\mathrm{\&delta;}}_{i\_j\_k}=\frac{1}{W_f\&CenterDot;H_f}{\mathrm{\&Sigma;}}_{x=1}^{W_f}{\mathrm{\&Sigma;}}_{y=1}^{H_f}(p_{i\_Y}(x,y)-q_Y(x,y)),$

其中，W_f和H_f分别表示融合块的宽度与高度，p_{i_Y}(x,y)表示第i个融合块亮度分量Y中坐标为(x,y)的像素值，q_Y(x,y)表示基准融合块亮度分量Y中坐标为(x,y)的像素值；1≤i≤6且i≠s,1≤j≤64,1≤k≤108；

（8b）依据步骤（8a）中所得的时域相似性测度δ_{i_j_k}，计算时域融合滤波权值ω_{i_j_k}：

${\mathrm{\&omega;}}_{i\_j\_k}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;16\\ (48-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/32,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;48\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>48\end{array}\right.;$

（8c）使用步骤（8b）中所得时域融合滤波权值ω_{i_j_k}，计算基准融合块的亮度分量Y的时域融合滤波结果

${\hat{q}}_Y(x,y)=\frac{1}{c_f}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}\&CenterDot;p_{i\_Y}(x,y)+q_Y(x,y)),$

其中，c_f表示归一化融合权值，

（8d）使用与亮度分量Y相同的时域融合滤波权值ω_{i_j_k}，分别计算基准融合块的红色差分量U的时域融合滤波结果

和绿色差分量V的时域融合滤波结果

${\hat{q}}_U(x,y)=\frac{1}{c_f}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}\&CenterDot;p_{i\_U}(x,y)+q_U(x,y)),$

${\hat{q}}_V(x,y)=\frac{1}{c_f}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}\&CenterDot;p_{i\_V}(x,y)+q_V(x,y)),$

其中，p_{i_U}(x,y)和p_{i_V}(x,y)分别表示融合块组中第i个融合块红色差分量U和绿色差分量V中坐标为(x,y)的像素值，q_U(x,y)和q_V(x,y)分别表示基准融合块红色差分量U和绿色差分量V中坐标为(x,y)的像素值；

（8e）使用所述基准融合块的亮度分量Y的时域融合滤波结果

红色差分量U的时域融合滤波结果

和绿色差分量V的时域融合滤波结果

构成时域融合滤波图像z_fu中与基准融合块相对应位置的融合结果，所有基准融合块的时域融合滤波结果共同构成时域融合滤波图像z_fu。

步骤9：局部运动区域判定。

（9a）使用步骤（7d）中所得融合块组，计算融合块组中每个融合块与基准融合块的对应位置亮度分量Y的像素差值d_{i_Y}(x,y)：

d_{i_Y}(x,y)=p_{i_Y}(x,y)-q_Y(x,y)，

其中，p_{i_Y}(x,y)表示融合块组中第i个融合块亮度分量Y中坐标为(x,y)的像素值，q_Y(x,y)表示基准融合块亮度分量Y中坐标为(x,y)的像素值；

（9b）统计d_{i_Y}(x,y)大于阈值24的像素个数占融合块总像素个数的百分比p_{i_j_k}，计算融合块组中该百分比p_{i_j_k}的标准差σ_j,k；

（9c）若所述标准差σ_j,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域。

步骤10：局部运动区域时域融合滤波修正。

（10a）对步骤（9）中所判定的局部运动区域，依据时域相似性测度δ_{i_j_k}，按照下式计算局部运动区域的时域融合滤波权值ω′_{i_j_k}:

${\mathrm{\&omega;}}_{i\_j\_k}^{\&prime;}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;6\\ (12-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/6,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;12\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>12\end{array}\right.;$

（10b）用局部运动区域的时域滤波权值ω′_{i_j_k}计算局部运动区域亮度分量Y的时域融合滤波结果

${\tilde{q}}_Y(x,y)=\frac{1}{c_f^{\&prime;}}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}^{\&prime;}\&CenterDot;p_{i\_Y}(x,y)+q_Y(x,y)),$

其中，

表示局部运动区域对应的基准融合块的亮度分量Y中坐标为(x,y)的像素值，

表示局部运动区域的归一化的融合权值，

（10c）使用与亮度分量Y相同的时域滤波权值计算局部运动区域的红色差分量U的时域融合滤波结果

和绿色差分量V的时域融合滤波结果 ${\tilde{q}}_V(x,y):$

${\tilde{q}}_U(x,y)=\frac{1}{c_f^{\&prime;}}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}^{\&prime;}\&CenterDot;p_{i\_U}(x,y)+q_U(x,y)),$

${\tilde{q}}_V(x,y)=\frac{1}{c_f^{\&prime;}}({\mathrm{\&Sigma;}}_i{\mathrm{\&omega;}}_{i\_j\_k}^{\&prime;}\&CenterDot;p_{i\_V}(x,y)+q_V(x,y));$

（10d）使用所述局部运动区域亮度分量Y的时域融合滤波结果

红色差分量U的时域融合滤波结果

和绿色差分量V的时域融合滤波结果

构成局部运动区域的时域融合结果。

（10e）使用对局部运动区域的融合结果，覆盖时域融合滤波图像z_fu中对应区域的融合结果，实现对时域融合滤波图像z_fu中局部运动区域的修正，得到最终的时域融合图像

步骤11：空域非局部均值滤波去噪。

（11a）采用步骤（2a）中所述的下采样方法，对步骤（10e）所得时域融合图像的亮度分量Y进行1级下采样，得到该时域融合图像

下采样图像Y_{1_fu}；

（11b）依照非局部均值滤波算法，设定所述下采样图像Y_{1_fu}中的相似块大小为5×5像素，搜索窗的范围为7×7像素；

（11c）以所述下采样图像Y_{1_fu}中的像素p为例，按照图3所示方式进行遍历搜索，计算像素p的非局部均值滤波权值；

图3中，A(p)表示以像素p为中心的相似块，B(p)表示以像素p为中心的搜索窗，A(q)表示以像素q为中心的搜索窗B(p)中的相似块，对B(p)内的像素按照图3中箭头方向进行逐像素遍历，计算像素p的非局部均值滤波权值w(p,q)：

$w(p,q)=\exp \left(\frac{{\left|\right|A\left(p\right)-A\left(q\right)\left|\right|}_2^2}{h^2}\right),$

其中，h表示滤波强度参数，对步骤9中所判定的局部运动区域，h的取值为20，对于其他区域，h的取值为10，

表示相似块A(p)和A(q)的欧式距离；

（11d）以时域融合图像中与像素p对应的像素p′为例，使用像素p的非局部均值滤波权值w(p,q)，计算像素p′亮度分量Y的空域加权滤波值

${\hat{p}}_Y^{\&prime;}=\frac{1}{W\left(p\right)}{\mathrm{\&Sigma;}}_{q\&Element;B\left(q\right)}{p}_Y^{\&prime;}\&CenterDot;w(p,q),$

其中，

表示像素p′的亮度分量，W(p)表示归一化滤波权值，按照如下公式计算：

$W\left(p\right)={\mathrm{\&Sigma;}}_{q\&Element;B\left(q\right)}w(p,q);$

（11）使用像素p的非局部均值滤波权值w(p,q)，分别计算像素p′的红色差分量U的空域加权滤波值

和绿色差分量V的空域加权滤波值

${\hat{p}}_U^{\&prime;}=\frac{1}{W\left(p\right)}{\mathrm{\&Sigma;}}_{q\&Element;B\left(q\right)}{p}_U^{\&prime;}\&CenterDot;w(p,q),$

${\hat{p}}_v^{\&prime;}=\frac{1}{W\left(p\right)}{\mathrm{\&Sigma;}}_{q\&Element;B\left(q\right)}{p}_v^{\&prime;}\&CenterDot;w(p,q),$

其中，

表示像素p′的红色差分量，表示像素的绿色差分量；

（11e）使用像素p′亮度分量Y的空域加权滤波值

红色差分量U的空域加权滤波值

以及绿色差分量V的空域加权滤波值

构成空域滤波图像z_f中与像素p′相对应的非局部均值滤波结果，由所有像素的非局部均值滤波结果共同构成空域滤波图像z_f。

步骤12：亮度校正和颜色校正。

（12a）计算步骤（11e）中所得的空域滤波图像z_f的平均亮度

（12b）根据空域滤波图像z_f的平均亮度

按照如下公式生成亮度校正曲线：

$\left\{\begin{array}{ccc}{y}_1=160.5e^{0.00195x}-156e^{-0.008996x}& 0\&le;{\stackrel{\&OverBar;}{Y}}_f<86& \\ y_2=(0.5(\frac{y_1}{x}-1)+1)x& 86\&le;{\stackrel{\&OverBar;}{Y}}_f<117& ,\\ y_3=x& 117\&le;{\stackrel{\&OverBar;}{Y}}_f\&le;255& \end{array}\right.,$

其中，x表示亮度校正曲线的水平坐标，y_l表示平均亮度z_f落入不同范围时亮度校正曲线的垂直坐标，l=1,2,3；

（12c）根据空域滤波图像z_f的平均亮度

选择相应亮度校正曲线，分别计算空域滤波图像z_f中像素p(x,y)亮度分量校正值

红色差分量校正值

和绿色差分量校正值

$\left\{\begin{array}{c}{\hat{p}}_Y(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;p_Y(x,y)\\ {\hat{p}}_U(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;(p_U(x,y)-128)+128\\ {\hat{p}}_V(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;\left(p_V(x,y)128\right)+128\end{array}\right.,$

其中，p_Y(x,y)表示像素p(x,y)的亮度分量，p_U(x,y)表示像素p(x,y)的红色差分量，p_V(x,y)表示像素p(x,y)的绿色差分量，α_l表示亮度校正曲线在x=pY(x,y)点的斜率，其计算公式为：

（12d）使用空域滤波图像z_f所有像素亮度分量、红色差分量以及绿色差分量的校正结果构成校正图像z_c。

步骤13：对比度增强。

（13a）对校正图像z_c进行高斯滤波，得到校正图像z_c的高斯滤波结果图像z_g；

（13b）利用高斯滤波结果图像z_g中的像素q_g(x,y)，计算校正图像z_c中像素q(x,y)的局部对比度增强结果

$\left\{\begin{array}{cc}\hat{q}(x,y)=q(x,y)& \frac{q(x,y)}{q_g(x,y)}\&GreaterEqual;1\\ \hat{q}(x,y)=q(x,y)\&CenterDot;{\left(\frac{q(x,y)}{q_g(x,y)}\right)}^{\mathrm{\&beta;}}& \frac{q(x,y)}{q_g(x,y)}<1\end{array}\right.,$

其中，β表示对比度增强因子，取值为0.3。

至此，获得最终的多帧去噪输出图像z_out。

以图4所示的低光照环境下采集的6帧自然图像的实验为例，得到的最终去噪图像是图7。

本发明的优点由以下仿真实现进一步说明。

1.仿真条件

本发明采用图4所示的低光照环境下采集的针对同一场景的6帧自然图像作为测试图像，其中图4(a)～图4(f)分别为输入的第1～6帧测试图像，软件平台为Visual Studio2010。

2.仿真内容与结果：

仿真1，使用现有的非局部均值滤波算法对图4所示的6帧图像分别进行非局部均值滤波，对滤波所得的6帧图像进行平均，所得的结果如图5所示。

仿真2，使用现有的块匹配三维协同滤波算法对图4所示的6帧图像分别进行滤波，对滤波所得的6帧图像进行平均，所得的结果如图6所示。

仿真3，使用本发明方法对图4所示的6帧图像分别进行非局部均值滤波，滤波所得的结果如图7所示。

为了便于比较所述实验方法的性能，对图4(a)所示的第一帧输入图像的局部区域进行原始分辨率显示，如图8所示；对图5中与图4(a)对应的局部区域进行原始分辨率显示，如图9所示；对图6中与图4(a)对应的区域进行原始分辨率显示，如图10所示；对图7中与图4(a)对应的局部区域进行原始分辨率显示，如图11所示。

从图9与图8的对比可以看出，现有的非局部均值滤波算法对低光照环境下采集的自然图像中噪声的抑制十分有限，并且在图像的边缘存在模糊；

从图10与图8的对比可以看出，相比于非局部均值滤波算法，现有的块匹配三维协同滤波算法对图像中噪声的抑制性能有了显著提升，但场景中的细节纹理区域，仍存在模糊；

从图11与图8的对比可以看出，相比于现有的非局部均值滤波算法和块匹配三维协同滤波算法，本发明方法能够有效抑制低光照环境下采集图像中的噪声，并且边缘与细节的保持度比较高，在去除噪声的同时可以有效的避免细节模糊。

Claims (8)

1.一种基于空域和时域联合滤波的多帧数字图像去噪方法，包括如下步骤：

（1）输入N帧低光照环境下采集的针对同一场景的自然图像z_i,i=1,2,…,N；

（2）计算各帧图像的平均水平梯度

和平均垂直梯度

将

与

之和作为图像清晰度测度G_i，选取G_i最大的一帧，作为参考图像z_r,1≤r≤N，并将该参考图像的清晰度测度记为G_max；

（3）计算各帧图像相对于参考图像z_r的全局运动矢量V_i；

（4）计算各帧图像清晰度测度G_i与参考图像清晰度测度G_max的比值R_i，将R_i大于阈值0.99所对应的第i帧图像记为候选参考图像z_r'，计算各帧图像相对于该候选参考图像z_r的候选全局运动矢量

（5）计算全局运动矢量V_i和候选全局运动矢量

的测度，并依据该测度从参考图像z_r和候选参考图像z_r'中选取基准图像

并将该基准图像对应的全局运动矢量作为基准全局运动矢量

（6）将输入的各帧图像都划分成大小相同的J个搜索块，利用上述基准全局运动矢量

计算各帧图像中每个搜索块相对于基准图像

中对应搜索块的局部运动矢量V_{i_j},1≤i≤N且i≠s,1≤j≤J；

（7）将各帧图像中的每个搜索块都划分为大小相同的Q个融合块，使用步骤（6）中所得局部运动矢量V_{i_j}，在基准图像

中找出与各帧图像中的每个融合块所对应的基准融合块，选取各帧图像中与同一基准融合块对应的融合块构成时域融合块组；

（8）分别计算每个时域融合块组中每个融合块与其对应基准融合块的像素差值，并将这些像素差值的平均值作为该时域融合块组中每个融合块与其对应基准融合块的时域相似性测度δ_{i_j_k},1≤i≤N且i≠s,1≤j≤J,1≤k≤Q；

（9）依据步骤（8）中所得的时域相似性测度δ_{i_j_k}，计算时域融合滤波权值ω_{i_j_k}：

${\mathrm{\&omega;}}_{i\_j\_k}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;16\\ (48-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/32,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;48\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>48\end{array}\right.;$

（10）用所述时域融合滤波权值ω_{i_j_k}对每个时域融合块组所对应的基准融合块进行时域加权融合滤波，得到时域融合滤波图像z_fu；

（11）依据步骤（8）中所得各个时域融合块组中每个融合块与其对应基准融合块的像素差值，统计像素差值大于阈值24的像素个数占融合块总像素个数的百分比，计算这些百分比的标准差σ_j,k，若σ_j,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域；

（12）对步骤（11）中所判定的局部运动区域，依据时域相似性测度δ_{i_j_k}，按照下式计算局部运动区域的时域融合滤波权值ω′_{i_j_k}:

${{\mathrm{\&omega;}}^{\&prime;}}_{i\_j\_k}=\left\{\begin{array}{cc}1& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;6\\ (12-|{\mathrm{\&delta;}}_{i\_j\_k}\left|\right)/6,& 16<\left|{\mathrm{\&delta;}}_{i\_j\_k}\right|\&le;12\\ 0& \left|{\mathrm{\&delta;}}_{i\_j\_k}\right|>12\end{array}\right.;$

（13）用局部运动区域的时域融合滤波权值ω′_{i_j_k}对所判定的局部运动区域所对应的基准融合块进行加权融合滤波，使用该融合结果覆盖时域融合滤波图像z_fu中对应位置的融合结果，得到最终的时域融合图像

并对该时域融合图像

进行单帧空域非局部均值滤波，得到空域滤波图像z_f；

（14）计算空域滤波图像z_f的平均亮度

根据该平均亮度

对空域滤波图像z_f进行亮度校正和色度校正，得到亮度和色度校正图像z_c，对z_c再进行局部对比度增强，得到最终的去噪结果图像z_out。

2.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤（2）所述的计算各帧图像的平均水平梯度

和平均垂直梯度

利用如下公式计算：

${\stackrel{\&OverBar;}{G}}_{i\_h}=\frac{1}{W\&CenterDot;H}{\mathrm{\&Sigma;}}_{x=1}^H{\mathrm{\&Sigma;}}_{y=1}^W{G}_{i\_h}(x,y),$

${\stackrel{\&OverBar;}{G}}_{i\_v}=\frac{1}{W\&CenterDot;H}{\mathrm{\&Sigma;}}_{x=1}^H{\mathrm{\&Sigma;}}_{y=1}^W{G}_{i\_v}(x,y),$

其中，W和H分别表示输入图像的宽度和高度，G_{i_h}(x,y)表示输入图像z_i中坐标为(x,y)的像素的水平梯度，G_{i_v}(x,y)表示输入图像z_i中坐标为(x,y)的像素的垂直梯度；水平梯度G_{i_h}(x,y)和垂直梯度G_{i_v}(x,y)的计算公式如下：

G_{i_h}(x,y)=S_h*A(x,y)，

G_{i_v}(x,y)=S_v*A(x,y)，

其中，A(x,y)表示图像中以坐标(x,y)为中心的3×3像素大小的图像块，S_h和S_v分别表示3×3大小的水平和垂直Sobel算子：

$S_h=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right],S_v=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right],$

“*”表示卷积运算。

3.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤（3）所述的计算各帧图像相对于参考图像z_r的全局运动矢量V_i，按照如下步骤进行：

（3.1）对输入的N帧图像z_i,i=1,2,…,N进行M级下采样，记z_i的各级下采样图像分别为z_{m_i}，记参考图像z_r的各级下采样图像为z_{m_r},1≤m≤M；

（3.2）使用输入图像z_i的第M级下采样图像z_{M_i}相对于参考图像z_r的第M级下采样图像z_{M_r}进行全图运动搜索，选取使均方误差最小的运动矢量作为下采样图像z_{M_i}相对于参考图像z_r的下采样图像z_{M_r}的第M级全局运动矢量V_{M_i}；

（3.3）以第M级全局运动矢量V_{M_i}所指向的参考图像z_r的第M-1级下采样图像z_{M-1_r}中的像素位置为中心，使用输入图像z_i的第M-1级下采样图像z_{M_i}相对于所述中心的3×3像素大小邻域内进行运动搜索，选取使均方误差最小的运动矢量对第M级全局运动矢量V_{M_i}修正，作为下采样图像z_{M-1_i}相对于参考图像z_r的下采样图像z_{M-1_i}的第M-1级全局运动矢量V_{M-1_i}；

（3.4）依照步骤（3.3）对全局运动矢量V_{M_i}进行逐级修正，直至获得输入图像z_i相对于参考图像z_r的全局运动矢量V_i。

4.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤（5）所述的选取基准图像

和基准全局运动矢量

按照如下步骤进行：

（5.1）计算各帧图像清晰度测度G_i与参考图像清晰度测度G_max的比值R_i，与阈值0.875相比较，若满足：

R_i<0.875，

则判定该比值R_i对应的第i帧图像z_i为模糊图像，并从输入的N帧图像中剔除该比值R_i对应的第i帧图像，保留剩余图像用于后续处理，记剩余图像的帧数为N'；

（5.2）根据全局运动矢量V_i的水平分量V_{i_h}与垂直分量V_{i_v}，计算全局运动矢量测度S_r：

$S_r={\mathrm{\&Sigma;}}_{i=1}^{N^{\&prime;}}\left(\right|V_{i\_h}|+|V_{i\_v}\left|\right);$

根据候选全局运动矢量的水平分量和垂直分量

计算全局运动矢量测度S_r'：

$S_{r^{\&prime;}}={\mathrm{\&Sigma;}}_{i=1}^{N^{\&prime;}}\left(\right|{\hat{V}}_{i\_h}|+|{\hat{V}}_{i\_v}\left|\right);$

（5.3）比较全局运动矢量测度S_r和S_r'，若满足条件：

S_r<S_r'+25，

则选取候选参考图像z_r'为基准图像

并选取与z_r'对应的全局运动矢量作为基准全局运动矢量

否则选取参考图像z_r为基准图像

并选取与z_r对应的全局运动矢量V_i作为基准全局运动矢量

5.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤（6）所述计算各帧图像中每个搜索块相对于基准图像

中对应搜索块的局部运动矢量V_{i_j}，按照如下步骤进行：

（6.1）将输入的N帧图像z_i,i=1,2,…,N都划分成大小相同的J个搜索块，记第i帧输入图像z_i的第j个搜索块为B_{i_j}，记基准图像

中的第j个搜索块为基准搜索块B_{s_j},1≤j≤J；

（6.2）以基准全局运动矢量

所指向基准搜索块B_{s_j}中的像素位置为中心，使用输入图像z_i的第j个搜索块B_{i_j}在该中心的5×5像素大小邻域内进行运动搜索，以均方误差最小为准则对基准全局运动矢量

进行修正，获得搜索块B_{i_j}相对于基准搜索块B_{s_j}的局部运动矢量V_{i_j}。

6.根据权利要求1所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中所述步骤（13）中对时域融合图像

进行单帧空域非局部均值滤波，按照如下步骤进行：

（13.1）对时域融合图像

的亮度分量Y_fu进行1级下采样，得到下采样图像Y_{1_fu}；

（13.2）依照非局部均值滤波算法，设定下采样图像Y_{1_fu}中的相似块大小为5×5像素，搜索窗的范围为7×7像素；

（13.3）设A(p)表示以p为中心像素的相似块，设B(p)表示以p为中心像素的搜索窗，设A(q)表示以q为中心像素的B(p)中的相似块，对B(p)内的像素进行遍历，计算像素p的非局部均值滤波权值：

$w(p,q)=\exp (-\frac{\left|\right|A\left(p\right)-A\left(q\right){\left|\right|}_2^2}{h^2}),$

其中，h表示滤波强度参数，

表示相似块A(p)和A(q)的欧式距离；

（13.4）设p'表示与像素p相对应的时域融合图像

中的像素，计算像素p'亮度分量的加权滤波值

${\hat{P}}_Y^{\&prime;}=\frac{1}{W\left(p\right)}{\mathrm{\&Sigma;}}_{q\&Element;B\left(p\right)}{\hat{p}}_Y^{\&prime;}\&CenterDot;w(p,q),$

其中，p'_Y表示像素p'的亮度分量，W(p)表示归一化的滤波权值，按照如下公式计算：

$W\left(p\right)={\mathrm{\&Sigma;}}_{q\&Element;B\left(p\right)}w(p,q);$

（13.5）使用亮度分量Y_fu的滤波权值对时域融合图像

的两个色差分量进行加权滤波，并用亮度分量和两个色差分量的滤波结果构成空域滤波图像z_f。

7.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，所述步骤（14）中利用平均亮度

对空域滤波图像z_f进行亮度校正和颜色校正，按如下步骤进行：

（14.1）根据空域滤波图像z_f的平均亮度

按照如下公式生成亮度校正曲线：

$\left\{\begin{array}{ccc}{y}_1=160.5e^{0.00195x}-156e^{-0.008996x}& 0\&le;{\stackrel{\&OverBar;}{Y}}_f<86& \\ y_2=(0.5(\frac{y_1}{x}-1)+1)x& 86\&le;{\stackrel{\&OverBar;}{Y}}_f<117& ,\\ y_3=x& 117\&le;{\stackrel{\&OverBar;}{Y}}_f\&le;255& \end{array}\right.,$

其中，x表示亮度校正曲线的水平坐标，y_l表示平均亮度

落入不同范围时亮度校正曲线的垂直坐标，l=1,2,3；

（14.2）根据空域滤波图像z_f的平均亮度

选择相应亮度校正曲线，分别计算空域滤波图像z_f中像素p(x,y)亮度分量校正值

红色差分量校正值

和绿色差分量校正值

$\left\{\begin{array}{c}{\hat{P}}_Y(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;p_Y(X,Y)\\ {\hat{p}}_U(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;(P_U(x,y)-128)+128\\ {\hat{P}}_V(x,y)={\mathrm{\&alpha;}}_l\&CenterDot;(P_V(x,y)-128)+128\end{array}\right.,$

其中，p_Y(x,y)表示像素p(x,y)的亮度分量，p_U(x,y)表示像素p(x,y)的红色差分量，p_V(x,y)表示像素p(x,y)的绿色差分量，α_l表示亮度校正曲线在x=p_Y(x,y)的斜率，其计算公式为： ${\mathrm{\&alpha;}}_l=\frac{y_l}{x},l=\mathrm{1,2,3};$

（14.3）使用空域滤波图像z_f所有像素亮度分量、红色差分量以及绿色差分量的校正结果构成校正图像z_c。

8.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中所述步骤（14）中对校正图像z_c的进行局部对比度增强，按照如下步骤进行：

（14a）对校正图像z_c进行高斯滤波，得到校正图像z_c的高斯滤波结果图像z_g；

（14b）利用高斯滤波结果图像z_g中的坐标为(x,y)的像素q_g(x,y)，计算校正图像z_c中像素q(x,y)的局部对比度增强结果

$\left\{\begin{array}{cc}\hat{q}(x,y)=q(x,y)& \frac{q(x,y)}{q_g(x,y)}\&GreaterEqual;1\\ \hat{q}(x,y)=q(x,y)\&CenterDot;{\left(\frac{q(x,y)}{q_g(x,y)}\right)}^{\mathrm{\&beta;}}& \frac{q(x,y)}{q_g(x,y)}<1\end{array}\right.,$

其中，β表示对比度增强因子，取值为0.3。

Images