導向濾波(Guided Filter)的解析與實現
現在從一個最簡單的情形來開始我們的討論。假設有一個原始影象 pp較遠的畫素則具有更小的權重。
無論是簡單平滑，還是高斯平滑，它們都有一個共同的弱點，即它們都屬於各向同性濾波。我們都知道，一幅自然的影象可以被看成是有（過渡平緩的，也就是梯度較小）區域和（過渡尖銳的，也就是梯度較大）邊緣（也包括影象的紋理、細節等）共同組成的。噪聲是影響影象質量的不利因素，我們希望將其濾除。噪聲的特點通常是以其為中心的各個方向上梯度都較大而且相差不多。邊緣則不同，邊緣相比於區域也會出現梯度的越變，但是邊緣只有在其法向方向上才會出現較大的梯度，而在切向方向上梯度較小。
因此，對於各向同性濾波（例如簡單平滑或高斯平滑）而言，它們對待噪聲和邊緣資訊都採取一直的態度。結果，噪聲被磨平的同時，影象中具有重要地位的邊緣、紋理和細節也同時被抹平了。這是我們所不希望看到的。研究人員已經提出了很多Edge-perserving的影象降噪（平滑）演算法，例如雙邊濾波、自適應（維納）平滑濾波（請參見文獻【1】）、基於PM方程的各向異性濾波以及基於TV-norm的降噪演算法等。本文將考慮在文獻【2】中提出的另外一種Edge-perserving的影象濾波（平滑）演算法——導向濾波（Guided Filter）。當然，通過閱讀文獻【2】，我們也知道導向濾波的應用不止有Edge-perserving的影象平滑，還包括影象去霧、影象Matting等等。
歡迎關注白馬負金羈的部落格 http://blog.csdn.net/baimafujinji，為保證公式、圖表得以正確顯示，強烈建議你從該地址上檢視原版博文。本部落格主要關注方向包括：數字影象處理、演算法設計與分析、資料結構、機器學習、資料探勘、統計分析方法、自然語言處理。
演算法原理
導向濾波之所以叫這個名字，因為在演算法框架中，要對pp本身的時候，導向濾波就變成了一個Edge-perserving的濾波器，因此可以用於影象的平滑降噪。
導向濾波濾波的示意圖如下所示（圖片來自作者原文【2】）。注意這也是導向濾波所依賴的一個重要假設——導向濾波器在導向影象II中，儘管丙處的梯度仍然大於乙處的梯度，進而大於甲處的梯度，但是倍數關係可能會扭曲。例如，丙處的梯度是乙處的1.5倍，而是甲處的6倍（即乙處的梯度是甲處的4倍），這時你就會想象，甲處是區域，而乙處和丙處就變成了邊緣。可見非線性關係會使得引導影象對於邊緣和區域的指示作用發生錯亂。
現在已知的是II），corr為相關，var為方差，cov為協方差。
基於MATLAB的演算法實現
下面來在MATLAB中具體實現一下導向濾波演算法（程式碼來自Dr. Kaiming He，使用時請尊重原作者權利）。
function q = guidedfilter(I, p, r, eps)
% - guidance image: I (should be a gray-scale/single channel image)
% - filtering input image: p (should be a gray-scale/single channel image)
% - local window radius: r - 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
上述程式碼的實現完全遵照上一小節最後給出的演算法流程圖。這裡需要略作解釋的地方是函式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- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
從上述程式碼可以看到，當引數r=1時，此時的濾波器視窗是3×33×3。將這樣大小的一個視窗扣在原矩陣中心，剛好可以覆蓋所有矩陣，此時求和為9，即把窗口裡覆蓋到的值全部加和。此外當把視窗中心挪動到左上角的畫素時，因為視窗覆蓋區域裡只有4個數字，所以結果為4。所以你可以看出這裡的Box Filter只是做了求和處理，並沒有歸一化，所以並不是真正的均值濾波（簡單平滑）。必須結合後面的一句
mean_I = boxfilter(I, r) ./ N;- 1
才算是完成了均值濾波。而這整個過程就相當於文獻【6】中介紹的函式imboxfilt。但是你會發現它們二者在執行的時候最終的結果（主要是位於影象四周邊緣的數值）會有細微的差異。這是因為MATLAB中的imboxfilt函式在處理位於影象四周邊緣的畫素時，需要虛擬地為原影象補齊濾波視窗覆蓋但是沒有值的區域。
下面給出上述boxfilter函式的實現程式碼。
function imDst = boxfilter(imSrc, r)
% BOXFILTER O(1) time box filtering using cumulative sum
%
% - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
% - Running time independent of r;
% - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
% - But much faster.
[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));
%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);
%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
上述程式碼基於積分圖實現，如果你對積分圖不是很瞭解，可以參考文獻【7】，這裡不再贅述。
下面我們來實驗一下上述導引濾波用於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);- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
執行上述程式碼，結果如下所示。可見效果還是很不錯的。但是我們可以來做一下事後分析，看看導向濾波是如果實現edge-perserving的平滑濾波效果的。當I=pI=p考慮兩種情況：
- 情況1：高方差區域，即表示影象II中基本保持固定，此時有σ2k<<ϵ