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

• 圖像降噪算法——時域降噪算法
• • 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 I_{t} It​中的一個圖像塊 T T T,設其尺寸為 θ p s × θ p s \theta_{p s} \times \theta_{p s} θps​×θps​，坐标位置為 x = ( x , y ) T \mathbf{x}=(x, y)^{T} x=(x,y)T，我們需要使用梯度下降法在幀 I t + 1 I_{t+1} It+1​尋找到最佳比對的尺寸同樣為 θ p s × θ p s \theta_{p s} \times \theta_{p s} θps​×θps​的圖像塊，實際上就是尋找一個運動向量 u = ( u , v ) \mathbf{u}=(u, v) u=(u,v)使得前後幀的圖像塊平方差最小： u = argmin ⁡ u ′ ∑ x [ I t + 1 ( x + u ′ ) − T ( x ) ] 2 \mathbf{u}=\operatorname{argmin}_{\mathbf{u}^{\prime}} \sum_{x}\left[I_{t+1}\left(\mathbf{x}+\mathbf{u}^{\prime}\right)-T(\mathbf{x})\right]^{2} u=argminu′​x∑​[It+1​(x+u′)−T(x)]2例如LK光流就是将其作為一個非線性優化問題，運動向量 u = ( u , v ) \mathbf{u}=(u, v) u=(u,v)作為優化變量，使用疊代法進行求解。在DIS稠密光流中對該方法進行了修改。DIS稠密光流将優化的目标修改為目前運動估計 u \mathbf{u} u基礎上的一個更新向量 Δ u \Delta \mathbf{u} Δu： Δ u = argmin ⁡ Δ u ′ ∑ x [ I t + 1 ( x + u + Δ u ′ ) − T ( x ) ] 2 \Delta \mathbf{u}=\operatorname{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} Δu=argminΔu′​x∑​[It+1​(x+u+Δu′)−T(x)]2進一步地将上式的目标損失函數修改為 ∑ x [ T ( x − Δ u ) − I t + 1 ( x + u ) ] 2 \sum_{x}\left[T(\mathbf{x}-\Delta \mathbf{u})-I_{t+1}(\mathbf{x}+\mathbf{u})\right]^{2} ∑x​[T(x−Δu)−It+1​(x+u)]2再進行求解，修改成這樣有什麼好處呢？這樣的修改友善我們實作LK光流的逆向組成算法，為了進一步了解逆向組成算法的流程，我們先對目标損失函數進行進一步分解：

我們先将上式改寫為： ∑ x [ T ( W ( x ; Δ u ) ) − I t + 1 ( W ( x ; u ) ) ] 2 \sum_{x}\left[T(\mathbf{W}(\mathbf{x} ; \Delta \mathbf{u}))-I_{t+1}(\mathbf{W}(\mathbf{x} ; \mathbf{u}))\right]^{2} x∑​[T(W(x;Δu))−It+1​(W(x;u))]2其中 W ( x ; u ) \mathbf{W}(\mathbf{x} ; \mathbf{\mathbf{u}}) W(x;u)描述的是一種映射關系 W ( x ; u ) = ( x − u y − v ) \mathbf{W}(\mathbf{x} ; \mathbf{\mathbf{u}})=\left(\begin{array}{l} x-u \\ y-v \end{array}\right) W(x;u)=(x−uy−v​)我們對上式進行一階泰勒展開有： ∑ x [ T ( W ( x ; 0 ) ) + ∇ T ∂ W ∂ u Δ u − I t + 1 ( W ( x ; u ) ) ] 2 \sum_{x}\left[T(\mathbf{W}(\mathbf{x} ; \mathbf{0}))+\nabla T \frac{\partial \mathbf{W}}{\partial \mathbf{u}} \Delta \mathbf{u}-I_{t+1}(\mathbf{W}(\mathbf{x} ; \mathbf{u}))\right]^{2} x∑​[T(W(x;0))+∇T∂u∂W​Δu−It+1​(W(x;u))]2其中 ∇ T \nabla T ∇T為圖像塊 T T T梯度圖像， ∂ W / ∂ u {\partial \mathbf{W}}/{\partial \mathbf{u}} ∂W/∂u為 u = 0 \mathbf{u}=\mathbf{0} u=0時的雅克比。為對上式求導，并令導數等于0，得到下式： 2 ∑ x [ ∇ T ∂ W ∂ u ] T [ T ( W ( x ; 0 ) ) + ∇ T ∂ W ∂ u Δ u − I t + 1 ( W ( x ; u ) ) ] = 0 2 \sum_{x}\left[\nabla T \frac{\partial \mathbf{W}}{\partial \mathbf{u}}\right]^{T}\left[T(\mathbf{W}(\mathbf{x}; \mathbf{0}))+\nabla T \frac{\partial \mathbf{W}}{\partial \mathbf{u}} \Delta \mathbf{u}-I_{t+1}(\mathbf{W}(\mathbf{x} ; \mathbf{u}))\right]=0 2x∑​[∇T∂u∂W​]T[T(W(x;0))+∇T∂u∂W​Δu−It+1​(W(x;u))]=0上式對 Δ u \Delta \mathbf{u} Δu進行求解，得到 Δ u \Delta \mathbf{u} Δu的解析解為： Δ u = H − 1 ∑ x [ ∇ T ∂ W ∂ u ] T [ T ( W ( x ; 0 ) ) − I t + 1 ( W ( x ; 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}; \mathbf{0}))-I_{t+1}(\mathbf{W}(\mathbf{x}; \mathbf{u}))] Δu=H−1x∑​[∇T∂u∂W​]T[T(W(x;0))−It+1​(W(x;u))]其中 H = ∑ x [ ∇ T ∂ W ∂ u ] T [ ∇ T ∂ W ∂ 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] H=x∑​[∇T∂u∂W​]T[∇T∂u∂W​]由于 ∂ W / ∂ u {\partial \mathbf{W}}/{\partial \mathbf{u}} ∂W/∂u為 u = 0 \mathbf{u}=\mathbf{0} u=0時的雅克比， ∇ T \nabla T ∇T為梯度圖像塊 T T T的梯度圖像，這兩者都是定值，是以矩陣 H H H也是定值，這樣就免去了我們在疊代優化 Δ u \Delta \mathbf{u} Δu時反複求解 H H H矩陣的問題。對上面這段推導有點懵逼的同學可以參考下論文《Lucas-Kanade 20 Years On: A Unifying Framework》，這篇論文将光流法總結為四種方法，按照更新量的不同分為additional和compositional兩種，按照更新方式不同分為forward和inverse兩種，上面推導的這種inverse ccompositional(逆向組成)算法是其中一種。

