图像降噪算法——时域降噪算法

图像降噪算法——时域降噪算法

最近在工作上接触到了时域降噪相关的算法，这里进行一个简单的总结。我们可以这样理解降噪算法，降噪的过程其实就是通过观测减小噪声方差的过程，对于单帧降噪，做法可以是通过在图片上寻找类似的像素，通过加权平均来获得减少噪声的方差，最经典的就是非局部均值滤波算法，而对于时域降噪，这个概念就更加直接，例如我在第一帧我看到这个像素值是200，第二帧我看到这个像素值是100，第三帧这个像素值是180，第四帧这个像素是190，那么这个像素就就更有可能是180而不是100。

在时域降噪算法中，最重要的两个步骤分别是对齐和融合，下面我分享两篇论文来对这个过程进行解析.

1.《MeshFLow Video Denoising》

这是2016年的一篇paper，这篇论文实现的Mesh FLow是基于Homography FLow实现的，Homography FLow的实现可以参考论文《Fast Burst Image Denoising》，所谓Homograhpy Flow就是将图像帧划分为一个个节点，然后通过求解单应矩阵描述不同帧节点与节点之间的变换关系，然后通过单应矩阵来对齐像素，数学原理和SLAM的前端是类似的。Mesh Flow的流程图如下：

根据Mesh FLow的运动模型对齐步骤如下：
（1）预处理：输入图片后将图片进行网格划分，把每网格的角点当做节点；
（2）特征点追踪：通过KLT光流法对前后帧之间的特征点进行追踪；
（3）运动求解：首先通过追踪到的特征求取前后帧的全局单应矩阵，通过全局单应矩阵将前一阵节点的位置投影到当前帧，然后将当前帧上距离节点小于一定阈值的特征点的运动向量存储该到节点的运动向量Buffer，然后对运动向量的Buffer进行中值滤波，将中值滤波后的结果当做节点的运动向量；
（4）运动滤波：对节点和距离该节点小于一定阈值的节点的运动向量放到Buffer中进行中值滤波，将中值滤波的结果当做该节点的运动向量；
（5）构建映射：遍历所有网格，利用网格的四个节点的运动向量求解该网格运动的单应矩阵，然后利用该单应矩阵构建网格内若有像素的映射图，即根据给像素的位置查询映射图就可以知道该像素在下一帧的位置；
（6）像素对齐：根据映射图对齐像素位置，由于该算法中是利用多帧对齐，所以中间还有一步是将映射图叠加，构建不同帧之间的累加运动的映射图。

像素对齐之后就是进行融合，该算法融合的步骤非常简单，就是对齐后不同帧之间像素灰度值差小于一定阈值的像素进行平均作为当前帧的灰度值，这种简单的融合在灰度变化明显的区域容易造成”坏点“。

该算法较为简单，而且有开源代码，开源代码中是将所有图片全部读回来再一起处理，改成实时读图片的接口也不难，另外个人觉得可以优化的几点是：
（1）在特征点追踪是不一定要使用光流法，KLT光流的鲁棒性有限，可以试试ORB特征点说不定在噪声环境下能有更好的下过；
（2）融合可以结合空域滤波进行优化，例如某一个像素在对齐之后前后帧的像素都不满足像素差阈值可以直接对该像素进行空域滤波，这样可以减少”坏点“的存在；
（3）算法是多帧融合，相机运动较快时”鬼影“现象是不可避免的。

2. 《Real-Time Video Denoising On Mobile Phones》

