雙邊濾波

雙邊濾波很有名，使用廣泛，簡單的說就是一種同時考慮了畫素空間差異與強度差異的濾波器，因此具有保持影象邊緣的特性。

先看看我們熟悉的高斯濾波器

其中W是權重，i和j是畫素索引，K是歸一化常量。公式中可以看出，權重只和畫素之間的空間距離有關係，無論影象的內容是什麼，都有相同的濾波效果。

再來看看雙邊濾波器，它只是在原有高斯函式的基礎上加了一項，如下

其中 I 是畫素的強度值，所以在強度差距大的地方（邊緣），權重會減小，濾波效應也就變小。總體而言，在畫素強度變換不大的區域，雙邊濾波有類似於高斯濾波的效果，而在影象邊緣等強度梯度較大的地方，可以保持梯度。

引導濾波

引導濾波是近三年才出現的濾波技術，知道的人還不多。它與雙邊濾波最大的相似之處，就是同樣具有保持邊緣特性。在引導濾波的定義中，用到了區域性線性模型，至於該模型，可以暫時用下圖簡單的理解

該模型認為，某函式上一點與其鄰近部分的點成線性關係，一個複雜的函式就可以用很多區域性的線性函式來表示，當需要求該函式上某一點的值時，只需計算所有包含該點的線性函式的值並做平均即可。這種模型，在表示非解析函式上，非常有用。

同理，我們可以認為影象是一個二維函式，而且沒法寫出解析表示式，因此我們假設該函式的輸出與輸入在一個二維視窗內滿足線性關係，如下

其中，q是輸出畫素的值，I是輸入影象的值，i和k是畫素索引，a和b是當視窗中心位於k時該線性函式的係數。其實，輸入影象不一定是待濾波的影象本身，也可以是其他影象即引導影象，這也是為何稱為引導濾波的原因。對上式兩邊取梯度，可以得到

即當輸入影象I有梯度時，輸出q也有類似的梯度，現在可以解釋為什麼引導濾波有邊緣保持特性了。

下一步是求出線性函式的係數，也就是線性迴歸，即希望擬合函式的輸出值與真實值p之間的差距最小，也就是讓下式最小

這裡p只能是待濾波影象，並不像I那樣可以是其他影象。同時，a之前的係數（以後都寫為e）用於防止求得的a過大，也是調節濾波器濾波效果的重要引數。通過最小二乘法，我們可以得到

其中，是I在視窗w_k中的平均值，是I在視窗w_k中的方差，是視窗w_k中畫素的數量，是待濾波影象p在視窗w_k中的均值。

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

這裡，w_k是所有包含畫素i的視窗，k是其中心位置。

當把引導濾波用作邊緣保持濾波器時，往往有 I = p ，如果e=0，顯然a=1, b=0是E(a,b)為最小值的解，從上式可以看出，這時的濾波器沒有任何作用，將輸入原封不動的輸出。如果e>0，在畫素強度變化小的區域（或單色區域），有a近似於（或等於）0，而b近似於（或等於），即做了一個加權均值濾波；而在變化大的區域，a近似於1，b近似於0，對影象的濾波效果很弱，有助於保持邊緣。而e的作用就是界定什麼是變化大，什麼是變化小。在視窗大小不變的情況下，隨著e的增大，濾波效果越明顯。

在濾波效果上，引導濾波和雙邊濾波差不多，在一些細節上，引導濾波較好。引導濾波最大的優勢在於，可以寫出時間複雜度與視窗大小無關的演算法（打算在之後的文章中討論），因此在使用大視窗處理圖片時，其效率更高。

關於引導濾波更多的討論和應用，可以參看下面的論文

上一篇文章已經說了引導濾波的基本理論，而且我們也知道引導濾波可以寫出時間複雜度與視窗大小無關的演算法，現在就來使用C++並藉助OpenCV實現這一演算法。

實現這種演算法的關鍵思想是盒式濾波（box filter），而且必須是通過積分圖來實現的盒式濾波，否則不可能與視窗大小無關，好在OpenCV的boxFilter函式滿足這個要求。

再看看引導濾波的公式

先計算a_k的分子，Ip 在視窗w_k中的和，再除以視窗中畫素的個數，剛好就是盒式濾波，因此我們可以將輸入的引導影象 I 和濾波影象 p 相乘，並對相乘後的影象做box filtering，即得第一項的結果。後面的和分別為 I 和 p 在視窗w_k中均值，因此分別對 I 和 p 進行box filtering，再將box filtering之後的結果相乘即可。實際上，a_k的分子就是 Ip 在視窗w_k中的協方差。

接下來計算a_k的分母部分。是引導圖 I 在視窗w_k中的方差，學過概率論與數理統計的朋友應該知道，方差和期望（均值）之間是有關係的，如下式

因此在計算 I 的方差時，我們可以先計算 I*I 的均值，再減去 I 均值的平方即的平方。在方法上，計算 I*I 的均值和計算 Ip 的均值是一樣的。最後，對計算出來的方差影象，加上常量e（每個元素都加e），分母就計算完了，自然，a_k在所有視窗中的值也就得到了。b_k的計算太簡單了，大家都懂的。

注意，我們的計算都是對整個影象的，以影象為單位進行計算，所以最後算出的也是兩張圖，a_k的圖（左邊）和b_k的圖（右邊），如下

