圖像濾波（Image Denoising）

? ? ? ? 在討論所有問題之前，首先要對濾波這個問題進行良好的數(shù)學建模，一般來講光學圖像濾波都會采用如下的模型

y_j = z(\boldsymbol x_j) + e \quad j = 1,2,\cdots

其中\boldsymbol x_j代表了空間位置信息，y_j=y(\boldsymbol x_j)代表了對應(yīng)\boldsymbol x_j位置上的觀察到的像素值信息，z_j=z(\boldsymbol x_j)代表了對應(yīng)\boldsymbol x_j位置上真實的像素值信息，是需要求解的對象，常常被稱為ground truth，e是均值為0的隨機噪聲。

? ? ?? 如果對同一目標z_j進行多次觀測，為了從這一系列觀測到的值\boldsymbol y = [y_1, y_2, \cdots, y_N]^T中恢復(fù)出真實的信號，根據(jù)貝葉斯的理論框架，如果已知獨立同分布的噪聲e的概率密度函數(shù)p(e)，有表達式

p(z_j|\boldsymbol y) = \frac{p(\boldsymbol y|z_j)p(z_j)}{p(\boldsymbol y)}\propto p(\boldsymbol y|z_j)p(z_j)

濾波問題就是在求解如下的方程

\hat z_j = \arg\mathop{\max}_{z_j} p(z_j | \boldsymbol y)

在均勻先驗p(z_j) = \text{const}的假設(shè)下，如果e是零均值高斯白噪聲\mathcal{N}(0, \sigma^2 \boldsymbol I)，存在如下的表達式

\begin{aligned} \hat z_j =& \arg\mathop{\max}_{z_j} p(z_j | \boldsymbol y) = \arg\mathop{\max}_{z_j} p(\boldsymbol y | z_j)p(z_j)\\ =& \arg\mathop{\max}_{z_j} e^{-\frac{1}{2\sigma^2}(\boldsymbol y - z_j \boldsymbol 1)^T (\boldsymbol y - z_j \boldsymbol 1)} = \arg\mathop{\min}_{z_j} (\boldsymbol y - z_j \boldsymbol 1)^T (\boldsymbol y - z_j \boldsymbol 1) \\ =& \arg\mathop{\min}_{z_j} \sum_{i=1}^N (y_i - z_j)^2 \end{aligned}

很顯然，上式的解就是常見的均值濾波器。它建立在對于同一個目標的多次觀測的假設(shè)基礎(chǔ)上。在圖像處理中，一般認為空間上臨近的像素其相同的可能性大，因此往往在目標點的空間領(lǐng)域內(nèi)進行均值操作。但是上式的要求是嚴格的，它要求周圍像素點和中心像素點擁有相同的真值。可以想象，這種方法在勻質(zhì)區(qū)域上擁有較好的效果，事實上它是勻質(zhì)區(qū)域在最小方差無偏意義上最好的濾波器，達到了該問題的克拉美羅界。對于非勻質(zhì)區(qū)域，周圍像素點不再和中心像素點保持絕對一致，該濾波器無偏差地濾掉了屬于信號的高頻信息（比如邊緣信息），引起了邊緣模糊等現(xiàn)象。因此后來引入了擁有權(quán)重的濾波器，對這些周圍領(lǐng)域像素點賦予不同的權(quán)重，和中心像素點越相似的權(quán)重越大，越不相似的權(quán)重越小。這些濾波器可以統(tǒng)一地用下面的式子來表達

\hat z_j = \arg\mathop{\min}_{z_j} \sum_{i=1}^N (y_i - z_j)^2 K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j)

其中的權(quán)重K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j)是一個對稱的正值的函數(shù)，它代表了像素之間的相似性和距離。很多濾波算法都可以用上式來表達，它們之間的區(qū)別就在于這個權(quán)重函數(shù)K的不同。前面的均值濾波器也是其中的一種特殊情況，其中K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j) = \frac{1}{N} = \text{const}。此外還有常見的高斯濾波器（Gaussian Filter，GF），它把權(quán)重函數(shù)寫成只和位置\boldsymbol x相關(guān)的函數(shù)，認為在空間上距離\boldsymbol x_j越近相似度越大，距離越遠相似度越小，其中相似度權(quán)重和距離的關(guān)系采用高斯函數(shù)來描述，如果距離度量采用歐氏距離，其權(quán)重可以用下式表達

