圖像降噪算法——時域降噪算法

• 圖像降噪算法——時域降噪算法
• • 1.《MeshFLow Video Denoising》
  • 2. 《Real-Time Video Denoising On Mobile Phones》
  • 3.《Analysis of Optical Flow Algorithms for Denosing》

圖像降噪算法——時域降噪算法

最近在工作上接觸到了時域降噪相關的算法，這里進行一個簡單的總結。我們可以這樣理解降噪算法，降噪的過程其實就是通過觀測減小噪聲方差的過程，對于單幀降噪，做法可以是通過在圖片上尋找類似的像素，通過加權平均來獲得減少噪聲的方差，最經典的就是非局部均值濾波算法，而對于時域降噪，這個概念就更加直接，例如我在第一幀我看到這個像素值是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))+?T\frac{?\mathbf{W}}{?\mathbf{u}}\Delta \mathbf{u}?I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))\right]}^2$$其中$?T$為圖像塊$T$梯度圖像，$?\mathbf{W}/?\mathbf{u}$為$\mathbf{u}=0$時的雅克比。為對上式求導，并令導數等于0，得到下式：$$2\sum _x{\left[?T\frac{?\mathbf{W}}{?\mathbf{u}}\right]}^T\left[T(\mathbf{W}(\mathbf{x};0))+?T\frac{?\mathbf{W}}{?\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[?T\frac{?\mathbf{W}}{?\mathbf{u}}\right]}^T[T(\mathbf{W}(\mathbf{x};0))?I_{t+1}(\mathbf{W}(\mathbf{x};\mathbf{u}))]$$其中$$H=\sum _x{\left[?T\frac{?\mathbf{W}}{?\mathbf{u}}\right]}^T\left[?T\frac{?\mathbf{W}}{?\mathbf{u}}\right]$$由于$?\mathbf{W}/?\mathbf{u}$為$\mathbf{u}=0$時的雅克比，$?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)}?{\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})?{\mathcal{L}}_c^l(\mathbf{p})+{\mathcal{I}}_p(\mathbf{p})?{\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?{\mathcal{L}}_c^l(\mathbf{p})+w_p^l?\left({\mathcal{L}}_a^l(\mathbf{p})+\mathcal{I}?{\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})?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}}?\left(1?{\exp }^{?n^l(\mathbf{p})?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在這個領域也都有實現，還有基于卡爾曼濾波實現的，但是這些偏學術性會更強，之后有時間再去探究，這篇博客就到這里，有問題歡迎交流～

此外，這里我寫一個各種算法的總結目錄圖像降噪算法——圖像降噪算法總結，對圖像降噪算法感興趣的同學歡迎參考

總結

以上是生活随笔為你收集整理的图像降噪算法——时域降噪算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。