在圖中可以看到，在邊緣部分或變化劇烈的部分，a的值接近於1（白色），b的值接近為0（黑色），而在變化平坦的區域，a的值接近0（黑色），b的值為平坦區域畫素的均值。這與上一篇文章中所說的規律是一致的。

下面看第二個公式

輸出值q又與兩個均值有關，分別為a和b在視窗w_i中的均值（不是w_k），所以還是box filtering，我們將上一步得到兩個影象都進行盒式濾波，得到兩個新圖：a_i和b_i，然後用a_i乘以引導影象 I ，再加上b_i，即得最終濾波之後的輸出，如下（左邊為原圖，右邊為濾波之後的影象，其中濾波視窗半徑為8，e的值為500）：

下面是整個演算法的程式碼，僅供參考

1. void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon)
2. {
3. CV_Assert(radius >= 2 && epsilon > 0);
4. CV_Assert(source.data != NULL && source.channels() == 1);
5. CV_Assert(guided_image.channels() == 1);
6. CV_Assert(source.rows == guided_image.rows && source.cols == guided_image.cols);
7. Mat guided;
8. if (guided_image.data == source.data)
9. {
10. //make a copy
11. guided_image.copyTo(guided);
12. }
13. else
14. {
15. guided = guided_image;
16. }
17. //將輸入擴充套件為32位浮點型，以便以後做乘法
18. Mat source_32f, guided_32f;
19. makeDepth32f(source, source_32f);
20. makeDepth32f(guided, guided_32f);
21. //計算I*p和I*I
22. Mat mat_Ip, mat_I2;
23. multiply(guided_32f, source_32f, mat_Ip);
24. multiply(guided_32f, guided_32f, mat_I2);
25. //計算各種均值
26. Mat mean_p, mean_I, mean_Ip, mean_I2;
27. Size win_size(2*radius + 1, 2*radius + 1);
28. boxFilter(source_32f, mean_p, CV_32F, win_size);
29. boxFilter(guided_32f, mean_I, CV_32F, win_size);
30. boxFilter(mat_Ip, mean_Ip, CV_32F, win_size);
31. boxFilter(mat_I2, mean_I2, CV_32F, win_size);
32. //計算Ip的協方差和I的方差
33. Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
34. Mat var_I = mean_I2 - mean_I.mul(mean_I);
35. var_I += epsilon;
36. //求a和b
37. Mat a, b;
38. divide(cov_Ip, var_I, a);
39. b = mean_p - a.mul(mean_I);
40. //對包含畫素i的所有a、b做平均
41. Mat mean_a, mean_b;
42. boxFilter(a, mean_a, CV_32F, win_size);
43. boxFilter(b, mean_b, CV_32F, win_size);
44. //計算輸出 (depth == CV_32F)
45. output = mean_a.mul(guided_32f) + mean_b;
46. }

void guidedFilter(Mat& source, Mat& guided_image, Mat& output, int radius, float epsilon)
{
	CV_Assert(radius >= 2 && epsilon > 0);
	CV_Assert(source.data != NULL && source.channels() == 1);
	CV_Assert(guided_image.channels() == 1);
	CV_Assert(source.rows == guided_image.rows && source.cols == guided_image.cols);

	Mat guided;
	if (guided_image.data == source.data)
	{
		//make a copy
		guided_image.copyTo(guided);
	}
	else
	{
		guided = guided_image;
	}

	//將輸入擴充套件為32位浮點型，以便以後做乘法
	Mat source_32f, guided_32f;
	makeDepth32f(source, source_32f);
	makeDepth32f(guided, guided_32f);

	//計算I*p和I*I
	Mat mat_Ip, mat_I2;
	multiply(guided_32f, source_32f, mat_Ip);
	multiply(guided_32f, guided_32f, mat_I2);

	//計算各種均值
	Mat mean_p, mean_I, mean_Ip, mean_I2;
	Size win_size(2*radius + 1, 2*radius + 1);
	boxFilter(source_32f, mean_p, CV_32F, win_size);
	boxFilter(guided_32f, mean_I, CV_32F, win_size);
	boxFilter(mat_Ip, mean_Ip, CV_32F, win_size);
	boxFilter(mat_I2, mean_I2, CV_32F, win_size);

	//計算Ip的協方差和I的方差
	Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
	Mat var_I = mean_I2 - mean_I.mul(mean_I);
	var_I += epsilon;

	//求a和b
	Mat a, b;
	divide(cov_Ip, var_I, a);
	b = mean_p - a.mul(mean_I);

	//對包含畫素i的所有a、b做平均
	Mat mean_a, mean_b;
	boxFilter(a, mean_a, CV_32F, win_size);
	boxFilter(b, mean_b, CV_32F, win_size);

	//計算輸出 (depth == CV_32F)
	output = mean_a.mul(guided_32f) + mean_b;
}

1. void makeDepth32f(Mat& source, Mat& output)
2. {
3. if (source.depth() != CV_32F ) > FLT_EPSILON)
4. source.convertTo(output, CV_32F);
5. else
6. output = source;
7. }

void makeDepth32f(Mat& source, Mat& output)
{
	if (source.depth() != CV_32F ) > FLT_EPSILON)
		source.convertTo(output, CV_32F);
	else
		output = source;
}