K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j) = e^{-\frac{||\boldsymbol x_i - \boldsymbol x_j||^2}{\sigma_x^2}}

除了高斯函數(shù)之外，還有很多類似的只與空間位置相關(guān)的核函數(shù)。遺憾的是這樣的描述方式只考慮了空間上的關(guān)系，缺少了基于內(nèi)容的修正，因此并沒有解決邊緣、角點這樣的高頻信號的模糊問題，因此后來提出了雙邊濾波器，其權(quán)重函數(shù)如下式所示

K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j) = e^{-\frac{||\boldsymbol x_i - \boldsymbol x_j||^2}{\sigma_x^2}-\frac{||y_i - y_j||^2}{\sigma_y^2}}

這種雙邊濾波器考慮到了像素值之間的相似程度，對于含有邊緣的圖像能夠起到一定的保邊作用。但是雙邊濾波器對于低信噪比下的圖像，性能會惡化。這是因為噪聲的干擾，處于同一邊的像素點之間的相似程度也會很低。其根本原因在于位置\boldsymbol x的值是精確的，但是像素y的值是近似的，直接采用單個像素的y值會受到噪聲的影響，為了增加濾波器對于噪聲的抵抗力，需要用周圍領(lǐng)域來估計和y相關(guān)的參數(shù)，因此后面又有學者提出了局部自適應(yīng)回歸核LARK，其權(quán)重核函數(shù)如下所示

\begin{aligned} K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j) =& e^{-(\boldsymbol x_i - \boldsymbol x_j)^T \boldsymbol C_i (\boldsymbol x_i - \boldsymbol x_j)} \\ \boldsymbol C_i =& \sum_{y_k \in window} \begin{bmatrix} \left(\frac{\partial y_k}{\partial x_1}\right)^2 & \frac{\partial y_k}{\partial x_1}\frac{\partial y_k}{\partial x_2} \\ \frac{\partial y_k}{\partial x_1}\frac{\partial y_k}{\partial x_2} & \left(\frac{\partial y_k}{\partial x_2}\right)^2 \\ \end{bmatrix} \end{aligned}

這種方法首先估計窗內(nèi)每個點的梯度協(xié)方差矩陣，然后求平均，這一步的本質(zhì)就是在求解窗內(nèi)邊緣的方向，以便于制定和圖像結(jié)構(gòu)相關(guān)的濾波器核函數(shù)。相比于前面的均勻濾波和高斯濾波，LARK權(quán)重系數(shù)在每個方向上下降的速度是不一樣的，和梯度的大小成正相關(guān)，因此它考慮了圖像本身的結(jié)構(gòu)信息。相比于雙邊濾波器，LARK利用周圍領(lǐng)域來估計和y相關(guān)的系數(shù)，抗噪聲能力更強。LARK可以認為是一種從數(shù)據(jù)中學習魯棒的結(jié)構(gòu)探針的方法。

為了進一步提升性能，需要將更多的相似像素加入到這個估計式子中來，因此有必要將周圍領(lǐng)域擴大到更大的范圍，比如說全圖，這就是非鄰域平均算法（Non-Local Mean,NLM）的出發(fā)點?？梢杂孟旅娴氖阶觼肀磉_

K(\boldsymbol x_i, \boldsymbol x_j, y_i, y_j) = \lim\limits_{\sigma_x \to \infty} e^{-\frac{||\boldsymbol x_i - \boldsymbol x_j||^2}{\sigma_x^2}-\frac{||y_i - y_j||^2}{\sigma_y^2}} = e^{-\frac{||y_i - y_j||^2}{\sigma_y^2}}

