CN101964863A - 一种自适应的时空域视频图像降噪方法 - Google Patents

Abstract

一种自适应的时空域视频图像降噪方法，包括以下个步骤：第一步：对于当前帧的图像，要进行空域滤波，以去除噪声；第二步：再把空域滤波后的像素与前一帧的对应像素加权平均，进行时域滤波；第三步：更新滤波参数；在一个起始时刻，对滤波参数进行初始化，利用该滤波参数，对第一个宏块进行滤波；每个宏块包含亮度分量，蓝色色度分量，红色色度分量三个分量，要分别对每个分量进行滤波；滤波分为两步，空域滤波和时域滤波；当前宏块滤波结束后，需要对滤波参数进行更新，以进行下一个宏块的滤波，更新过程分为三步，第一步是利用固定系数计算单个像素的差值，第二步是利用单个像素方差计算宏块方差，第三步则是利用方差计算下一宏块的滤波参数。

Description

一种自适应的时空域视频图像降噪方法 

技术领域

本发明属于数字视频处理领域，特别涉及一种自适应的时空域视频图像降噪方法。 

背景技术

利用计算机和网络技术的优势，在人们的日常生活中，数字图像和视频已经越发普遍，去噪可以从数字图像和视频中去除噪声，增强其可压缩性。噪声通常存在于未被压缩的图像和视频里，其是不需要的或不想要的元素或伪影。噪声可能有很多源，例如，在电子硬件里出现的散粒噪声，热噪声，信道噪声等。噪声的存在不仅降低了视频的可视质量，而且降低了视频编码的效果，并降低了媒体文件的可压缩性。它增加了视频的熵，使得需要更多的比特用来编码图像和视频，同时也降低了进一步提高比特率的运动估计精度。为了提高编码效率和可压缩性，在编码视频之前必须进行视频去噪。视频去噪的目的是尽可能准确地估计真实的图像信号。 

在目前存在的方法当中，有一些是直接从图像降噪方法中扩展过来的，如空域滤波方法[1，2]，小波收缩方法[3，4]，以及偏差分方程方法[5]。这些方法由于没有考虑到时空维度上分辨率及连续性的不同，易产生运动模糊和伪影。鉴于此问题，一些算法在空间滤波之前使用了运动补偿方法[6，7]，在[8]中使用了一种带有运动补偿的自适应的时空双向滤波器。另外还有一些基于统计的方法，如[9]中使用了基于贝叶斯的方法来进行去噪。这些算法都没有考虑与编码器的结合，并且硬件上也不易实现。 

发明内容

所要解决的问题： 

针对以上不足本发明提供了一种有效的去噪滤波的一种自适应的时空域视频图像降噪方法。 

技术方案： 

一种自适应的时空域视频图像降噪方法，包括以下个步骤： 

第一步：对于当前帧的图像，要进行空域滤波，以去除噪声； 

第二步：再把空域滤波后的像素与前一帧的对应像素加权平均，进行时域滤波； 

第三步：更新滤波参数； 

在一个起始时刻，对滤波参数进行初始化，利用该滤波参数，对第一个宏块进行滤波；每个宏块包含亮度分量，蓝色色度分量，红色色度分量三个分量，要分别对每个分量进行滤波；滤波分为两步，空域滤波和时域滤波； 

所有邻域像素都使用未滤波前的像素值；空域滤波的公式为 

$p^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)}{8}+p(x,y)$

令d＝abs(p(x+dx，y+dy)-p(x，y))，则 

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\≤\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s＜d\≤2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

时域滤波分两步，第一步是运动估计，该步骤的目的在于为当前宏块像素寻找最佳的时域匹配点，匹配准则为SAD，其计算公式为 

$\mathrm{SAD}\left(\mathrm{mv}\right)=\underset{i}{\mathrm{\Σ}}\mathrm{abs}(p\left(i\right)-p^{\′}(i+\mathrm{mv}))$

式中p(i)表示当前帧中位置i处的像素值，p′(i+mv)表示参考帧中位置i+mv处的像素值；对于某个范围内的mv，对每一个都计算一个SAD值，最小的SAD值所对应的mv为最佳运动位移，其对应的参考块为当前块的最佳匹配块；如果该方法应用于视频编码器中，运动估计已在编码器中实现，不需额外实现； 

时域滤波的第二步是将当前像素与匹配像素进行加权平均，从而得到最终的滤波像素，其计算公式如下 

$p^{\′\′}(x,y,t)=p^{\′}(x,y,t)+\frac{D(p^{\′}(x,y,t),p^{\′}(x^{\′},y^{\′},t-1),\mathrm{factor}\_t)}{2}$

