導(dǎo)向濾波(Guided Filter)的解析與實現(xiàn)

現(xiàn)在從一個最簡單的情形來開始我們的討論。假設(shè)有一個原始圖像 pp，其中含有一些噪聲，欲將這些噪聲濾出，最簡單的、最基本的方法，大家可能會想到采用一些低通濾波器，例如簡單平滑（也稱Box Filter）或者高斯平滑等。濾波之后的圖像為qq，如下圖所示，圖像qq中第ii個像素是由圖像pp中以第ii個像素為中心的一個窗口ww中的像素確定的。

具體而言，在簡單平滑中，圖像qq中第ii個像素是由圖像pp中以第ii個像素為中心的一個窗口ww中的所有像素取平均而得來的，即

qi=∑j∈wiWij?pj

其中，Wij=1/nWij=1/n，nn是窗口ww中的像素數(shù)目。也就是說，在簡單平滑中，以像素ii為中心的一個窗口ww中的像素具有等同的權(quán)值。但是在高斯平滑中，權(quán)值WijWij將服從二維的高斯分布，結(jié)果導(dǎo)致離像素ii更接近的像素將具有更高的權(quán)重，反之離像素ii較遠的像素則具有更小的權(quán)重。

無論是簡單平滑，還是高斯平滑，它們都有一個共同的弱點，即它們都屬于各向同性濾波。我們都知道，一幅自然的圖像可以被看成是有（過渡平緩的，也就是梯度較?。﹨^(qū)域和（過渡尖銳的，也就是梯度較大）邊緣（也包括圖像的紋理、細節(jié)等）共同組成的。噪聲是影響圖像質(zhì)量的不利因素，我們希望將其濾除。噪聲的特點通常是以其為中心的各個方向上梯度都較大而且相差不多。邊緣則不同，邊緣相比于區(qū)域也會出現(xiàn)梯度的越變，但是邊緣只有在其法向方向上才會出現(xiàn)較大的梯度，而在切向方向上梯度較小。

因此，對于各向同性濾波（例如簡單平滑或高斯平滑）而言，它們對待噪聲和邊緣信息都采取一直的態(tài)度。結(jié)果，噪聲被磨平的同時，圖像中具有重要地位的邊緣、紋理和細節(jié)也同時被抹平了。這是我們所不希望看到的。研究人員已經(jīng)提出了很多Edge-perserving的圖像降噪（平滑）算法，例如雙邊濾波、自適應(yīng)（維納）平滑濾波（請參見文獻【1】）、基于PM方程的各向異性濾波以及基于TV-norm的降噪算法等。本文將考慮在文獻【2】中提出的另外一種Edge-perserving的圖像濾波（平滑）算法——導(dǎo)向濾波（Guided Filter）。當(dāng)然，通過閱讀文獻【2】，我們也知道導(dǎo)向濾波的應(yīng)用不止有Edge-perserving的圖像平滑，還包括圖像去霧、圖像Matting等等。

歡迎關(guān)注白馬負(fù)金羈的博客 http://blog.csdn.net/baimafujinji，為保證公式、圖表得以正確顯示，強烈建議你從該地址上查看原版博文。本博客主要關(guān)注方向包括：數(shù)字圖像處理、算法設(shè)計與分析、數(shù)據(jù)結(jié)構(gòu)、機器學(xué)習(xí)、數(shù)據(jù)挖掘、統(tǒng)計分析方法、自然語言處理。

算法原理

導(dǎo)向濾波之所以叫這個名字，因為在算法框架中，要對pp進行濾波而得到qq，還得需要一個引導(dǎo)圖像II。此時，濾波器的數(shù)學(xué)公式為

qi=∑j∈wiWij(I)?pj

注意，這里的Wij(I)Wij(I)就表示由引導(dǎo)圖像II來確定加權(quán)平均運算中所采用的權(quán)值。注意，引導(dǎo)圖像可以是單獨的一幅圖像，也可以輸入的圖像pp本身。當(dāng)引導(dǎo)圖像就是pp本身的時候，導(dǎo)向濾波就變成了一個Edge-perserving的濾波器，因此可以用于圖像的平滑降噪。

導(dǎo)向濾波濾波的示意圖如下所示（圖片來自作者原文【2】）。注意這也是導(dǎo)向濾波所依賴的一個重要假設(shè)——導(dǎo)向濾波器在導(dǎo)向圖像II和濾波輸出qq之間在一個二維窗口內(nèi)是一個局部線性模型，aa和bb是當(dāng)窗口中心位于kk時該線性函數(shù)的系數(shù)，即qi=akIi+bi,?i∈wkqi=akIi+bi,?i∈wk。導(dǎo)引圖像與qq之間存在線性關(guān)系，這樣設(shè)定是因為我們希望導(dǎo)引圖像提供的是信息主要用于指示哪些是邊緣哪些是區(qū)域，所以在濾波時，如果導(dǎo)引圖告訴我們這里是區(qū)域，那么就將其磨平。如果導(dǎo)引圖告訴我們這里是邊緣，這在最終的濾波結(jié)果里就要設(shè)法保留這些邊緣信息。只有當(dāng)II和qq之間是線性關(guān)系的這種引導(dǎo)關(guān)系才有意義。