在完成圖像塊的光流追蹤之後，目前幀中的像素可能被參考幀中的多個圖像塊覆寫，我們對所有覆寫該像素的圖像塊光流進行權重平均就獲得該點的光流值： U s ( x ) = 1 Z ∑ i N s λ i , x max ⁡ ( 1 , ∥ d i ( x ) ∥ 2 ) ⋅ u i \mathbf{U}_{s}(\mathbf{x})=\frac{1}{Z} \sum_{i}^{N_{s}} \frac{\lambda_{i, \mathbf{x}}}{\max \left(1,\left\|d_{i}(\mathbf{x})\right\|_{2}\right)} \cdot \mathbf{u}_{i} Us​(x)=Z1​i∑Ns​​max(1,∥di​(x)∥2​)λi,x​​⋅ui​其中 u i \mathbf{u}_{i} ui​為圖像塊的光流， U s ( x ) \mathbf{U}_{s}(\mathbf{x}) Us​(x)為像素 x \mathbf{x} x處的光流大小， d i ( x ) = I t + 1 ( x + u i ) − T ( x ) d_{i}(\mathbf{x})=I_{t+1}\left(\mathbf{x}+\mathbf{u}_{i}\right)-T(\mathbf{x}) di​(x)=It+1​(x+ui​)−T(x)，在實作過程中還有一步Variational refinement的過程，但是這一步不是必要的，在OpenCV的實作中這一步也是一個可選項，在此就不進行詳細介紹，說道OpenCV，在OpenCV4.0中內建了DIS光流算法，如果大家想複現這個時域降噪算法，可以繼承OpenCV4.0中的DIS稠密光流類來進一步實作。

以上将了該時域降噪算法中對齊部分使用的DIS稠密光流算法，接下來該算法是如何進行融合的，有兩點是值得注意的：

（1）通過高斯金字塔求解DIS光流法，通過拉普拉斯金字塔進行融合

還不了解拉普拉斯金字塔的同學可以自行百度下，拉普拉斯金字塔可以無損地恢複原始圖像，相對直接使用高斯金字塔進行融合顯然能保留更多細節

（2）使用前一幀的降噪結果和目前幀進行融合，通過設計融合權重保證正确性