式中，p′(x′，y′，t-1)代表p′(x，y，t)在运动补偿时对应的参考帧的像素；D函数的定义为： 

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\≤\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s＜d\≤2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

当前宏块滤波结束后，需要对滤波参数进行更新，以进行下一个宏块的滤波，更新过程分为三步，第一步是利用固定系数计算单个像素的差值，第二步是利用单个像素方差计算宏块方差，第三步则是利用方差计算下一宏块的滤波参数； 

在计算单个像素方差步骤中，使用了一个固定的滤波系数来，计算公式为 

$d^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),20)}{8}$

计算宏块方差使用了单个像素差值的绝对值，计算公式为 

$\mathrm{var}=\frac{\mathrm{\Σ\Σabs}\left(d^{\′}(x,y)\right)}{\mathrm{pixelcount}}$

最后通过var更新下一个宏块的滤波参数 

factor_s＝1.0+luma_var*luma_var/1.4 

factor_t＝factor_s/2 

有益效果： 

本发明提供了一种有效的去噪滤波，用来自适应性地降低在同质区和高纹理区内的空间噪声，而又不会牺牲图片细节，也不会需要太繁重的计算能力，从而使其能够应用于实时处理。 

附图说明

图1为一种自适应的时空域视频图像降噪方法的原理框图； 

图2为本发明的空域滤波中邻域像素示意图； 

图3为本发明的滤波参数更新流程图； 

图4为本发明集成于编码器的示意图； 

图5-1为针对噪声方差为10的coastguard序列使用与不使用滤波算法的编码率失真对比； 

图5-2为针对噪声方差为10的foreman序列使用与不使用滤波算法的编码率失真对比； 

图6-1为针对噪声方差为15的coastguard序列使用与不使用滤波算法的编码率失真对比； 

图6-2为针对噪声方差为15的foreman序列使用与不使用滤波算法的编码率失真对比； 

具体实施方式

本发明主要包括以下三个步骤： 

第一步：对于当前帧的图像，要进行空域滤波，以去除噪声 

第二步：再把空域滤波后的像素与前一帧的对应像素加权平均，进行时域滤波。 

第三步：更新滤波参数。 

如图1所示，以宏块为单位，不断迭代实现。在一个起始时刻，对滤波参数进行初始化，利用该滤波参数，对第一个宏块进行滤波。每个宏块包含亮度分量，蓝色色度分量，红色色度分量三个分量，要分别对每个分量进行滤波。滤波分为两步，空域滤波和时域滤波。 

空域滤波示意图如图2所示，中间方块代表当前要滤波的像素点，四周八个方块代表滤波时所使用的邻域像素。所有邻域像素都使用未滤波前的像素值。空域滤波的公式为 

$p^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)}{8}+p(x,y)$

令d＝abs(p(x+dx，y+dy)-p(x，y))，则 

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\≤\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s＜d\≤2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

时域滤波分两步，第一步是运动估计，该步骤的目的在于为当前宏块像素寻找最佳的时域匹配点，匹配准则为SAD，其计算公式为 

$\mathrm{SAD}\left(\mathrm{mv}\right)=\underset{i}{\mathrm{\Σ}}\mathrm{abs}(p\left(i\right)-p^{\′}(i+\mathrm{mv}))$

式中p(i)表示当前帧中位置i处的像素值，p′(i+mv)表示参考帧中位置i+mv处的像素值。对于某个范围内的mv，对每一个都计算一个SAD值，最小的SAD值所对应的mv为最佳运动位移，其对应的参考块为当前块的最佳匹配块。如果该方法应用于视频编码器中，运动估计已在编码器中实现，不需额外实现。 

时域滤波的第二步是将当前像素与匹配像素进行加权平均，从而得到最终的滤波像素，其计算公式如下 

$p^{\′\′}(x,y,t)=p^{\′}(x,y,t)+\frac{D(p^{\′}(x,y,t),p^{\′}(x^{\′},y^{\′},t-1),\mathrm{factor}\_t)}{2}$

式中，p′(x′，y′t-1)代表p′(x，y，t)在运动补偿时对应的参考帧的像素。D函数的定义与空域滤波中相同。 

当前宏块滤波结束后，需要对滤波参数进行更新，以进行下一个宏块的滤波，更新流程如图3所示，更新过程分为三步，第一步是利用固定系数计算单个像素的差值，第二步是利用单个像素方差计算宏块方差，第三步则是利用方差计算下一宏块的滤波参数。 

在计算单个像素方差步骤中，使用了一个固定的滤波系数来，计算公式为 

$d^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),20)}{8}$

