o(1)復雜度之雙邊濾波算法的原理、流程、實現及效果。

 一、引言    

     雙邊濾波在圖像處理領域中有着廣泛的應用，比如去噪、去馬賽克、光流估計等等，最近，比較流行的Non-Local算法也可以看成是雙邊濾波的一種擴展。自從Tomasi et al等人提出該算法那一天起，如何快速的實現他，一直是人們討論和研究的焦點之一，在2011年及2012年Kunal N. Chaudhury等人發表的相關論文中，提出了基於三角函數關系的值域核算法，能有效而又准確的實現高效雙邊算法。本文主要對此論文提出的方法加以闡述。

     雙邊濾波的邊緣保持特性主要是通過在卷積的過程中組合空域函數和值域核函數來實現的，典型的核函數為高斯分布函數，如下所示：

                                 （1）

    其中:

            　　　　　　　　　　　　　　（2）

    為歸一化的作用。

    σ_s為空域高斯函數的標准差，σ_r為值域高斯函數的標准差，Ω表示卷積的定義域。 可見，在圖像的平坦區域，f(y)-f(x)的值變化很小，對應的值域權重接近於1，此時空域權重起主要作用，相當於直接對此區域進行高斯模糊，在邊緣區域，f(y)-f(x)會有較大的差異，此時值域系數會下降，從而導致此處整個核函數的分布的下降，而保持了邊緣的細節信息。

     直接的編碼實現上述過程是相當耗時的，其時間復雜度為O(σ_s²) ，因此嚴重的限制住了該算法的推廣和實際使用。不斷有學者提出了解決的辦法，其中Porikli基於一些假定對此過程進行了優化，比如我就實現過其中一種：空域函數為均值函數，值域為任何其他函數，此時可以用直方圖技術進行處理，可減少計算量，但我的實踐表明該算法那速度還是慢，並且效果也不好。

     在2011的論文《Fast O(1) bilateral filtering using trigonometric range kernels》中，作者提出了用Raised cosines函數來逼近高斯值域函數，並利用一些特性把值域函數分解為一些列函數的疊加，從而實現函數的加速。下面我們重點描述下該過程。

二、推導

     1、一些基礎理論和常識。

    (1) Cos函數在[-Pi/2,Pi/2]之間為非負、對稱、在半周期內單調遞增以及且有峰值的函數；

    (2) 歐拉公式：  exp(ix)=cos(x)+isin(x);

    (3) 分配律： exp(a+b)=exp(a)*exp(b)；

    (4) 圖像的動態范圍：[0,T]，比如對於灰度圖像即為[0,255];

     2、一些有用的論證

     (1) 對於式子：

                                                                                      （3）

         其中s是自變量，取值范圍[-T,T]，令γ= Pi / 2T，則γs的值在[-Pi/2,Pi/2]內。此時，可以證明：

                            （4）

       (2) 當N足夠大時，有下式成立：

           　　　　　　　　　　　　　　　　（5）

　　　　如果令ρ=γσ，則上式就變為：

           　　　　　　 （6）

       同樣，上面成立的條件也必須有：

          

       當γs的值在[-Pi/2,Pi/2]時，因此只需要 即可，此時要求 ；

         式6中，最右側部分即為高斯函數，此時說明，可以用 Raised cosines函數來近似的模擬高斯函數，我們用一段matlab函數來驗證該結果：

 

clc;
T=255;
Delta =80;
Gamma = pi/(2*T);
Rho= Gamma * Delta;
Color = ['b','g','r','c','m','y','k'];
x=-T:T;
y1=exp(-x.^2/(2*Delta*Delta));
plot(x,y1,'--b','LineWidth',2);
hold on;
for k=2:7
    y2= cos(Gamma .* x/ (Rho * sqrt(k))).^(k);
    plot(x,y2,Color(k));
