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

具體而言，在簡(jiǎn)單平滑中，圖像qq中第ii個(gè)像素是由圖像pp中以第ii個(gè)像素為中心的一個(gè)窗口ww中的所有像素取平均而得來(lái)的，即

qi=∑j∈wiWij?pj

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

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

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

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

算法原理

導(dǎo)向?yàn)V波之所以叫這個(gè)名字，因?yàn)樵谒惴蚣苤校獙?duì)pp進(jìn)行濾波而得到qq，還得需要一個(gè)引導(dǎo)圖像II。此時(shí)，濾波器的數(shù)學(xué)公式為

qi=∑j∈wiWij(I)?pj

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

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

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

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

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

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

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

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。

此外，在計(jì)算每個(gè)窗口的線性系數(shù)時(shí)，我們可以發(fā)現(xiàn)一個(gè)像素會(huì)被多個(gè)窗口包含，也就是說(shuō)，每個(gè)像素都由多個(gè)線性函數(shù)所描述。因此，如之前所說(shuō)，要具體求某一點(diǎn)的輸出值時(shí)，只需將所有包含該點(diǎn)的線性函數(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)向?yàn)V波算法的流程，fmeanfmean為一個(gè)窗口半徑為rr的均值濾波器（對(duì)應(yīng)的窗口大小為2r+12r+1），corr為相關(guān)，var為方差，cov為協(xié)方差。

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

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

上述代碼的實(shí)現(xiàn)完全遵照上一小節(jié)最后給出的算法流程圖。這里需要略作解釋的地方是函數(shù)boxfilter，它是基于積分圖算法實(shí)現(xiàn)的Box Filter。關(guān)于Box Filter，你也可以參考文獻(xiàn)【6】以了解更多。首先來(lái)看看它到底做了些什么（因?yàn)檫@個(gè)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時(shí)，此時(shí)的濾波器窗口是3×33×3。將這樣大小的一個(gè)窗口扣在原矩陣中心，剛好可以覆蓋所有矩陣，此時(shí)求和為9，即把窗口里覆蓋到的值全部加和。此外當(dāng)把窗口中心挪動(dòng)到左上角的像素時(shí)，因?yàn)榇翱诟采w區(qū)域里只有4個(gè)數(shù)字，所以結(jié)果為4。所以你可以看出這里的Box Filter只是做了求和處理，并沒(méi)有歸一化，所以并不是真正的均值濾波（簡(jiǎn)單平滑）。必須結(jié)合后面的一句

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

1

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

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

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

下面我們來(lái)實(shí)驗(yà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é)果如下所示。可見(jiàn)效果還是很不錯(cuò)的。但是我們可以來(lái)做一下事后分析，看看導(dǎo)向?yàn)V波是如果實(shí)現(xiàn)edge-perserving的平滑濾波效果的。當(dāng)I=pI=p時(shí)，導(dǎo)向?yàn)V波就變成了邊緣保持的濾波操作，此時(shí)原來(lái)求出的aa和bb的表達(dá)式就變成了：

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

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

考慮兩種情況：

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

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

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

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

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