计算宏块方差使用了单个像素差值的绝对值，计算公式为 

$\mathrm{var}=\frac{\mathrm{\Σ\Σabs}\left(d^{\′}(x,y)\right)}{\mathrm{pixelcount}}$

最后通过var更新下一个宏块的滤波参数 

factor_s＝1.0+luma_var*luma_var/1.4 

factor_t＝factor_s/2 

下面以本发明应用于AVC编码器的场景为例，说明本发明的具体实施方式。 

图4是一个宏块的编码过程，空心框中为标准编码器中的编码流程，实心框中表示了本发明的滤波过程，该图描述了本发明如何集成于标准的视频编码器中。 

当一个宏块的数据被送进来时，首先由空域滤波模块对其进行降噪处理。每个宏块包含Y，Cb，Cr三个分量，要分别对每个分量进行滤波。 

空域滤波示意图如图2所示，中间方块代表当前要滤波的像素点，四周八个方块代表滤波时所使用的邻域像素。所有邻域像素都使用未滤波前的像素值。空域滤波的公式为 

$p^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)}{8}+p(x,y)$

令d＝abs(p(x+dx，y+dy)-p(x，y))，则 

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\≤\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s＜d\≤2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

在起始宏块中，factor_s设为一个固定值，在以后的宏块中，其值由参数更新模块提供。 

接着，编码器以空域滤波后的数据进行帧内预测和运动估计，得到预测数据。运动估计得到的预测数据用于时域滤波，如图所示。将当前像素与预测值进行加权平均，从而得到最终的滤波像素，其计算公式如下 

$p^{\′\′}(x,y,t)=p^{\′}(x,y,t)+\frac{D(p^{\′}(x,y,t),p^{\′}(x^{\′},y^{\′},t-1),\mathrm{factor}\_t)}{2}$

式中，p′(x′，y′，t-1)代表p′(x，y，t)在运动补偿时对应的参考帧的像素。D函数的定义与空域滤波中一样。factor_t为factor_s的1/2。 

编码器在模式决策后，将时域滤波后的值与最终预测值相减得到残差，然后进行变换，量化，熵编码，得到输出码流。 

在编码器编码的同时，参数更新模块也利用原始数据计算下一个宏去噪时所使用的参数。更新流程图如图3所示，更新过程分为三步，第一步是利用固定系数计算单个像素的差值，第二步是利用单个像素方差计算宏块方差，第三步则是利用方差计算下一宏块的滤波参数。 

在计算单个像素方差步骤中，使用了一个固定的滤波系数来，计算公式为 

$d^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),20)}{8}$

计算宏块方差使用了单个像素差值的绝对值，对于亮度块，计算公式为 

$\mathrm{var}=\frac{\mathrm{\Σ\Σabs}\left(d^{\′}(x,y)\right)}{256}$

对于色度块，计算公式为 

$\mathrm{var}=\frac{\mathrm{\Σ\Σabs}\left(d^{\′}(x,y)\right)}{64}$

最后通过var更新下一个宏块的滤波参数 

factor_s＝1.0+luma_var*luma_var/1.4 

factor_t＝factor_s/2 

性能测试 

1、实验条件：选择两个CIF序列，在其中加入噪声方差为10的高斯噪声，在滤波和不滤波的条件下使用AVS在不同QP下各编码100帧图像，将编码图像与未加噪图像进行比较计算PSNR，从而得到率失真曲线。QP分别为8，12，16，20，24，27。 

图5-1、图5-2的英文名称的含义如下： 

Unfilter：不使用时空滤波。 

Filter：使用时空滤波。 

2、实验条件：选择两个CIF序列，在其中加入噪声方差为15的高斯噪声，在滤波和不滤波的条件下使用AVS在不同QP下各编码100帧图像，将编码图像与未加噪图像进行比较计算PSNR，从而得到率失真曲线。QP分别为8，12，16，20，24，27。 

图6-1、图6-2的英文名称的含义如下： 

Unfilter：不使用时空滤波。 

Filter：使用时空滤波。 

[1]J.Boulanger，C.Kervrann，and P.Bouthemy.Space-time adaptation for patch-based imagesequence restoration.Technical report，2006. 

[2]A.Buades，B.Coll，and J.-M.Morel.Denoising image sequences does not require motionestimation.Technical report，Preprint CMLA 2005-18，2005. 

[3]N.Rajpoot，Z.Yao，and R.Wilson.Adaptive wavelet restoration of noisy video sequences.IEEE International Conference on Image Processing(ICIP)，2004. 