你可以想象一種情況來理解這種假設(shè)。II中的甲處和乙處都是區(qū)域，而丙處是邊緣。那甲處的梯度和乙處的梯度應(yīng)該不會相差太大，例如乙處的梯度是甲處梯度的1.5倍。但是由于丙處是邊緣，所以丙處的梯度會比較大，例如可能是甲處的6倍，也就是乙處的4倍。如果II和qq之間滿足線性關(guān)系，那么甲乙丙的梯度大小和倍數(shù)關(guān)系就都不會被扭曲。否則，如果二者之間的關(guān)系是非線性的，那么可能的結(jié)果是在qq中，盡管丙處的梯度仍然大于乙處的梯度，進而大于甲處的梯度，但是倍數(shù)關(guān)系可能會扭曲。例如，丙處的梯度是乙處的1.5倍，而是甲處的6倍（即乙處的梯度是甲處的4倍），這時你就會想象，甲處是區(qū)域，而乙處和丙處就變成了邊緣?？梢姺蔷€性關(guān)系會使得引導(dǎo)圖像對于邊緣和區(qū)域的指示作用發(fā)生錯亂。

現(xiàn)在已知的是II和pp，要求的是qq。而如果能求得參數(shù)aa和bb，顯然就能通過II和qq之間的線性關(guān)系來求出II。另一方面，還可以知道pp是qq受到噪聲污染而產(chǎn)生的退化圖像，假設(shè)噪聲是nn，則有qi=pi?niqi=pi?ni。根據(jù)無約束圖像復(fù)原的方法（可以參考【4】），這時可以設(shè)定最優(yōu)化目標(biāo)為min||n||min||n||，等價地有minn2minn2，即min∑i∈wk(qi?pi)2min∑i∈wk(qi?pi)2，于是有

argmin∑i∈wk(akIi+bk?pi)2

argmin∑i∈wk(akIi+bk?pi)2

于是便得到了一個最小二乘問題。但是我們也知道普通最小二乘有時候引起一些麻煩，因此要適當(dāng)?shù)匾霊土P項，即采用正則化的手段。原作者在處理這里時所采用的方法借鑒了從普通線性回歸改進到嶺回歸時所采用的方法（如果你對正則化、嶺回歸、或者最小二乘法還不夠清楚，那么你可以參考文獻【5】），即求解下面這個最優(yōu)化（最小化）目標(biāo)所對應(yīng)的參數(shù)aa和bb。

E(ak,bk)=∑i∈wk[(akIi+bk?pi)2+?a2k]

求解上述最優(yōu)化問題（跟最小二乘法的推導(dǎo)過程一致），便會得到：

ak=1|w|∑i∈wkIipi?μipkˉσ2k+?bk=pkˉ?akμk

其中，μkμk是II中窗口wkwk中的平均值，σ2kσk2是II中窗口wkwk中的方差，|w||w|是窗口wkwk中像素的數(shù)量，pkˉpkˉ是待濾波圖像pp在窗口wkwk中的均值，即pkˉ=1|w|∑i∈wkpipkˉ=1|w|∑i∈wkpi。

此外，在計算每個窗口的線性系數(shù)時，我們可以發(fā)現(xiàn)一個像素會被多個窗口包含，也就是說，每個像素都由多個線性函數(shù)所描述。因此，如之前所說，要具體求某一點的輸出值時，只需將所有包含該點的線性函數(shù)值平均即可，如下

qi=1|w|∑k:i∈wk(akIi+bk)=aiˉIi+biˉ

qi=1|w|∑k:i∈wk(akIi+bk)=aiˉIi+biˉ

其中，aiˉ=1|w|∑k∈wiakaiˉ=1|w|∑k∈wiak，biˉ=1|w|∑k∈wibkbiˉ=1|w|∑k∈wibk。

下面給出導(dǎo)向濾波算法的流程，fmeanfmean為一個窗口半徑為rr的均值濾波器（對應(yīng)的窗口大小為2r+12r+1），corr為相關(guān)，var為方差，cov為協(xié)方差。

基于MATLAB的算法實現(xiàn)

下面來在MATLAB中具體實現(xiàn)一下導(dǎo)向濾波算法（代碼來自Dr. Kaiming He，使用時請尊重原作者權(quán)利）。