这是2018年的一篇Paper，这篇论文是基于DIS稠密光流实现的，是Google Pixel 2相机上硬件实现的算法，相较上一篇论文实现更复杂，该算法流程图如下：
该算法的对齐部分是根据DIS稠密光流算法，这里先介绍DIS稠密光流，DIS稠密光流参考论文《Fast Optical Flow using Dense Inverse Search》。
DIS稠密光流法首先将图像进行网格划分成图像块，然后构建金字塔，然后对图像金字塔逐层进行图像块的光流追踪，每一层光流由上一层的光流进行初始化。下面详细介绍下DIS稠密光流中使用的逆向组合算法。对于参考帧$I_t$中的一个图像块$T$,设其尺寸为${\theta }_{ps}\times {\theta }_{ps}$，坐标位置为$\mathbf{x}=(x,y{)}^T$，我们需要使用梯度下降法在帧$I_{t+1}$寻找到最佳匹配的尺寸同样为${\theta }_{ps}\times {\theta }_{ps}$的图像块，实际上就是寻找一个运动向量$\mathbf{u}=(u,v)$使得前后帧的图像块平方差最小：$$\mathbf{u}={\mathrm{argmin}}_{{\mathbf{u}}^{\prime }}\sum _x{\left[I_{t+1}\left(\mathbf{x}+{\mathbf{u}}^{\prime }\right)-T(\mathbf{x})\right]}^2$$例如LK光流就是将其作为一个非线性优化问题，运动向量$\mathbf{u}=(u,v)$作为优化变量，使用迭代法进行求解。在DIS稠密光流中对该方法进行了修改。DIS稠密光流将优化的目标修改为当前运动估计$\mathbf{u}$基础上的一个更新向量$\Delta \mathbf{u}$：$$\Delta \mathbf{u}={\mathrm{argmin}}_{\Delta {\mathbf{u}}^{\prime }}\sum _x{\left[I_{t+1}\left(\mathbf{x}+\mathbf{u}+\Delta {\mathbf{u}}^{\prime }\right)-T(\mathbf{x})\right]}^2$$进一步地将上式的目标损失函数修改为${\sum }_x{\left[T(\mathbf{x}-\Delta \mathbf{u})-I_{t+1}(\mathbf{x}+\mathbf{u})\right]}^2$再进行求解，修改成这样有什么好处呢？这样的修改方便我们实现LK光流的逆向组成算法，为了进一步理解逆向组成算法的流程，我们先对目标损失函数进行进一步分解：
我们先将上式改写为：$$\sum _x{\left[T(\mathbf{W}(\mathbf{x};\Delta \mathbf{u}))-I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))\right]}^2$$其中$\mathbf{W}(\mathbf{x};\mathbf{u})$描述的是一种映射关系$$\mathbf{W}(\mathbf{x};\mathbf{u})=\left(\begin{array}{l}x-u\\ y-v\end{array}\right)$$我们对上式进行一阶泰勒展开有：$$\sum _x{\left[T(\mathbf{W}(\mathbf{x};0))+\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\Delta \mathbf{u}-I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))\right]}^2$$其中$\nabla T$为图像块$T$梯度图像，$\partial \mathbf{W}/\partial \mathbf{u}$为$\mathbf{u}=0$时的雅克比。为对上式求导，并令导数等于0，得到下式：$$2\sum _x{\left[\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\right]}^T\left[T(\mathbf{W}(\mathbf{x};0))+\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\Delta \mathbf{u}-I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))\right]=0$$上式对$\Delta \mathbf{u}$进行求解，得到$\Delta \mathbf{u}$的解析解为：$$\Delta \mathbf{u}=H^{-1}\sum _x{\left[\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\right]}^T[T(\mathbf{W}(\mathbf{x};0))-I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))]$$其中$$H=\sum _x{\left[\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\right]}^T\left[\nabla T\frac{\partial \mathbf{W}}{\partial \mathbf{u}}\right]$$由于$\partial \mathbf{W}/\partial \mathbf{u}$为$\mathbf{u}=0$时的雅克比，$\nabla T$为梯度图像块$T$的梯度图像，这两者都是定值，因此矩阵$H$也是定值，这样就免去了我们在迭代优化$\Delta \mathbf{u}$时反复求解$H$矩阵的问题。对上面这段推导有点懵逼的同学可以参考下论文《Lucas-Kanade 20 Years On: A Unifying Framework》，这篇论文将光流法总结为四种方法，按照更新量的不同分为additional和compositional两种，按照更新方式不同分为forward和inverse两种，上面推导的这种inverse ccompositional(逆向组成)算法是其中一种。
在完成图像块的光流追踪之后，当前帧中的像素可能被参考帧中的多个图像块覆盖，我们对所有覆盖该像素的图像块光流进行加权平均就获得该点的光流值：$${\mathbf{U}}_s(\mathbf{x})=\frac{1}{Z}\sum _i^{N_s}\frac{{\lambda }_{i,\mathbf{x}}}{\max \left(1,{\Vert d_i(\mathbf{x})\Vert }_2\right)}\cdot {\mathbf{u}}_i$$其中${\mathbf{u}}_i$为图像块的光流，${\mathbf{U}}_s(\mathbf{x})$为像素$\mathbf{x}$处的光流大小，$d_i(\mathbf{x})=I_{t+1}\left(\mathbf{x}+{\mathbf{u}}_i\right)-T(\mathbf{x})$，在实现过程中还有一步Variational refinement的过程，但是这一步不是必要的，在OpenCV的实现中这一步也是一个可选项，在此就不进行详细介绍，说道OpenCV，在OpenCV4.0中集成了DIS光流算法，如果大家想复现这个时域降噪算法，可以继承OpenCV4.0中的DIS稠密光流类来进一步实现。