end

    

 

      從左圖的曲線分布可見，N=2(綠色），N=3(紅色）兩條曲線明顯不符合我們的定義域的要 求，分別出現了非單調遞增和負值得情況。之后隨着N的不斷增大，曲線越來越接近高斯分布曲線（藍色曲線)。 這從實際的角度證明了公式6的正確性。

　　3、推導

      將公式（4）帶入公式（6）中，得到：

                （7）

       將公式（7）帶入原始的雙邊濾波的公式（1）中，有：

    （8）

其中：

            

     式（8）的第三行利用了前面基礎理論中的第三條：分配率。注意式中的積分是對y積分，因此可以把f(x)相關部分提取到積分符號的外圍。

     同樣，對於η,可以推導得到：

                （9）

     注意（8)和(9)兩式中后的后半部分，比如 這個，可以看出這個實際上就是對cos(β)這幅圖像進行高斯模糊，而高斯模糊可以通過FFT或者回溯算法快速O(1)實現，這樣兩式8/9兩式就分別轉換為對f(y)cos(β)、f(y)sin(β)、cos(β)以及sin(β)圖像進行一系列高斯模糊的過程了。

     至此，所有的推導工作完成，那么我們在理一下算法的執行步驟吧：

     （1）：已知條件：輸入圖像f(x),動態范圍[-T,T],空域和值域方差σ_s、σ_r。

     （2）：設定γ = Pi / 2T,ρ=γσ_r,N=1/(γσ_r*γσ_r) ;

      (3)  : for (0≤n≤N),獲取f(y)cos(β)、f(y)sin(β)、cos(β)以及sin(β)所對應的圖像數據（浮點類型）；

     （4）：利用O(1)高斯模糊算法對上述四個圖像數據進行標准差為σ_s的高斯模糊並累計。

     （5）：對累加后的數據進行除法操作，獲得最終的結果圖像。     

    注意：式8和式9中的乘法最后會有虛部的數據出現，在處理時可以直接丟棄掉。

三、實現及效果

    以上算法在論文Fast O(1) bilateral filtering using trigonometric range kernels中有着較為詳細的論述，論文中還給出了JAVA代碼實現的鏈接，但是該鏈接已經失效，需要JAVA代碼做參考的可從此處下載：BilateralFilter-src.rar，其中的BilateralFilter_.jar可在ImageJ中作為插件加載，而這篇論文的對應代碼在解壓后的bilateralfilterinstant文件夾中。注意，這個ImageJ的插件寫的似乎有問題，運行時點plugins-->BilateralFilters-->Bilateral Filter Instant 后彈出的對話框中，不要勾選多線程才能對輸入的任意參數進行處理，否則圖像無任何反映。

     在第三步和第四步的處理，N+1次循環之間時沒有任何關系的，因此，只要內存許可，各循環之間可以並行的執行，這對於現在的2核和4核的CPU的有一定的意義，不過相比GPU來說，可能意義更大吧。

     按上述過程編制代碼，測試效果測試如下:

    （1） 去除噪音

      

                          lena圖增加標准差為20的高斯噪音　　　　　　　　　　　　　　　　　　　                　使用σ_s=20、σ_r=40雙邊濾波后的結果

     （2）普通圖像的邊緣保持結果

    

　　　　　　　　　　　　 帶有噪音的美女圖1　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　使用σ_s=10、σ_r=35雙邊濾波后的結果   

 

  

　　 　　帶有噪音的美女圖2　　　　　　　　　　　　　　　　　　　　使用σ_s=15、σ_r=25雙邊濾波后的結果     　　　　　　 　使用σ_s=10、σ_r=35雙邊濾波后的結果   

     美女1 細膩的頭發在濾波后依然有着飄然的感覺，而面部的斑點也已經悄然消失。美女2的面部也得到了恰當的護理，而其湛藍的雙眼依舊是那么有神。

四、效率分析及進一步優化

     很明顯，上述算法的執行時間直接依賴於N的大小，而從相關推導公式中看，N的值直接取決於T和σ_r的大小，T越大,N越大，σ_r越大，N越小。當T固定時，比如為255時，取不同的σ_r對應的N值如下表所示：