上述代碼的實現(xiàn)完全遵照上一小節(jié)最后給出的算法流程圖。這里需要略作解釋的地方是函數(shù)boxfilter，它是基于積分圖算法實現(xiàn)的Box Filter。關(guān)于Box Filter，你也可以參考文獻【6】以了解更多。首先來看看它到底做了些什么（因為這個Box Filter 和通常意義上的均值濾波并不完全一樣）。

A =

1 1 1

1 1 1

1 1 1

>> B = boxfilter(A, 1);

>> B

B =

4 6 4

6 9 6

4 6 4

從上述代碼可以看到，當(dāng)參數(shù)r=1時，此時的濾波器窗口是3×33×3。將這樣大小的一個窗口扣在原矩陣中心，剛好可以覆蓋所有矩陣，此時求和為9，即把窗口里覆蓋到的值全部加和。此外當(dāng)把窗口中心挪動到左上角的像素時，因為窗口覆蓋區(qū)域里只有4個數(shù)字，所以結(jié)果為4。所以你可以看出這里的Box Filter只是做了求和處理，并沒有歸一化，所以并不是真正的均值濾波（簡單平滑）。必須結(jié)合后面的一句

mean_I = boxfilter(I, r) ./ N;

1

才算是完成了均值濾波。而這整個過程就相當(dāng)于文獻【6】中介紹的函數(shù)imboxfilt。但是你會發(fā)現(xiàn)它們二者在執(zhí)行的時候最終的結(jié)果（主要是位于圖像四周邊緣的數(shù)值）會有細微的差異。這是因為MATLAB中的imboxfilt函數(shù)在處理位于圖像四周邊緣的像素時，需要虛擬地為原圖像補齊濾波窗口覆蓋但是沒有值的區(qū)域。

下面給出上述boxfilter函數(shù)的實現(xiàn)代碼。

上述代碼基于積分圖實現(xiàn)，如果你對積分圖不是很了解，可以參考其他文獻，這里不再贅述。

下面我們來實驗一下上述導(dǎo)引濾波用于edge-perserving的平滑濾波效果。

I = double(imread('cat.bmp')) / 255;

p = I;

r = 4; % try r=2, 4, or 8

eps = 0.2^2; % try eps=0.1^2, 0.2^2, 0.4^2

O = guidedfilter(I, p, r, eps);

subplot(121), imshow(I);

subplot(122), imshow(O);

執(zhí)行上述代碼，結(jié)果如下所示?？梢娦Ч€是很不錯的。但是我們可以來做一下事后分析，看看導(dǎo)向濾波是如果實現(xiàn)edge-perserving的平滑濾波效果的。當(dāng)I=pI=p時，導(dǎo)向濾波就變成了邊緣保持的濾波操作，此時原來求出的aa和bb的表達式就變成了：

ak=σ2kσ2k+?bk=(1?ak)μk

ak=σk2σk2+?bk=(1?ak)μk

考慮兩種情況：

情況1：高方差區(qū)域，即表示圖像II在窗口wkwk中變化比較大，此時我們有σ2k>>?σk2>>?，于是有ak≈1ak≈1和bk≈0bk≈0。

情況2：平滑區(qū)域（方差不大），即圖像II在窗口wkwk中基本保持固定，此時有σ2k<

也就是說在方差比較大的區(qū)域，保持值不變，在平滑區(qū)域，使用臨近像素平均（也就退化為普通均值濾波）。（這個思想跟文獻【1】里設(shè)計的自適應(yīng)降噪濾波器有異曲同工之妙，當(dāng)然也不完全相同！）

上面的給出的是對灰度圖像進行導(dǎo)向濾波的代碼。你可能會好奇彩色圖像該如何使用導(dǎo)向濾波。一個比較直接的方法就是將引導(dǎo)濾波分別應(yīng)用到RGB三個顏色通道中，然后在組合成結(jié)果圖像。更多細節(jié)可以參考作者原文。不僅如此，在本文開始我們也談到，導(dǎo)向濾波不僅可以應(yīng)用與圖像平滑，還可以用于圖像去霧、圖像Matting等多個領(lǐng)域，限于篇幅我們無法一一詳述，有興趣的讀者可以參考作者原文以了解導(dǎo)向濾波的其他應(yīng)用。

最后，需要說明的是在新版的MATLAB中（R2014及以后），已經(jīng)內(nèi)置了用于導(dǎo)向濾波的函數(shù)imguidedfilter（而這個函數(shù)實現(xiàn)的其實是Fast Guided Filter），也就是說在實際開發(fā)中我們已經(jīng)不再需要編寫上面那樣的代碼，而是只要簡單調(diào)用MATLAB的內(nèi)置函數(shù)就可以了。