本发明公开了一种基于空域和时域联合滤波的多帧数字图像去噪方法。其方案是：输入低光照环境下采集的同一场景的多帧数字图像；选取多帧图像中最清晰的一张作为基准图像；对多帧数字图像进行全局和局部配准；以基准图像为参考，在空域和时间域中建立与基准图像中局部区域相似的相似组，并采用相似组对基准图像中的局部区域进行协同滤波；利用滤波后图像中的亮度信息以及像素分布信息，对滤波后图像依次进行颜色校正和对比度增强，得到去噪后的图像。本发明能有效抑制数字图像中的噪声和运动模糊等因素对于图像质量的影响，并在抑制噪声的同时有效的保留图像中的细节纹理，可用于提高低光照环境下数字图像采集设备中的图像质量。
1.一种基于空域和时域联合滤波的多帧数字图像去噪方法，包括如下步骤：(1)输入N帧低光照环境下采集的针对同一场景的自然图像zi,i＝1,2,…,N；
(2)计算各帧图像的平均水平梯度 和平均垂直梯度 将 与 之和作为图像清晰度测度Gi，选取Gi最大的一帧，作为参考图像zr,1≤r≤N，并将该参考图像的清晰度测度记为Gmax；
(3)计算各帧图像相对于参考图像zr的全局运动矢量Vi；
(4)计算各帧图像清晰度测度Gi与参考图像清晰度测度Gmax的比值Ri，将Ri大于阈值
0.99所对应的第i帧图像记为候选参考图像zr'，计算各帧图像相对于该候选参考图像zr'的候选全局运动矢量(5)计算全局运动矢量Vi和候选全局运动矢量 的测度，并依据该测度从参考图像zr和候选参考图像zr'中选取基准图像 并将该基准图像 对应的全局运动矢量作为基准全局运动矢量(6)将输入的各帧图像都划分成大小相同的J个搜索块，利用上述基准全局运动矢量 计算各帧图像中每个搜索块相对于基准图像 中对应搜索块的局部运动矢量Vi_j,1≤i≤N且i≠s,1≤j≤J；
(7)将各帧图像中的每个搜索块都划分为大小相同的Q个融合块，使用步骤(6)中所得局部运动矢量Vi_j，在基准图像 中找出与各帧图像中的每个融合块所对应的基准融合块，选取各帧图像中与同一基准融合块对应的融合块构成时域融合块组；
(8)分别计算每个时域融合块组中每个融合块与其对应基准融合块的像素差值，并将这些像素差值的平均值作为该时域融合块组中每个融合块与其对应基准融合块的时域相似性测度δi_j_k,1≤i≤N且i≠s,1≤j≤J,1≤k≤Q；
(9)依据步骤(8)中所得的时域相似性测度δi_j_k，计算时域融合滤波权值ωi_j_k：(10)用所述时域融合滤波权值ωi_j_k对每个时域融合块组所对应的基准融合块进行时域加权融合滤波，得到时域融合滤波图像zfu；
(11)依据步骤(8)中所得各个时域融合块组中每个融合块与其对应基准融合块的像素差值，统计像素差值大于阈值24的像素个数占融合块总像素个数的百分比，计算这些百分比的标准差σj,k，若σj,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域；
(12)对步骤(11)中所判定的局部运动区域，依据时域相似性测度δi_j_k，按照下式计算局部运动区域的时域融合滤波权值ω'i_j_k：(13)用局部运动区域的时域融合滤波权值ω'i_j_k对所判定的局部运动区域所对应的基准融合块进行加权融合滤波，使用该融合结果覆盖时域融合滤波图像zfu中对应位置的融合结果，得到最终的时域融合图像 并对该时域融合图像 进行单帧空域非局部均值滤波，得到空域滤波图像zf；
(14)计算空域滤波图像zf的平均亮度 根据该平均亮度 对空域滤波图像zf进行亮度校正和色度校正，得到亮度和色度校正图像zc，对zc再进行局部对比度增强，得到最终的去噪结果图像zout。
2.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤(2)所述的计算各帧图像的平均水平梯度 和平均垂直梯度 利用如下公式计算：其中，W和H分别表示输入图像的宽度和高度，Gi_h(x,y)表示输入图像zi中坐标为(x,y)的像素的水平梯度，Gi_v(x,y)表示输入图像zi中坐标为(x,y)的像素的垂直梯度；
水平梯度Gi_h(x,y)和垂直梯度Gi_v(x,y)的计算公式如下：Gi_h(x,y)＝Sh*A(x,y)，
其中，A(x,y)表示图像中以坐标(x,y)为中心的3×3像素大小的图像块，Sh和Sv分别表示3×3大小的水平和垂直Sobel算子：“*”表示卷积运算。
3.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤(3)所述的计算各帧图像相对于参考图像zr的全局运动矢量Vi，按照如下步骤进行：(3.1)对输入的N帧图像zi,i＝1,2,…,N进行M级下采样，记zi的各级下采样图像分别为zm_i，记参考图像zr的各级下采样图像为zm_r,1≤m≤M；
(3.2)使用输入图像zi的第M级下采样图像zM_i相对于参考图像zr的第M级下采样图像zM_r进行全图运动搜索，选取使均方误差最小的运动矢量作为下采样图像zM_i相对于参考图像zr的下采样图像zM_r的第M级全局运动矢量VM_i；
(3.3)以第M级全局运动矢量VM_i所指向的参考图像zr的第M-1级下采样图像zM-1_r中的像素位置为中心，使用输入图像zi的第M-1级下采样图像zM-1_i相对于所述中心的3×3像素大小邻域内进行运动搜索，选取使均方误差最小的运动矢量对第M级全局运动矢量VM_i修正，作为下采样图像zM-1_i相对于参考图像zr的下采样图像zM-1_r的第M-1级全局运动矢量VM-1_i；
(3.4)依照步骤(3.3)对第M级全局运动矢量VM_i进行逐级修正，直至获得输入图像zi相对于参考图像zr的全局运动矢量Vi。
4.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤(5)所述的选取基准图像 和基准全局运动矢量 按照如下步骤进行：(5.1)计算各帧图像清晰度测度Gi与参考图像清晰度测度Gmax的比值Ri，与阈值0.875相比较，若满足：Ri<0.875，
则判定该比值Ri对应的第i帧图像zi为模糊图像，并从输入的N帧图像中剔除该比值Ri对应的第i帧图像，保留剩余图像用于后续处理，记剩余图像的帧数为N'；
(5.2)根据全局运动矢量Vi的水平分量Vi_h与垂直分量Vi_v，计算全局运动矢量测度Sr：根据候选全局运动矢量 的水平分量 和垂直分量 计算候选全局运动矢量测度Sr'：(5.3)比较全局运动矢量测度Sr和候选全局运动矢量测度Sr'，若满足条件：Sr
则选取候选参考图像zr'为基准图像 并选取与zr'对应的全局运动矢量 作为基准全局运动矢量 否则选取参考图像zr为基准图像 并选取与zr对应的全局运动矢量Vi作为基准全局运动矢量
5.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中步骤(6)所述计算各帧图像中每个搜索块相对于基准图像 中对应搜索块的局部运动矢量Vi_j，按照如下步骤进行：(6.1)将输入的N帧图像zi,i＝1,2,…,N都划分成大小相同的J个搜索块，记第i帧输入图像zi的第j个搜索块为Bi_j，记基准图像 中的第j个搜索块为基准搜索块Bs_j,1≤j≤J；
(6.2)以基准全局运动矢量 所指向基准搜索块Bs_j中的像素位置为中心，使用输入图像zi的第j个搜索块Bi_j在该中心的5×5像素大小邻域内进行运动搜索，以均方误差最小为准则对基准全局运动矢量 进行修正，获得搜索块Bi_j相对于基准搜索块Bs_j的局部运动矢量Vi_j。
6.根据权利要求1所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中所述步骤(13)中对时域融合图像 进行单帧空域非局部均值滤波，按照如下步骤进行：(13.1)对时域融合图像 的亮度分量Yfu进行1级下采样，得到下采样图像Y1_fu；
(13.2)依照非局部均值滤波算法，设定下采样图像Y1_fu中的相似块大小为5×5像素，搜索窗的范围为7×7像素；
(13.3)设A(p)表示以p为中心像素的相似块，设B(p)表示以p为中心像素的搜索窗，设A(q)表示以q为中心像素的B(p)中的相似块，对B(p)内的像素进行遍历，计算像素p的非局部均值滤波权值：其中，h表示滤波强度参数， 表示相似块A(p)和A(q)的欧式距离；
(13.4)设p'表示与像素p相对应的时域融合图像 中的像素，计算像素p'亮度分量的加权滤波值其中，p'Y表示像素p'的亮度分量，W(p)表示归一化的滤波权值，按照如下公式计算：W(p)＝∑q∈B(p)w(p,q)；
(13.5)使用亮度分量Yfu的滤波权值对时域融合图像 的红色色差分量和绿色色差分量进行加权滤波，并用亮度分量、红色色差分量和绿色色差分量的滤波结果构成空域滤波图像zf。
7.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，所述步骤(14)中利用平均亮度 对空域滤波图像zf进行亮度校正和色度校正，按如下步骤进行：(14.1)根据空域滤波图像zf的平均亮度 按照如下公式生成亮度校正曲线：其中，x表示亮度校正曲线的水平坐标，yl表示平均亮度 落入不同范围时亮度校正曲线的垂直坐标，l＝1,2,3；
(14.2)根据空域滤波图像zf的平均亮度 选择相应亮度校正曲线，分别计算空域滤波图像zf中像素p(x,y)亮度分量校正值 红色差分量校正值 和绿色差分量校正值其中，pY(x,y)表示像素p(x,y)的亮度分量，pU(x,y)表示像素p(x,y)的红色差分量，pV(x,y)表示像素p(x,y)的绿色差分量，αl表示亮度校正曲线在x＝pY(x,y)的斜率，其计算公式为：(14.3)使用空域滤波图像zf所有像素亮度分量、红色差分量以及绿色差分量的校正结果构成校正图像zc。
8.根据权利要求1中所述的基于空域和时域联合滤波的多帧数字图像去噪方法，其中所述步骤(14)中对校正图像zc的进行局部对比度增强，按照如下步骤进行：(14a)对校正图像zc进行高斯滤波，得到校正图像zc的高斯滤波结果图像zg；
(14b)利用高斯滤波结果图像zg中的坐标为(x,y)的像素qg(x,y)，计算校正图像zc中像素q(x,y)的局部对比度增强结果其中，β表示对比度增强因子，取值为0.3。
基于空域和时域联合滤波的多帧数字图像去噪方法
技术领域
[0001] 本发明属于图像处理技术领域，特别是涉及图像去噪方法，可用于消费类数字图像采集设备、医学影像、天文影像以及环境变化评估等领域的数字图像预处理。
背景技术
[0002] 图像去噪是图像处理领域的一个热点问题，也是一个极具挑战的研究方向。图像中的噪声会阻碍人们对图像内容的理解，采用图像去噪可以有效的抑制图像中噪声对图像质量的影响，提高人们对图像内容的认知程度，以便对图像作进一步的处理。
[0003] 根据数字图像的特点及其统计特性，多年来已经有很多学者提出了多种不同的去噪算法，按照其实现方法大致可以分为空域和频域两类，且其中大多都是基于图像局部信息的平滑处理。基于图像局部信息的平滑处理会使图像丢失很多细节信息，去噪效果十分不理想。2005年A.Buades，B.Coll等人对双边滤波去噪算法进行了改进，提出了一种非局部均值滤波的去噪方法。该方法打破了以往采用的图像空间域中“局部平滑”的思想，转而利用自然图像中广泛存在的“非局部”的空域相似性，通过在整幅图像中搜索相似块进行加权平均，以达到抑制噪声的效果。2007年，K.Dabov，Alessandro Foi等人在非局部均值图像去噪方法的基础上，提出一种块匹配三维协同滤波的去噪方法。该方法结合了图像空域的非局部相似性和频域中的稀疏性，在获取空域相似块的基础上，对多个空域相似块构成的分组在频域中进行三维协同滤波，可以在有效抑制图像噪声的同时尽可能的保留图像中的细节信息，是目前公认的性能较好的图像去噪方法。然而，对于单帧的图像去噪问题而言，在未知噪声统计特性的前提下从噪声图像中恢复出原始的无噪声目标图像，必须在一定的约束条件之下才可以获得近似的最优解。已有的方法大都将对图像中噪声的统计特性进行理想的先验性假设作为求解近似最优解的约束条件，由于与实际数字图像采集设备所获取的自然场景图像中的噪声统计特性存在较大差异，因而去噪效果并不理想。
[0004] 相比于单帧图像，使用针对同一场景的多帧图像可以提供更为丰富的时域信息，进而为图像去噪研究提供了新的思路。对于低光照环境下基于多帧的图像去噪，为了避免由于长时间曝光所导致的各帧图像中的运动模糊，要求所采集的针对同一场景的多帧图像必须是在较短的曝光时间内获取的，而较短的曝光时间势必会由于进光量减少导致图像中的噪声强度显著增强。这种情况下，如何有效的利用多帧图像中的空域和时域相似性来抑制噪声的影响，是一项非常复杂而又富有挑战性的工作。此外，由于拍摄者手部的抖动所造成的多帧图像之间存在相对位移，以及场景中的局部运动所导致的图像局部模糊，都为使用多帧图像去噪设置了障碍。直接将已有的用于单帧的图像去噪方法扩展到多帧图像去噪，由于没有充分利用多帧图像所提供的丰富的时域信息，同时不能很好的解决上述的多帧图像之间的相对位移和场景中局部运动对图像质量造成的影响，因而很难获得理想的去噪效果。
发明内容
[0005] 本发明的目的在于针对上述已有技术的不足，提出一种基于空域和时域联合滤波的多帧数字图像去噪方法，在抑制噪声的同时消除低光照环境下多帧图像之间的相对位移和场景中局部运动所导致的图像模糊，提高图像的去噪效果，增强图像的主观质量。
[0006] 本发明目的技术方案是：通过对低光照环境下采集的多帧图像进行全局和局部配准，分析和判别场景中的局部运动区域，利用局部运动区域的清晰度测度和多帧图像的时域相关性消除运动模糊；结合多帧图像的空域和时域相似性对图像进行空域和时域的融合滤波，抑制图像中的噪声；结合场景的亮度特性对去噪图像进行色度校正和进行对比度增强，实现多帧图像去噪。其实现步骤包括如下：
[0007] (1)输入N帧低光照环境下采集的针对同一场景的自然图像zi,i＝1,2,…,N；
[0008] (2)计算各帧图像的平均水平梯度 和平均垂直梯度 将 与 之和作为图像清晰度测度Gi，选取Gi最大的一帧，作为参考图像zr,1≤r≤N，并将该参考图像的清晰度测度记为Gmax；
[0009] (3)计算各帧图像相对于参考图像zr的全局运动矢量Vi；
[0010] (4)计算各帧图像清晰度测度Gi与参考图像清晰度测度Gmax的比值Ri，将Ri大于阈值0.99所对应的第i帧图像记为候选参考图像zr′，计算各帧图像相对于该候选参考图像zr′的候选全局运动矢量
[0011] (5)计算全局运动矢量Vi和候选全局运动矢量 的测度，并依据该测度从参考图像zr和候选参考图像zr′中选取基准图像 1≤s≤N，并将该基准图像 对应的全局运动矢量作为基准全局运动矢量
[0012] (6)将输入的各帧图像都划分成大小相同的J个搜索块，利用上述基准全局运动矢量 计算各帧图像中每个搜索块相对于基准图像 中对应搜索块的局部运动矢量Vi_j,1≤i≤N且i≠s,1≤j≤J；
[0013] (7)将各帧图像中的每个搜索块都划分为大小相同的Q个融合块，使用步骤(6)中所得局部运动矢量Vi_j，在基准图像 中找出与各帧图像中的每个融合块所对应的基准融合块，选取各帧图像中与同一基准融合块对应的融合块构成时域融合块组；
[0014] (8)分别计算每个时域融合块组中每个融合块与其对应基准融合块的像素差值，并将这些像素差值的平均值作为该时域融合块组中每个融合块与其对应基准融合块的时域相似性测度δi_j_k,1≤i≤N且i≠s,1≤j≤J,1≤k≤Q；
[0015] (9)依据步骤(8)中所得的时域相似性测度δi_j_k，计算时域滤波权值ωi_j_k：
[0016]
[0017] (10)用所述时域滤波权值ωi_j_k对每个时域融合块组所对应的基准融合块进行时域加权融合滤波，得到时域融合滤波图像zfu；
[0018] (11)依据步骤(8)中所得各个时域融合块组中每个融合块与其对应基准融合块的像素差值，统计像素差值大于阈值24的像素个数占融合块总像素个数的百分比，计算这些百分比的标准差σj,k，若σj,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域；
[0019] (12)对步骤(11)中所判定的局部运动区域，依据时域相似性测度δi_j_k，按照下式计算局部运动区域的时域滤波权值ω′i_j_k：
[0020]
[0021] (13)用局部运动区域的时域滤波权值ω′i_j_k对所判定的局部运动区域所对应的基准融合块进行加权融合滤波，使用该融合结果覆盖时域融合滤波图像zfu中对应位置的融合结果，得到最终的时域融合图像 并对该时域融合图像 进行单帧空域非局部均值滤波，得到空域滤波图像zf；
[0022] (14)计算空域滤波图像zf的平均亮度 根据该平均亮度 对空域滤波图像zf进行亮度校正和色度校正，得到亮度和色度校正图像zc，对zc再进行局部对比度增强，得到最终的去噪结果图像zout。
[0023] 本发明与现有方法相比具有以下优点：
[0024] 1.本发明由于在对多帧图像进行全局和局部配准的过程中，依据图像清晰度测度和全局运动矢量测度进行基准图像的选取，并使用该基准图像作为参考进行全局和局部配准，提高了配准的精度，为后续时域融合滤波提供了准确的时域相关信息，有助于提升图像去噪的效果；
[0025] 2.本发明由于在对多帧图像进行时域融合的过程中，对场景中的局部运动区域进行了判定，并对局部运动区域采用了与场景中其他区域不同的时域融合方法，可在保证时域去噪效果的同时消除场景中由于局部运动所引起的图像局部模糊，提高了去噪图像的主观质量；
[0026] 3.本发明由于对空域滤波图像进行了亮度和色度校正，降低了光照不足所引起的图像亮度和色度偏差，并对亮度和色度校正后图像进行局部对比度增强，消除了空时域滤波所导致的对比度下降，进一步提升了去噪图像的主观质量。
附图说明
[0027] 图1是本发明的流程图；
[0028] 图2是本发明中使用的图像下采样示意图；
[0029] 图3是本发明中进行中空域非局部均值滤波权值搜索示意图；
[0030] 图4是本发明测试使用的低光照环境下采集的6帧自然图像；
[0031] 图5是现有非局部均值滤波去噪算法的结果图像；
[0032] 图6是现有块匹配三维协同滤波去噪算法的结果图像；
[0033] 图7是本发明仿真得到的去噪结果图像；
[0034] 图8是输入的第1帧自然图像的局部原始分辨率显示。
[0035] 图9是现有非局部均值滤波去噪算法的结果图像的局部原始分辨率显示。
[0036] 图10是现有块匹配三维协同滤波去噪算法的结果图像的局部原始分辨率显示。
[0037] 图11是本发明仿真所得去噪结果图像的局部原始分辨率显示。
具体实施方式
[0038] 下面结合附图对本发明的实施例做详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和过程，但本发明的保护范围不限于下述的实施例。
[0039] 参照图1，本发明的实现步骤如下：
[0040] 步骤1：输入多帧图像。
[0041] (1a)输入低光照环境下采集的针对同一场景的N帧自然图像zi,i＝1,2,…,N，本实施例中取N＝6，但不局限于N＝6的情况；
[0042] (1b)输入各帧图像的分辨率为W×H像素，本实施例中W＝3264，H＝2448，但不局限于输入图像分辨率为3264×2448的情况；
[0043] (1c)本实施例中，输入图像的格式为YUV，对于其他图像格式，例如RGB、YCbCr、HSL等均可转化为YUV格式后进行输入。
[0044] 步骤2：选取参考图像。
[0045] (2a)对步骤1中输入的6帧图像zi,i＝1,2,…,6的亮度分量Y进行2级下采样：
[0046] 图2给出本步骤进行下采样的原理示意图，图2中a,b,c,d表示待采样图像中大小为2×2的像素块，x表示下采样图像中由a,b,c,d下采样所得到的像素，x的计算公式如下：
[0047]
[0048] 依照图2，第1级下采样过程如下：
[0049] 输入图像zi的亮度分量Y，对该亮度分量Y中每个大小为2×2的像素块进行下采样，使用式1)计算该像素块的下采样值，所有2×2像素块的下采样值构成该亮度分量Y的第1级下采样图像z1_i；
[0050] 第2级的下采样过程如下：
[0051] 输入图像zi亮度分量Y的第1级下采样图像z1_i，对该下采样图像z1_i中每个大小为2×2的像素块使用式1)计算该像素块的下采样值，所有2×2像素块的下采样值构成输入图像zi亮度分量Y的第2级下采样图像z2_i；
[0052] 所述对输入图像zi的亮度分量Y的下采样不限于本实施例中所使用的方法，也可以采用已有的其他下采样方法，例如从每2×2个像素块中选择1个像素作为下采样图像中对应位置的像素；
[0053] (2b)分别计算下采样图像z2_i中每个像素的水平梯度和垂直梯度，以像素坐标为(x,y)的像素为例，其水平梯度Gi_h(x,y)和垂直梯度Gi_v(x,y)的计算公式如下：
[0054]
[0055] 其中，A(x,y)表示下采样图像z2_i中以像素坐标(x,y)为中心大小为3×3的图像块，Sh和Sv分别表示3×3大小的水平Sobel算子Sh和垂直Sobel算子Sv：
[0056]
[0057] “*”表示卷积运算；
[0058] 所述计算水平梯度和垂直梯度的Sobel算子不限于本实施例中所使用的Sobel算子，也可以使用其他形式的Sobel算子；
[0059] (2c)分别计算各帧的下采样图像z2_i的平均水平梯度 和垂直梯度[0060]
[0061] 其中，W2_i和H2_i分别表示下采样图像z2_i的宽度和高度；
[0062] (2d)分别计算各帧的下采样图像z2_i的图像清晰度测度Gi：
[0063]
[0064] 选取下采样图像z2_i中Gi最大的一帧，选取其对应的输入图像zi作为参考图像zr,1≤r≤N，并记该参考图像zr的清晰度测度为Gmax；
[0065] 所述计算图像清晰度测度Gi不限于本实施例中所使用平均水平梯度 和垂直梯度 之和，也可以使用其他形式的图像清晰度测度Gi，例如使用平均水平梯度 和垂直梯度 的平方和等。
[0066] 步骤3：计算全局运动矢量Vi。
[0067] (3a)对步骤1中所述的输入图像zi,i＝1,2,…,6的亮度分量Y，利用与步骤(2a)相同的方法进行M＝4级下采样，得到输入图像zi亮度分量Y的下采样图像zm_i，以及参考图像zr亮度分量Y的下采样图像zm_r,1≤m≤4,1≤r≤6；
[0068] 所述的下采样级数M不限于本实施例中所选取的M＝4，M的取值与输入图像的大小有关，可依照第M级下采样图像的分辨率不大于256×256像素确定下采样级数M；
[0069] (3b)以第4级参考图像的下采样图像z4_r为参考，使用第4级下采样图像z4_i相对于该参考的中心位置进行水平和垂直位移，并计算位移后下采样图像z4_i与参考图像的下采样图像z4_r的均方误差D4_i：
[0070]
[0071] -W4_i≤u4_i≤W4_i，-H4_i≤v4_i≤H4_i，
[0072] 其中，W4_i和H4_i分别表示下采样图像z4_i的宽度与高度，p(x,y)表示参考图像的下采样图像z4_r中坐标为(x,y)的像素，q(x,y)表示下采样图像z4_i中坐标为(x,y)的像素，u4_i表示水平位移，v4_i表示垂直位移；
[0073] (3c)选取使均方误差D4_i最小的水平位移u4_i和垂直位移v4_i，分别作为全局运动矢量V4_i的水平分量V4_i_h和垂直分量V4_i_v；
[0074] 所述运动矢量选取准则，不限于本实施例中所采用的均方误差最小准则，也可以采用其他准则，例如绝对差之和(SAD)最小准则等；
[0075] (3d)以步骤(3c)中全局运动矢量V4_i所指向的第3级参考图像的下采样图像z3_r中的像素位置为中心，使用第3级下采样图像z3_i在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z3_i相对于参考图像的下采样图像z3_r的均方误差D3_i：
[0076]
[0077] -W3_i≤u3_i≤W3_i，-H3_i≤v3_i≤H3_i，
[0078] 其中，W3_i和H3_i分别表示下采样图像z3_i的宽度与高度，p(x,y)表示参考图像的下采样图像z3_r中坐标为(x,y)的像素，q(x,y)表示下采样图像z3_i中坐标为(x,y)的像素，u3_i表示水平位移，v3_i表示垂直位移；
[0079] 选取使均方误差D3_i最小的水平位移u3_i和垂直位移v3_i，对全局运动矢量V4_i的水平分量V4_i_h和垂直分量V4_i_v进行修正，得到全局运动矢量V3_i；
[0080] 所述的全局运动矢量修正范围不限于本实施例中所采用的5×5像素大小的邻域，也可以采用不同尺度的邻域，例如7×7像素大小的邻域等；
[0081] (3e)以步骤(3d)中全局运动矢量V3_i所指向的第2级参考图像的下采样图像z2_r中的像素位置为中心，使用第2级下采样图像z2_i在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z2_i相对于参考图像的下采样图像z2_r的均方误差D2_i：
[0082]
[0083] -W2_i≤u2_i≤W2_i，-H2_i≤v2_i≤H2_i，
[0084] 其中，W2_i和H2_i分别表示下采样图像z2_i的宽度与高度，p(x,y)表示参考图像的下采样图像z2_r中坐标为(x,y)的像素，q(x,y)表示下采样图像z2_i中坐标为(x,y)的像素，u2_i表示水平位移，v2_i表示垂直位移；
[0085] 选取使均方误差D2_i最小的水平位移u2_i和垂直位移v2_i，对全局运动矢量V3_i的水平分量V3_i_h和垂直分量V3_i_v进行修正，得到全局运动矢量V2_i；
[0086] (3f)以步骤(3d)中全局运动矢量V2_i所指向的第1级参考图像的下采样图像z1_r中的像素位置为中心，使用第1级下采样图像z1_i在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后下采样图像z1_i相对于参考图像的下采样图像z1_r的均方误差D1_i：
[0087]
[0088] -W1_i≤u1_i≤W1_i，-H1_i≤v1_i≤H1_i，
[0089] 其中，W1_i和H1_i分别表示下采样图像z1_i的宽度与高度，p(x,y)表示参考图像的下采样图像z1_r中坐标为(x,y)的像素，q(x,y)表示下采样图像z1_i中坐标为(x,y)的像素，u1_i表示水平位移，v1_i表示垂直位移；
[0090] 选取使均方误差D1_i最小的水平位移u1_i和垂直位移v1_i，对全局运动矢量V2_i的水平分量V2_i_h和垂直分量V2_i_v进行修正，得到全局运动矢量V1_i；
[0091] (3g)以步骤(3f)中全局运动矢量V1_i所指向的参考图像zr的亮度分量Y中的像素位置为中心，使用输入图像zi的亮度分量Y在该中心5×5像素大小的邻域内进行水平和垂直位移，计算位移后输入图像zi的亮度分量Y相对于参考图像zr的亮度分量Y的均方误差Di：
[0092]
[0093] -W≤ui≤W，-H≤vi≤H，
[0094] 其中，W和H分别表示输入图像zi的宽度与高度，p(x,y)表示参考图像zr的亮度分量Y中坐标为(x,y)的像素，q(x,y)表示输入图像zi的亮度分量Y中坐标为(x,y)的像素，ui表示水平位移，vi表示垂直位移；
[0095] 选取使均方误差Di最小的水平位移ui和垂直位移vi，对全局运动矢量V1_i的水平分量V1_i_h和垂直分量V1_i_v进行修正，得到全局运动矢量Vi。
[0096] 步骤4：确定候选参考图像。
[0097] (4a)计算除参考图像zr外，输入各帧图像的清晰度测度Gi与参考图像zr的清晰度测度Gmax之比Ri，若满足：Ri<0.875，则判定该Ri对应的第i帧输入图像zi为模糊图像，并从输入图像中剔除该Ri对应的第i帧图像，保留剩余图像用于后续处理，并记剩余图像的帧数为N′，否则，保留所有输入图像用于后续处理；
[0098] 以下实施例以保留所有输入图像用于后续处理为例，对于剩余图像帧数为N′的情况，可将后续处理中的帧数6替换为N′即可；
[0099] 所述的模糊图像判定准则不限于本实施例所使用的准则，也可以采用其他的判定准则，例如直接使用清晰度测度Gi的大小作为模糊图像的判定准则；
[0100] (4b)对步骤(4a)中所保留的各帧图像，若存在清晰度测度Gi与参考图像zr的清晰度测度Gmax之比Ri满足：Ri>0.99，则选取该Ri对应的第i帧图像作为候选参考图像zr′，否则，选择参考图像zr作为基准图像 1≤s≤N，并选择全局运动矢量Vi作为基准全局运动矢量 并转入步骤6；
[0101] (4c)以候选参考图像zr′为输入图像中的第1帧，参考图像zr为输入图像中的第2帧为例，即r′＝1,r＝2，从步骤(3g)所得的全局运动矢量Vi中，选取候选参考图像zr′相对参考图像zr的全局运动矢量Vr′：
[0102] Vr′＝V1＝[V1_h,V1_v]T，
[0103] 其中，V1_h和V1_v分别表示全局运动矢量V1的水平分量和垂直分量；
[0104] (4d)根据步骤(4c)所得候选参考图像zr′相对参考图像zr的全局运动矢量Vr′，计算参考图像zr相对于候选参考图像zr′的候选全局运动矢量
[0105]
[0106] (4e)计算输入的第i帧相对于候选参考图像zr′的候选全局运动矢量[0107]
[0108] 使用 i＝3,4,5,6和 r＝2构成输入各帧相对于候选参考图像zr′的候选全局运动矢量 i＝1,2,",6且i≠r′。
[0109] 步骤5：选取基准图像和基准全局运动矢量。
[0110] (5a)计算各帧图像相对于参考图像zr的全局运动矢量测度Sr：
[0111]
[0112] 其中，Vi_h表示第i帧输入图像相对于参考图像zr的全局运动矢量Vi的水平分量，Vi_v表示第i帧输入图像相对于参考图像zr的全局运动矢量Vi的垂直分量；
[0113] 所述全局运动矢量测度不限于本实施例中所使用的全局运动矢量的水平分量与垂直分量的绝对值之和，也可以使用其他全局运动矢量测度，例如全局运动矢量的水平分量与垂直分量的平方和等；
[0114] (5b)计算各帧图像相对于候选参考图像zr′的全局运动矢量测度Sr′：
[0115]
[0116] 其中， 表示第i帧输入图像相对于候选参考图像zr′的全局运动矢量 的水平分量， 表示第i帧输入图像相对于候选参考图像zr′的全局运动矢量 的垂直分量；
[0117] (5c)比较全局运动矢量测度Sr和Sr′，若满足条件：Sr
[0118] 步骤6：计算局部运动矢量。
[0119] (6a)将步骤(2a)中所述的1级下采样图像z1_i划分为J＝64个大小相同的搜索块，将基准图像 的1级下采样图像 划分为J＝64个大小相同的基准搜索块；
[0120] 所述的划分搜索块的个J＝64数不限于本实施例中所使用的64，可根据输入图像的分辨率进行合理选择；
[0121] (6b)以基准全局运动矢量 所指向的第j个基准搜索块中的像素位置为参考，以该参考位置为中心，对第i帧输入图像的下采样图像z1_i的第j个搜索块进行水平位移和垂直位移，计算所述第j个搜索块相对于第j个基准搜索块的均方误差D1_i_j：
[0122]
[0123] -7≤u1_i_j≤7，-7≤v1_i_j≤7，
[0124] 其中，WS1和HS1分别表示所述搜索块的宽度与高度，p(x,y)表示基准搜索块中坐标为(x,y)的像素，q(x,y)表示所述搜索块中坐标为(x,y)的像素，u1_i_j表示水平位移，v1_i_j表示垂直位移；
[0125] (6c)选取使均方误差D1_i_j最小的水平位移u1_i_j和垂直位移v1_i_j，对基准全局运动矢量 的水平分量 和垂直分量 进行修正，得到所述搜索块相对于基准搜索块的局部运动矢量V1_i_j：
[0126]
[0127] 其中，V1_i_j_h和V1_i_j_v分别表示局部运动矢量V1_i_j的水平分量和垂直分量；
[0128] (6d)将输入图像zi的亮度分量Y划分为J＝64个大小相同的搜索块，将基准图像 的亮度分量Y划分为J＝64个大小相同的基准搜索块；
[0129] (6e)以局部运动矢量V1_i_j所指向的基准图像 的亮度分量Y的第j个基准搜索块中的像素位置为参考，以该参考位置为中心，对输入图像zi的第j个搜索块进行水平位移和垂直位移，计算所述第j个搜索块相对于第j个基准搜索块的均方误差Di_j：
[0130]
[0131] -5≤ui_j≤5，-5≤vi_j≤5，
[0132] 其中，WS和HS分别表示所述搜索块的宽度与高度，p(x,y)表示所述基准搜索块中坐标为(x,y)的像素，q(x,y)表示所述搜索块中坐标为(x,y)的像素，ui_j表示水平位移，vi_j表示垂直位移；
[0133] (6f)选取使均方误差Di_j最小的水平位移ui_j和垂直位移vi_j，对局部运动矢量V1_i_j的水平分量V1_i_j_h和垂直分量V1_i_j_v进行修正，得到所述搜索块相对于基准搜索块的局部运动矢量Vi_j：T T
[0134] Vi_j＝[Vi_j_h,Vi_j_v]＝[V1_i_j_h+ui_j,V1_i_j_v+vi_j]，
[0135] 其中，Vi_j_h和Vi_j_v分别表示局部运动矢量Vi_j的水平分量和垂直分量。
[0136] 步骤7：选取融合块组。
[0137] (7a)将步骤(6d)中每个搜索块划分为Q＝108个大小相同的融合块，并将步骤(6d)中每个基准搜索块划分为Q＝108个大小相同的基准融合块；
[0138] 所述的划分融合块的个数Q不限于本实施例中所使用的108，可根据输入图像的分辨率进行合理选择；
[0139] (7b)以第k个基准融合块为参考，在各帧输入图像zi中找出使用步骤(6f)中所得局部运动矢量Vi_j所指向的融合块，用所找出的各帧图像中的融合块构成一个融合块组；1≤k≤108；
[0140] (7c)使用步骤(2d)中所述方法，计算融合块组中每个融合块的清晰度测度Gi,j,k，将清晰度测度Gi,j,k小于阈值TBlur的融合块判定为模糊融合块，并将该模糊融合块从所述融合块组中剔除，阈值TBlur的计算公式如下：
[0141] TBlur＝μg-3σg，
[0142] 其中，μg表示融合块组中所有融合块清晰度测度的均值，σg融合块组中所有融合块清晰度测度的方差；
[0143] 以下实施例按照融合块组中不存在模糊融合块的情况为例，对于存在模糊融合块的情况，可在针对融合块组进行处理时，忽略对模糊融合块的处理，并对融合块组中融合块的个数进行相应调整；
[0144] (7d)经步骤(7c)剔除模糊融合块后的融合块组，即为所选取的融合块组。
[0145] 步骤8：时域融合滤波去噪。
[0146] (8a)计算步骤(7d)中所得融合块组中，每个融合块与对应的基准融合块的时域相似性测度δi_j_k：
[0147]
[0148] 其中，Wf和Hf分别表示融合块的宽度与高度，pi_Y(x,y)表示第i个融合块亮度分量Y中坐标为(x,y)的像素值，qY(x,y)表示基准融合块亮度分量Y中坐标为(x,y)的像素值；1≤i≤6且i≠s,1≤j≤64,1≤k≤108；
[0149] (8b)依据步骤(8a)中所得的时域相似性测度δi_j_k，计算时域融合滤波权值ωi_j_k：
[0150]
[0151] (8c)使用步骤(8b)中所得时域融合滤波权值ωi_j_k，计算基准融合块的亮度分量Y的时域融合滤波结果
[0152]
[0153] 其中，cf表示归一化融合权值，
[0154] (8d)使用与亮度分量Y相同的时域融合滤波权值ωi_j_k，分别计算基准融合块的红色差分量U的时域融合滤波结果 和绿色差分量V的时域融合滤波结果
[0155]
[0156]
[0157] 其中，pi_U(x,y)和pi_V(x,y)分别表示融合块组中第i个融合块红色差分量U和绿色差分量V中坐标为(x,y)的像素值，qU(x,y)和qV(x,y)分别表示基准融合块红色差分量U和绿色差分量V中坐标为(x,y)的像素值；
[0158] (8e)使用所述基准融合块的亮度分量Y的时域融合滤波结果 红色差分量U的时域融合滤波结果 和绿色差分量V的时域融合滤波结果 构成时域融合滤波图像zfu中与基准融合块相对应位置的融合结果，所有基准融合块的时域融合滤波结果共同构成时域融合滤波图像zfu。
[0159] 步骤9：局部运动区域判定。
[0160] (9a)使用步骤(7d)中所得融合块组，计算融合块组中每个融合块与基准融合块的对应位置亮度分量Y的像素差值di_Y(x,y)：
[0161] di_Y(x,y)＝pi_Y(x,y)-qY(x,y)，
[0162] 其中，pi_Y(x,y)表示融合块组中第i个融合块亮度分量Y中坐标为(x,y)的像素值，qY(x,y)表示基准融合块亮度分量Y中坐标为(x,y)的像素值；
[0163] (9b)统计di_Y(x,y)大于阈值24的像素个数占融合块总像素个数的百分比pi_j_k，计算融合块组中该百分比pi_j_k的标准差σj,k；
[0164] (9c)若所述标准差σj,k大于阈值0.09，则判定该时域融合块组及其对应的基准融合块为场景中的局部运动区域。
[0165] 步骤10：局部运动区域时域融合滤波修正。
[0166] (10a)对步骤(9)中所判定的局部运动区域，依据时域相似性测度δi_j_k，按照下式计算局部运动区域的时域融合滤波权值ω′i_j_k：
[0167]
[0168] (10b)用局部运动区域的时域滤波权值ω′i_j_k，计算局部运动区域亮度分量Y的时域融合滤波结果
[0169]
[0170] 其中， 表示局部运动区域对应的基准融合块的亮度分量Y中坐标为(x,y)的像素值，c′f表示局部运动区域的归一化的融合权值，
[0171] (10c)使用与亮度分量Y相同的时域滤波权值ω′i_j_k，计算局部运动区域的红色差分量U的时域融合滤波结果 和绿色差分量V的时域融合滤波结果
[0172]
[0173]
[0174] (10d)使用所述局部运动区域亮度分量Y的时域融合滤波结果 红色差分量U的时域融合滤波结果 和绿色差分量V的时域融合滤波结果 构成局部运动区域的时域融合结果。
[0175] (10e)使用对局部运动区域的融合结果，覆盖时域融合滤波图像zfu中对应区域的融合结果，实现对时域融合滤波图像zfu中局部运动区域的修正，得到最终的时域融合图像[0176] 步骤11：空域非局部均值滤波去噪。
[0177] (11a)采用步骤(2a)中所述的下采样方法，对步骤(10e)所得时域融合图像 的亮度分量Y进行1级下采样，得到该时域融合图像 下采样图像Y1_fu；
[0178] (11b)依照非局部均值滤波算法，设定所述下采样图像Y1_fu中的相似块大小为5×5像素，搜索窗的范围为7×7像素；
[0179] (11c)以所述下采样图像Y1_fu中的像素p为例，按照图3所示方式进行遍历搜索，计算像素p的非局部均值滤波权值；
[0180] 图3中，A(p)表示以像素p为中心的相似块，B(p)表示以像素p为中心的搜索窗，A(q)表示以像素q为中心的搜索窗B(p)中的相似块，对B(p)内的像素按照图3中箭头方向进行逐像素遍历，计算像素p的非局部均值滤波权值w(p,q)：
[0181]
[0182] 其中，h表示滤波强度参数，对步骤9中所判定的局部运动区域，h的取值为20，对于其他区域，h的取值为10， 表示相似块A(p)和A(q)的欧式距离；
[0183] (11d)以时域融合图像 中与像素p对应的像素 为例，使用像素p的非局部均值滤波权值w(p,q)，计算像素p′亮度分量Y的空域加权滤波值
[0184]
[0185] 其中，p′Y表示像素p′的亮度分量，W(p)表示归一化滤波权值，按照如下公式计算：
[0186] W(p)＝∑q∈B(q)w(p,q)；
[0187] (11)使用像素p的非局部均值滤波权值w(p,q)，分别计算像素p′的红色差分量U的空域加权滤波值 和绿色差分量V的空域加权滤波值
[0188]
[0189]
[0190] 其中，p′U表示像素p′的红色差分量，p′ V表示像素的绿色差分量；
[0191] (11e)使用像素p′亮度分量Y的空域加权滤波值 红色差分量U的空域加权滤波值 以及绿色差分量V的空域加权滤波值 构成空域滤波图像zf中与像素p′相对应的非局部均值滤波结果，由所有像素的非局部均值滤波结果共同构成空域滤波图像zf。
[0192] 步骤12：亮度校正和色度校正。
[0193] (12a)计算步骤(11e)中所得的空域滤波图像zf的平均亮度
[0194] (12b)根据空域滤波图像zf的平均亮度 按照如下公式生成亮度校正曲线：
[0195]
[0196] 其中，x表示亮度校正曲线的水平坐标，yl表示平均亮度zf落入不同范围时亮度校正曲线的垂直坐标，l＝1,2,3；
[0197] (12c)根据空域滤波图像zf的平均亮度 选择相应亮度校正曲线，分别计算空域滤波图像zf中像素p(x,y)亮度分量校正值 红色差分量校正值 和绿色差分量校正值
[0198]
[0199] 其中，pY(x,y)表示像素p(x,y)的亮度分量，pU(x,y)表示像素p(x,y)的红色差分量，pV(x,y)表示像素p(x,y)的绿色差分量，αl表示亮度校正曲线在x＝pY(x,y)点的斜率，其计算公式为：
[0200] (12d)使用空域滤波图像zf所有像素亮度分量、红色差分量以及绿色差分量的校正结果构成校正图像zc。
[0201] 步骤13：对比度增强。
[0202] (13a)对校正图像zc进行高斯滤波，得到校正图像zc的高斯滤波结果图像zg；
[0203] (13b)利用高斯滤波结果图像zg中的像素qg(x,y)，计算校正图像zc中像素q(x,y)的局部对比度增强结果
[0204]
[0205] 其中，β表示对比度增强因子，取值为0.3。
[0206] 至此，获得最终的多帧去噪输出图像zout。
[0207] 以图4所示的低光照环境下采集的6帧自然图像的实验为例，得到的最终去噪图像是图7。
[0208] 本发明的优点由以下仿真实现进一步说明。
[0209] 1.仿真条件
[0210] 本发明采用图4所示的低光照环境下采集的针对同一场景的6帧自然图像作为测试图像，其中图4(a)～图4(f)分别为输入的第1～6帧测试图像，软件平台为Visu al Studio 2010。
[0211] 2.仿真内容与结果：
[0212] 仿真1，使用现有的非局部均值滤波算法对图4所示的6帧图像分别进行非局部均值滤波，对滤波所得的6帧图像进行平均，所得的结果如图5所示。
[0213] 仿真2，使用现有的块匹配三维协同滤波算法对图4所示的6帧图像分别进行滤波，对滤波所得的6帧图像进行平均，所得的结果如图6所示。
[0214] 仿真3，使用本发明方法对图4所示的6帧图像分别进行非局部均值滤波，滤波所得的结果如图7所示。
[0215] 为了便于比较所述实验方法的性能，对图4(a)所示的第一帧输入图像的局部区域进行原始分辨率显示，如图8所示；对图5中与图4(a)对应的局部区域进行原始分辨率显示，如图9所示；对图6中与图4(a)对应的区域进行原始分辨率显示，如图10所示；对图7中与图4(a)对应的局部区域进行原始分辨率显示，如图11所示。
[0216] 从图9与图8的对比可以看出，现有的非局部均值滤波算法对低光照环境下采集的自然图像中噪声的抑制十分有限，并且在图像的边缘存在模糊；
[0217] 从图10与图8的对比可以看出，相比于非局部均值滤波算法，现有的块匹配三维协同滤波算法对图像中噪声的抑制性能有了显著提升，但场景中的细节纹理区域，仍存在模糊；
[0218] 从图11与图8的对比可以看出，相比于现有的非局部均值滤波算法和块匹配三维协同滤波算法，本发明方法能够有效抑制低光照环境下采集图像中的噪声，并且边缘与细节的保持度比较高，在去除噪声的同时可以有效的避免细节模糊。