σr	5	10	20	30	40	60	80	1000
N	1053	263	66	29	16	7	4	3

    上表中，可以發現，當σ_r小於20時，所需要的循環次數N極具增加，當N超過100時，可以認為這個算法已經不再具有優化的意義了，可能比原始的算法還慢。這個從原理上也很好解釋。當σ_r比較小時，高斯函數的曲線在中心線附近急劇下降，從而需要更多的三角函數來逼近他。

     因此，進一步的優化需要從T的取值以及N的方面予以考慮。

     我們觀察到，在對值域高斯核函數進行計算時，使用的是f(y)-f(x),而不是f(y),通常情況下，一副圖像的f(y)中所有值的最大值可能為255，但是某個區域內的f(y)-f(x)在全圖中的值卻不一定為255，如果這個值小於255，那么在上述算法中使用的T值就可以縮小，從而減少N的大小，比如當σ_r為10時，T取255，由上表可知N=263，而如果T為210，則大約只需要 263*（210/255）^2 = 178 次循環，幾乎減少了100次。

     那么一副圖像的實際的T值打開是什么個范疇呢，使用標准的Lena圖做了個測試，在不同的取樣半徑下獲得的結果如下表：

r	1	3	5	10	15	20	30
T	153	205	208	210	211	215	215

    這說明在很多情況下T值確實小於255，即使取樣半徑比較大。

    這里則涉及到一個平衡問題。計算實際的T值一般情況下會獲得小於255的結果，這有利於減小N，從而降低程序的執行時間，但是這個計算的過程本身也需要時間，如果這個時間大於其帶來的好處，則這個改進就是退步。幸好，在很久以前，關於指定半徑內的最大值算法就已經有了O(1)的快速算法， 其執行時間一般要小於進行一次本例中這種循環的時間。所以這個改進是值得的。

     那么我在看看小 σ_r時大N的問題，當σ_r比較小時，我們觀察其分布曲線，如下圖：

 

       　　　　　　　　　　　　　　σ_r = [1 3 5 10 15]時的曲線

     由上圖可以看出小σ_r時，曲線在中心線附近迅速衰減，理論表明這個距離為[-3σ_r,3σ_r],在此之外的值可以忽略不計，因此，那些對最終結果沒有什么貢獻的循環就完全可以舍棄，這部分的理論推導可以詳見論文 Acceleration of the shiftable O(1) algorithm for bilateral filtering and non-local means 。

      我們知道，Non-Local算法在很大程度是雙邊模糊的擴展，只是其值域的相似度函數更加復雜，不是簡單的f(y)-(f(x)那么簡單了，而是和f(y)和f(x)的領域有關，因此直接的Non-Local實現理論上比雙邊濾波還要耗時，上面介紹的這種優化方式在后面這篇論文里提到也是可以用於Non-Local的，有興趣的朋友可以自己去研究下。

五：小結和展望

      可以看到，本文的這種優化方式實際上是利用Cos函數去逼近高斯函數，在代碼層次上，需要(N+1)*4此高斯模糊，而由上面相關表格可以看到，N的數值一般都在10以上，因此，至少要執行44次高斯模糊，這還不包括獲取需要高斯模糊的數據部分以及最后的濾波結果獲取。即使高斯模糊再高效，比如對於600*400的彩色圖，5ms足也，那么執行雙邊模糊保守估計也要400ms左右，因此這個算法說實在的效率還是不行。

      這兩篇文章分別是2011年及2012年發表的，應該是代表着目前比較先進的技術，我在網上經常看到有人說雙邊濾波可是實時，實在是不曉得那些高人用的是什么理論，抑或是什么超級機器。

      同我之間的一篇博文中雙指數邊緣平滑濾波器用於磨皮算法的嘗試 提到的Beeps邊緣保留算法相比，這里的速度就要慢很多了，而兩者的效果相比基本上差不多，所以實在很糾結。

      希望有更好的方法用於該算法。

      我用C編寫了一個標准DLL，其C#調用示例可見: http://files.cnblogs.com/Imageshop/BilateralFilterTest.rar

 

 

 

 

*********************************作者： laviewpbt   時間： 2013.11.4   聯系QQ:  33184777  轉載請保留本行信息************************