首先進行像素對齊 L a l ( p ) = L p l ( p + A l ( p ) ) \mathcal{L}_{a}^{l}(\mathbf{p})=\mathcal{L}_{p}^{l}\left(\mathbf{p}+A^{l}(\mathbf{p})\right) Lal​(p)=Lpl​(p+Al(p))其中 L a l ( p ) \mathcal{L}_{a}^{l}(\mathbf{p}) Lal​(p)為第 l l l層對齊的拉普拉斯金字塔，融合公式為： L r l ( p ) = I c ( p ) ⋅ L c l ( p ) + I p ( p ) ⋅ L a l ( p ) \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}) Lrl​(p)=Ic​(p)⋅Lcl​(p)+Ip​(p)⋅Lal​(p)我們定義 L Δ l ( p ) = L c l ( p ) − L a l ( p ) \mathcal{L}_{\Delta}^{l}(\mathbf{p})=\mathcal{L}_{c}^{l}(\mathbf{p})-\mathcal{L}_{a}^{l}(\mathbf{p}) LΔl​(p)=Lcl​(p)−Lal​(p)，那麼最終的融合公式變為： L r l ( p ) = w c l ⋅ L c l ( p ) + w p l ⋅ ( L a l ( p ) + I ⋅ L Δ l ( 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) Lrl​(p)=wcl​⋅Lcl​(p)+wpl​⋅(Lal​(p)+I⋅LΔl​(p)) 其中 w c l w_{c}^{l} wcl​和 w p l w_{p}^{l} wpl​兩個系數用于控制降噪強度， w c l w_{c}^{l} wcl​越大細節保留越好， w p l w_{p}^{l} wpl​越大降噪強度越大，為了同時保留好的圖像細節以及降噪強度要求 w c l + w p l ≥ 1 w_{c}^{l}+w_{p}^{l} \geq 1 wcl​+wpl​≥1。

系數 I \mathcal{I} I存在的目的主要是為了消除由于圖像像素沒有完全對齊進行融合而造成的"鬼影"問題，系數 I \mathcal{I} I的取值範圍是 [ 0 , 1 ] [0,1] [0,1]，系數 I \mathcal{I} I越大則降噪強度越小，在本算法中，将系數 I \mathcal{I} I定義為 I = max ⁡ ( I e , I Δ ) \mathcal{I}=\max \left(\mathcal{I}_{e}, \mathcal{I}_{\Delta}\right) I=max(Ie​,IΔ​)，其中 I e \mathcal{I}_{e} Ie​是和對齊誤差有關， I Δ \mathcal{I}_{\Delta} IΔ​和像素誤差有關，其中 I e = max ⁡ ( 1 , A e l ( p ) ⋅ C e ) \mathcal{I}_{e}=\max \left(1, A_{e}^{l}(\mathbf{p}) \cdot C_{e}\right) Ie​=max(1,Ael​(p)⋅Ce​)其中 A e l ( p ) A_{e}^{l}(\mathbf{p}) Ael​(p)是在對齊階段的圖像塊比對誤差， C e C_{e} Ce​是可以設定的參數，當圖像中出現遮擋、運動的物體時， I e \mathcal{I}_{e} Ie​就會控制算法關閉這一區域的降噪。另外 I Δ ( p ) = ( 1 + exp ⁡ − ( ∣ L Δ l ( p ) ∣ − m ) ) − 1 \mathcal{I}_{\Delta}(\mathbf{p})=\left(1+\exp ^{-\left(\left|\mathcal{L}_{\Delta}^{l}(\mathbf{p})\right|-m\right)}\right)^{-1} IΔ​(p)=(1+exp−(∣LΔl​(p)∣−m))−1 m ( p ) = 1 + C middle ⋅ ( 1 − exp ⁡ − n l ( p ) ⋅ C noise l ) m(\mathbf{p})=1+C_{\text {middle}} \cdot\left(1-\exp ^{-n^{l}(\mathbf{p}) \cdot C_{\text {noise}}^{l}}\right) m(p)=1+Cmiddle​⋅(1−exp−nl(p)⋅Cnoisel​)其中 n l ( p ) n^{l}(\mathbf{p}) nl(p)為噪聲強度， C n o i s e l C_{\mathrm{noise}}^{l} Cnoisel​和 C middle  C_{\text {middle }} Cmiddle ​為可以設定系數，可以發現當 ∣ L Δ l ( p ) ∣ = m ( p ) \left|\mathcal{L}_{\Delta}^{l}(\mathbf{p})\right|=m(\mathbf{p}) ∣∣​LΔl​(p)∣∣​=m(p)是，系數 I Δ ( p ) \mathcal{I}_{\Delta}(\mathbf{p}) IΔ​(p)為0.5， I Δ ( p ) \mathcal{I}_{\Delta}(\mathbf{p}) IΔ​(p)随噪聲強度和像素誤差的關系圖如下：

可以看到，噪聲強度小的時候，像素誤差小降噪強度也越小（紅色區域），噪聲強度大的時候，像素誤差小也會有較大的降噪強度（藍色區域），這個算法從論文的效果看還是很牛逼的，但是自己複現效果怎麼樣就得看實作了。

3.《Analysis of Optical Flow Algorithms for Denosing》

這篇論文應該是2015年的一篇課程總結還是啥，裡面有很長的篇幅都是在劃水介紹一些基本知識，在後半部分主要是通過實驗對比了四種不同的光流算法對降噪效果的影響，并且提出一個結論說先空域降噪後再時域降噪效果會更好，對時域降噪不了解并且時間充裕的同學可以看下這篇論文，這裡就不展開了。

除了這幾篇論文外，時域降噪其實是一個大的領域，包括深度學習例如LSTM在這個領域也都有實作，還有基于卡爾曼濾波實作的，但是這些偏學術性會更強，之後有時間再去探究，這篇部落格就到這裡，有問題歡迎交流～

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