現在從一個最簡單的情形來開始我們的討論。假設有一個原始圖像 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中的像素數目。也就是說，在簡單平滑中，以像素ii為中心的一個窗口ww中的像素具有等同的權值。但是在高斯平滑中，權值WijWij將服從二維的高斯分布，結果導致離像素ii更接近的像素將具有更高的權重，反之離像素ii較遠的像素則具有更小的權重。

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

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

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

算法原理

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

qi=∑j∈wiWij(I)?pj

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

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

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

現在已知的是II和pp，要求的是qq。而如果能求得參數aa和bb，顯然就能通過II和qq之間的線性關系來求出II。另一方面，還可以知道pp是qq受到噪聲污染而產生的退化圖像，假設噪聲是nn，則有qi=pi?niqi=pi?ni。根據無約束圖像復原的方法（可以參考【4】），這時可以設定最優化目標為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

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

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

求解上述最優化問題（跟最小二乘法的推導過程一致），便會得到：

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

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

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

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。

下面給出導向濾波算法的流程，fmeanfmean為一個窗口半徑為rr的均值濾波器（對應的窗口大小為2r+12r+1），corr為相關，var為方差，cov為協方差。

基于MATLAB的算法實現

下面來在MATLAB中具體實現一下導向濾波算法（代碼來自Dr. Kaiming He，使用時請尊重原作者權利）。

上述代碼的實現完全遵照上一小節最后給出的算法流程圖。這里需要略作解釋的地方是函數boxfilter，它是基于積分圖算法實現的Box Filter。關于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

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

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

1

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

下面給出上述boxfilter函數的實現代碼。

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

下面我們來實驗一下上述導引濾波用于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);

執行上述代碼，結果如下所示。可見效果還是很不錯的。但是我們可以來做一下事后分析，看看導向濾波是如果實現edge-perserving的平滑濾波效果的。當I=pI=p時，導向濾波就變成了邊緣保持的濾波操作，此時原來求出的aa和bb的表達式就變成了：

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

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

考慮兩種情況：

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

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

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

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

最后，需要說明的是在新版的MATLAB中（R2014及以后），已經內置了用于導向濾波的函數imguidedfilter（而這個函數實現的其實是Fast Guided Filter），也就是說在實際開發中我們已經不再需要編寫上面那樣的代碼，而是只要簡單調用MATLAB的內置函數就可以了。