它完全忽視了距離對于權(quán)重的影響，只考慮圖像像素和結(jié)構(gòu)之間的相似性，這樣可以把更多的屬于同一目標的像素點加入到計算中來，進一步降低估計結(jié)果的方差，提高泛化性能。遺憾的是如果嚴格按照公式來，對每一點都要做這樣的全圖相似性計算操作，其計算復(fù)雜度急劇上升，因此在應(yīng)用中往往還是采用local的鄰域方案，只不過是一個相對更大的鄰域而已。如果Non-Local算法只根據(jù)單個像素的觀測值y來計算權(quán)重，同樣對于低信噪比下的結(jié)果會造成性能惡化。一般采用基于小塊的Non-Local算法來增強算法的魯棒性。但是對于那些不具有自相似性的圖像，會過度平滑復(fù)雜不規(guī)則結(jié)構(gòu)。

如果對式用矩陣的方式來描述，那么可以寫成

\begin{aligned} \hat z_j =& \arg\mathop{\min}_{z_j} (\boldsymbol y - z_j \boldsymbol 1)^T \begin{bmatrix} K(\boldsymbol x_1, \boldsymbol x_j, y_1, y_j) & & & \\ & K(\boldsymbol x_2, \boldsymbol x_j, y_2, y_j) & & \\ & & \ddots & \\ & & & K(\boldsymbol x_N, \boldsymbol x_j, y_N, y_j) \end{bmatrix} (\boldsymbol y - z_j \boldsymbol 1) \\ =& \arg\mathop{\min}_{z_j} (\boldsymbol y - z_j \boldsymbol 1)^T \boldsymbol K_j (\boldsymbol y - z_j \boldsymbol 1) \end{aligned}

上面方程的閉式解是

\hat z_j = (\boldsymbol 1^T \boldsymbol K_j \boldsymbol 1)^{-1} \boldsymbol 1^T \boldsymbol K_j \boldsymbol y = \boldsymbol w_j^T \boldsymbol y

上式是對單點濾波的閉式解，如果是對整個區(qū)域里的點濾波則可以寫成

\boldsymbol{\hat z} = \begin{bmatrix} \boldsymbol w_1^T \\ \boldsymbol w_2^T \\ \vdots \\ \boldsymbol w_n^T \end{bmatrix} \boldsymbol y = \boldsymbol W^T \boldsymbol y

可以看到前面這些空域方法的核心就是權(quán)重參數(shù)\boldsymbol W，雖然不同方法得到的參數(shù)不一樣，但它們都是周圍鄰域點的函數(shù)。如果對式稍作變換可以得到

\boldsymbol{\hat z} = \boldsymbol D (\boldsymbol D^T\boldsymbol W^T \boldsymbol D) \boldsymbol D^T \boldsymbol y = \boldsymbol D \boldsymbol W

這就是變換域濾波算法，先把信號\boldsymbol y通過正交矩陣\boldsymbol D^T變換到變換域上，然后在變換域通過矩陣\boldsymbol W進行濾波，再通過正交矩陣\boldsymbol D變換到原域上。常用的變換有傅里葉變換、小波變換、DCT變換等等。

前面提到的算法都是從式p(z_j|\boldsymbol y)\propto p(\boldsymbol y|z_j)p(z_j)出發(fā)的，它們有一個共同的特性，就是都在均勻先驗p(z_j) = \text{const}的假設(shè)下。這種假設(shè)認為濾波的圖像的像素值均勻分布在整個色彩空間，而事實上，實際中的圖像在整個色彩空間是有一定偏好的。如果利用這種先驗信息，那么可以得到更好的結(jié)果。

\begin{aligned} \boldsymbol{\hat z} =& \arg\mathop{\max}_{\boldsymbol z} p(\boldsymbol y | \boldsymbol z) p(\boldsymbol z) \\ =& \arg\mathop{\max}_{\boldsymbol z} e^{-\frac{1}{2\sigma^2}||\boldsymbol y - \boldsymbol z||^2} e^{\ln p(\boldsymbol z)} \\ =& \arg\mathop{\min}_{\boldsymbol z} ||\boldsymbol y - \boldsymbol z||^2 + \lambda \Phi(\boldsymbol z) \end{aligned}