[4]I.W.Selesnick and K.Y.Li.Video denoising using 2d and 3d dual-tree complex wavelettransforms.International Conference on Acoustics，Speech，and Signal Processing(ICASSP)，2004. 

[5]S.H.Lee and M.G.Kang.Spatio-temporal video filtering algorithmbased on 3-d anisotropicdiffusion equation.In The IEEE International Conference on Image Processing(ICIP)，1998. 

[6]R.Dugad and N.Ahuja.Video denoising by combining kalman and wiener estimates.InIEEE International Conference on Image Processing，1999. 

[7]V.Zlokolica，A.Pizurica，and W.Philips.Wavelet-domain video denoising based onreliability measures.IEEE Transactions on Circuits and Systems for Video Technology，16(8)：993-1007，2006. 

[8]E.P.Bennett and L.McMilla.Video enhancement using per-pixel virtual exposures.InACMSIGGRAPH，2005. 

[9]A.Kokaram and S.J.Godsill.Mcmc for joint noise reduction and missing data treatementin degraded video.IEEE Transactions on Signal Processing，Special Issue on MCMC，50：189-205，2002. 

Claims (1)

1.一种自适应的时空域视频图像降噪方法，其特征在于：包括以下个步骤：

第一步：对于当前帧的图像，要进行空域滤波，以去除噪声；

第二步：再把空域滤波后的像素与前一帧的对应像素加权平均，进行时域滤波；

第三步：更新滤波参数；

在一个起始时刻，对滤波参数进行初始化，利用该滤波参数，对第一个宏块进行滤波；每个宏块包含亮度分量，蓝色色度分量，红色色度分量三个分量，要分别对每个分量进行滤波；滤波分为两步，空域滤波和时域滤波；

所有邻域像素都使用未滤波前的像素值；空域滤波的公式为

$p^{\′}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\Σ}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\Σ}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)}{8}+p(x,y)$

令d＝abs(p(x+dx，y+dy)-p(x，y))，则

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\&le;\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s<d\&le;2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

时域滤波分两步，第一步是运动估计，该步骤的目的在于为当前宏块像素寻找最佳的时域匹配点，匹配准则为SAD，其计算公式为

$\mathrm{SAD}\left(\mathrm{mv}\right)=\underset{i}{\mathrm{\&Sigma;}}\mathrm{abs}(p\left(i\right)-p^{\&prime;}(i+\mathrm{mv}))$

式中p(i)表示当前帧中位置i处的像素值，p′(i+mv)表示参考帧中位置i+mv处的像素值；对于某个范围内的mv，对每一个都计算一个SAD值，最小的SAD值所对应的mv为最佳运动位移，其对应的参考块为当前块的最佳匹配块；如果该方法应用于视频编码器中，运动估计已在编码器中实现，不需额外实现；

时域滤波的第二步是将当前像素与匹配像素进行加权平均，从而得到最终的滤波像素，其计算公式如下

$p^{\&prime;\&prime;}(x,y,t)=p^{\&prime;}(x,y,t)+\frac{D(p^{\&prime;}(x,y,t),p^{\&prime;}(x^{\&prime;},y^{\&prime;},t-1),\mathrm{factor}\_t)}{2}$

式中，p′(x′，y′，t-1)代表p′(x，y，t)在运动补偿时对应的参考帧的像素；D函数的定义为：

$D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),\mathrm{factor}\_s)=\left\{\begin{array}{cc}d& d\&le;\mathrm{factor}\_s\\ d>>1& \mathrm{factor}\_s<d\&le;2*\mathrm{factor}\_s\\ d>>4& d>2*\mathrm{factor}\_s\end{array}\right.$

当前宏块滤波结束后，需要对滤波参数进行更新，以进行下一个宏块的滤波，更新过程分为三步，第一步是利用固定系数计算单个像素的差值，第二步是利用单个像素方差计算宏块方差，第三步则是利用方差计算下一宏块的滤波参数；

在计算单个像素方差步骤中，使用了一个固定的滤波系数来，计算公式为

$d^{\&prime;}(x,y)=\frac{\underset{\mathrm{dy}=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{\mathrm{dx}=-1}{\overset{1}{\mathrm{\&Sigma;}}}D(p(x+\mathrm{dx},y+\mathrm{dy}),p(x,y),20)}{8}$

计算宏块方差使用了单个像素差值的绝对值，计算公式为

$\mathrm{var}=\frac{\mathrm{\&Sigma;\&Sigma;abs}\left(d^{\&prime;}(x,y)\right)}{\mathrm{pixel\; count}}$

最后通过var更新下一个宏块的滤波参数

factor_s＝1.0+luma_var*luma_var/1.4

factor_t＝factor_s/2

Images