以上将了该时域降噪算法中对齐部分使用的DIS稠密光流算法，接下来该算法是如何进行融合的，有两点是值得注意的：
（1）通过高斯金字塔求解DIS光流法，通过拉普拉斯金字塔进行融合
还不了解拉普拉斯金字塔的同学可以自行百度下，拉普拉斯金字塔可以无损地恢复原始图像，相对直接使用高斯金字塔进行融合显然能保留更多细节
（2）使用前一帧的降噪结果和当前帧进行融合，通过设计融合权重保证正确性
首先进行像素对齐$${\mathcal{L}}_a^l(\mathbf{p})={\mathcal{L}}_p^l\left(\mathbf{p}+A^l(\mathbf{p})\right)$$其中${\mathcal{L}}_a^l(\mathbf{p})$为第$l$层对齐的拉普拉斯金字塔，融合公式为：$${\mathcal{L}}_r^l(\mathbf{p})={\mathcal{I}}_c(\mathbf{p})\cdot {\mathcal{L}}_c^l(\mathbf{p})+{\mathcal{I}}_p(\mathbf{p})\cdot {\mathcal{L}}_a^l(\mathbf{p})$$我们定义${\mathcal{L}}_{\Delta }^l(\mathbf{p})={\mathcal{L}}_c^l(\mathbf{p})-{\mathcal{L}}_a^l(\mathbf{p})$，那么最终的融合公式变为：$${\mathcal{L}}_r^l(\mathbf{p})=w_c^l\cdot {\mathcal{L}}_c^l(\mathbf{p})+w_p^l\cdot \left({\mathcal{L}}_a^l(\mathbf{p})+\mathcal{I}\cdot {\mathcal{L}}_{\Delta }^l(\mathbf{p})\right)$$ 其中$w_c^l$和$w_p^l$两个系数用于控制降噪强度，$w_c^l$越大细节保留越好，$w_p^l$越大降噪强度越大，为了同时保留好的图像细节以及降噪强度要求$w_c^l+w_p^l\ge 1$。
系数$\mathcal{I}$存在的目的主要是为了消除由于图像像素没有完全对齐进行融合而造成的"鬼影"问题，系数$\mathcal{I}$的取值范围是$[0,1]$，系数$\mathcal{I}$越大则降噪强度越小，在本算法中，将系数$\mathcal{I}$定义为$\mathcal{I}=\max \left({\mathcal{I}}_e,{\mathcal{I}}_{\Delta }\right)$，其中${\mathcal{I}}_e$是和对齐误差有关，${\mathcal{I}}_{\Delta }$和像素误差有关，其中$${\mathcal{I}}_e=\max \left(1,A_e^l(\mathbf{p})\cdot C_e\right)$$其中$A_e^l(\mathbf{p})$是在对齐阶段的图像块匹配误差，$C_e$是可以设定的参数，当图像中出现遮挡、运动的物体时，${\mathcal{I}}_e$就会控制算法关闭这一区域的降噪。另外$${\mathcal{I}}_{\Delta }(\mathbf{p})={\left(1+{\exp }^{-\left(|{\mathcal{L}}_{\Delta }^l(\mathbf{p})|-m\right)}\right)}^{-1}$$$$m(\mathbf{p})=1+C_{\text{middle}}\cdot \left(1-{\exp }^{-n^l(\mathbf{p})\cdot C_{\text{noise}}^l}\right)$$其中$n^l(\mathbf{p})$为噪声强度，$C_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}^l$和$C_{\text{middle }}$为可以设定系数，可以发现当$\left|{\mathcal{L}}_{\Delta }^l(\mathbf{p})\right|=m(\mathbf{p})$是，系数${\mathcal{I}}_{\Delta }(\mathbf{p})$为0.5，${\mathcal{I}}_{\Delta }(\mathbf{p})$随噪声强度和像素误差的关系图如下：

可以看到，噪声强度小的时候，像素误差小降噪强度也越小（红色区域），噪声强度大的时候，像素误差小也会有较大的降噪强度（蓝色区域），这个算法从论文的效果看还是很牛逼的，但是自己复现效果怎么样就得看实现了。

3.《Analysis of Optical Flow Algorithms for Denosing》

这篇论文应该是2015年的一篇课程总结还是啥，里面有很长的篇幅都是在划水介绍一些基本知识，在后半部分主要是通过实验对比了四种不同的光流算法对降噪效果的影响，并且提出一个结论说先空域降噪后再时域降噪效果会更好，对时域降噪不了解并且时间充裕的同学可以看下这篇论文，这里就不展开了。

除了这几篇论文外，时域降噪其实是一个大的领域，包括深度学习例如LSTM在这个领域也都有实现，还有基于卡尔曼滤波实现的，但是这些偏学术性会更强，之后有时间再去探究，这篇博客就到这里，有问题欢迎交流～

此外，这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结，对图像降噪算法感兴趣的同学欢迎参考