這種先驗包括非常流行的稀疏先驗假設(shè),稀疏先驗假設(shè)認為任意一個圖像\boldsymbol z都可以用一個完備的字典\boldsymbol D稀疏表達，它的先驗就是字典表達的稀疏性，也被稱為稀疏編碼（Sparse Coding，SC），如下式所示

\begin{aligned} \boldsymbol{\hat x} =& \arg\mathop{\min}_{\boldsymbol x} ||\boldsymbol y - \boldsymbol D \boldsymbol x||^2 + \lambda ||\boldsymbol x||_0 \\ \boldsymbol{\hat z} =& \boldsymbol D \boldsymbol{\hat x} \end{aligned}

其中\boldsymbol x是圖像\boldsymbol z在字典\boldsymbol D下的表達。

還有一種很流行的先驗是變分（Total Variation，TV）約束先驗。這種先驗的想法來源于Rudin在文獻中提到的一個經(jīng)驗性的觀察：含有噪聲的圖像比沒有噪聲的圖片，其變分絕對值之和要大很多。因此這可以作為一種信息對優(yōu)化加以約束。具體表達式如下

\boldsymbol{\hat z} = \arg\mathop{\min}_{\boldsymbol z} ||\boldsymbol y - \boldsymbol z||^2 + \lambda ||\triangledown \boldsymbol z||_1

在這些約束的作用下，使得濾波器更加傾向于修復(fù)符合先驗的圖像（概率大的\boldsymbol z），而不再是在全空間修復(fù)，對于特定問題可以超越一般的算法，先驗越符合實際，效果就越好。

注意到前面的這些方法，往往需要真值的一些特征參數(shù)，包括像素值、邊緣主方向等等。這些參數(shù)由于噪聲的干擾，本身也變成了一個隨機變量。前面已經(jīng)出現(xiàn)了用一個小patch來估計這些參數(shù)的方法，大大改善了估計的準確度。除此以外，還出現(xiàn)了循環(huán)方法，這些循環(huán)包括串聯(lián)濾波器，使得下一次的參數(shù)基于上一次濾波之后的更準確的結(jié)果，以達到比上一次更加準確的最終估計值（偏差不斷變大，方差不斷變?。?，也包括殘差濾波，不斷挖掘被濾掉部分中的高頻信息疊加到最后的結(jié)果中（偏差不斷變小，方差不斷變大）。它們不斷疊加的結(jié)果就是方差和偏差之間的權(quán)衡，在均方誤差（Mean Squared Error，MSE）最小的情況下，線性濾波器通過最優(yōu)迭代最終演變成維納濾波，達到線性均方最小誤差問題（Linear Minimun Mean Square Error，LMMSE）的克拉美羅界。

文獻Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering提出了一種新的濾波方法BM3D，它參考了NLM在很大的領(lǐng)域內(nèi)搜索相似塊，并且在變換域進行濾波，同時考慮了空間的位置關(guān)系和塊之間的特征相似性，并且采用了迭代的方法。使得下一次得到的維納濾波器參數(shù)估計更加準確，達到了經(jīng)典濾波算法的state of art。該方法非常復(fù)雜，以至于沒法簡單的用數(shù)學表達式來表達權(quán)重參數(shù)\boldsymbol W。具體的細節(jié)可以參考文獻。

這些傳統(tǒng)方法都有一個很大的缺陷，就是需要事先得知圖像的噪聲強度或者信噪比來估算其中的某些參數(shù)，這給實際應(yīng)用帶來了很大的限制。有時候甚至在同一張圖像上的噪聲強度并不是均勻的，而是有差別的，而這些方法都不能應(yīng)對這樣的情況。更重要的是，它們本身就是淺層的神經(jīng)網(wǎng)絡(luò)，表達能力有限。因此后面出現(xiàn)了采用神經(jīng)網(wǎng)絡(luò)的方法。

未完待續(xù)……

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容