閱讀Real-Time O(1) Bilateral Filtering 一文的相關感覺。

　　研究雙邊濾波有很長一段時間了，最近看了一篇Real-Time O(1) Bilateral Filtering的論文，標題很吸引人，就研讀了一番，通過幾天的攻讀，基本已理解其思想，現將這一過程作一簡單的小結。

     論文大於10MB，沒法上傳至博客園，能夠在這個連接下載：http://www.cs.cityu.edu.hk/~qiyang/publications/cvpr-09-qingxiong-yang.pdf。

     首先，先給出一個我本身的結論：這篇文章無啥新意，主要的算法思想都來自於另一篇論文，Fast Bilateral Filtering for the Display of High-Dynamic-Range Images，並且文中的部分實驗結果我認爲存在較大的水分，可是，其中提到的算法仍是比較快的。

     論文中對雙邊模糊的優化思路大概是這樣的：

     對於雙邊模糊，離散化後的表達式大概以下所示：

      

     f(s)是空域核函數，f(r)是值域核函數，  難以直接優化上式的緣由是 f(r)的存在。  

     論文中提出一種思路，若是上式中固定I(X)的值，則對於每個不一樣的I(y)值，上式的分子就至關於對f_R(I(x),I(y))*I(y)進行空域核卷積運算，分母則是對f_R(I(x),I(y))進行空域核卷積元算，而這種卷積運算有着快速的算法。 這樣，咱們在圖像的值域範圍內選定若干有表明性的I(X)值，分別進行卷積，而後對於圖像中的其餘的像素值，進行線性插值獲得。

     算法的主要貢獻也就在這裏，而這個想法是從Fast Bilateral Filtering for the Display of High-Dynamic-Range Images一文中獲得的，而且在此文中還提到了進行subsampleing進行進一步的優化，即這些抽樣卷積能夠在原圖的小圖中進行，而後最後的結果在原圖中經過雙線性插值獲取。

      關於直接採樣而後插值的算法源代碼能夠參考：http://files.cnblogs.com/Imageshop/qx_constant_time_bilateral_filter.rar

      下面爲其主要的實現代碼：

 1 int qx_constant_time_bilateral_filter::bilateral_filter(unsigned char **image_filtered,unsigned char **image,double sigma_range,unsigned char **texture)
 2 {
 3     unsigned char image_min,image_max; 
 4     int y,x,jk_0,jk_1;
 5     if(sigma_range>QX_DEF_THRESHOLD_ZERO) 
 6     {
 7         m_sigma_range=sigma_range;
 8         color_weighted_table_update(m_table,m_sigma_range*QX_DEF_CTBF_INTENSITY_RANGE,QX_DEF_CTBF_INTENSITY_RANGE);
 9     }
10     qx_timer timer;
11     timer.start();
12     if(texture==NULL)
13     {
14         vec_min_val(image_min,image[0],m_h*m_w);
15         vec_max_val(image_max,image[0],m_h*m_w);
16     }
17     else
18     {
19         vec_min_val(image_min,texture[0],m_h*m_w);
20         vec_max_val(image_max,texture[0],m_h*m_w);
21     }
22     m_nr_scale=qx_max(1,int(double(image_max-image_min)/(255*m_sigma_range)+0.5));
23     //printf("[qx_max,qx_min]:[%5.5f,%5.5f]\n",(float)image_max,(float)image_min);
24     //printf("[sigma_range: %1.3f]\n",m_sigma_range);
25     //printf("[nr_scale: %d]\n",m_nr_scale);
26     m_grayscale[0]=(double)image_min;
27     m_grayscale[m_nr_scale-1]=(double)image_max;
28     double delta_scale=double(image_max-image_min)/(m_nr_scale-1);
29     for(int i=1;i<m_nr_scale-1;i++) m_grayscale[i]=(double)image_min+delta_scale*i;
30     for(int i=0;i<m_nr_scale;i++)
31     {
32         double **jk;
33         if(i==0)
34         {
35             jk_0=0;
36             jk_1=1;
37             jk=m_jk[jk_0];
38         }
39         else 
40             jk=m_jk[jk_1];
41         for(y=0;y<m_h;y++)
42         {
43             for(x=0;x<m_w;x++)
44             {
45                 int index;
46                 if(texture==NULL) index=int(abs(m_grayscale[i]-image[y][x])+0.5f);
47                 else index=int(abs(m_grayscale[i]-texture[y][x])+0.5f); /*cross/joint bilateral filtering*/
48                 jk[y][x]=m_table[index]*image[y][x];
49                 m_wk[y][x]=m_table[index];
50             }
51         }
52         if(m_spatial_filter==QX_DEF_CTBF_BOX_BILATERAL_FILTER)
53         {
54             boxcar_sliding_window(jk,jk,m_box,m_h,m_w,m_radius);
55             boxcar_sliding_window(m_wk,m_wk,m_box,m_h,m_w,m_radius);
56         }
57         else if(m_spatial_filter==QX_DEF_CTBF_GAUSSIAN_BILATERAL_FILTER)
58         {
59             gaussian_recursive(jk,m_box,m_sigma_spatial*qx_min(m_h,m_w),0,m_h,m_w);
60             gaussian_recursive(m_wk,m_box,m_sigma_spatial*qx_min(m_h,m_w),0,m_h,m_w);
61         }
62         for(y=0;y<m_h;y++)
63         {
64             for(x=0;x<m_w;x++)
65             {
66                 jk[y][x]/=m_wk[y][x];
67             }
68         }
69         //image_display(jk,m_h,m_w);
70         if(i>0)
71         {
72             for(y=0;y<m_h;y++)
73             {
74                 for(x=0;x<m_w;x++)
75                 {
76                     double kf;
77                     if(texture==NULL) kf=double(image[y][x]-image_min)/delta_scale;
78                     else kf=double(texture[y][x]-image_min)/delta_scale; /*cross/joint bilateral filtering*/
79                     int k=int(kf); 
80                     if(k==(i-1))
81                     {
82                         double alpha=(k+1)-kf;
83                         image_filtered[y][x]=(unsigned char)qx_min(qx_max(alpha*m_jk[jk_0][y][x]+(1.f-alpha)*m_jk[jk_1][y][x],0.f)+0.5f,255.f);
84                     }
85                     else if(k==i&&i==(m_nr_scale-1)) image_filtered[y][x]=(unsigned char)(m_jk[jk_1][y][x]+0.5f);
86                 }
87             }
88             jk_1=jk_0;
89             jk_0=(jk_0+1)%2;
90         }
91     }
92     //timer.time_display("bilateral filter");
93     return(0);
94 }

 　　我這裏對其中的代碼進行簡單的描述：

      一、第1三、14行是獲取圖像的動態範圍，即具備最大亮度和最小亮度的像素值。

      二、 第22行的m_nr_scale是計算在原圖中的取樣數。26-29行中的m_grayscale是用來記錄取樣點的值的，好比若是動態範圍是[0,255],取樣數,5，則m_grayscale的值分別爲0、6四、12八、19二、255， 即對式（1）中的I(x)先固定爲這5個值，計算式（1）的結果。

      三、第32到第40行直接的這些代碼實際上是爲了節省內存的，由於若是取樣點爲5，那麼就須要5*2倍原圖大小內存的空間來存儲取樣點的卷積值，可是若是我按取樣點的大小順序計算，那麼每計算一個取樣點後（第一個除外，這就是70行的判斷），就能夠把原圖中夾子於這個取樣點及這個取樣點以前那個取樣數據之間的像素的結果值經過二者之間的線性插值獲取。這種方案就能夠只須要2*2倍原圖大小的內存。可是這種方案就涉及到插值的順序，32到40就是處理這樣的過程的，實際的寫法你能夠有不少種，上面的代碼寫的很爛的。

      四、52到61之間的代碼是看空域你是用什麼類型的卷積函數，這裏可使用任意的其餘的卷積函數，而至於的卷積函數也能夠時任意的，這個能夠參考代碼中的color_weighted_table_update函數內的代碼。

       五、第72到87行的代碼就是對其飛取樣點的數據進行插值的過程，注意一些邊緣的處理過程。

    用插值+SubSampleing的代碼能夠從這裏下載：http://files.cnblogs.com/Imageshop/qx_constant_time_bilateral_filter%28%E5%A2%9E%E5%BC%BA%E7%89%88%29.rar

    試驗結果（SigmaS=10,SigmaR=30,使用高斯卷積核函數)：

       

　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　 　上述算法的結果　　　　　　　　　　　　　　　　　　　　  標準的結果

　　上述的取樣數是按照第22行的m_nr_scale設置的，可見，視覺上彷佛二者之間沒有什麼差異。

     按照m_nr_scale的計算方式，若是SigmaR比較小，m_nr_scale值也會比較大，對於一些工程應用，每每SigmaR就是要取比較小的值才能保護住邊緣。所以，咱們嘗試修改m_nr_scale的值，實際的測試代表，將m_nr_scale的值再該小一半，也能得到很爲理想的效果，而速度確能夠提升一倍。

     另外，上述代碼還應對m_nr_scale的最小值作個限制，m_nr_scale必須至少大於等於2的，不然沒法插值的。

     在速度上，使用這種方式加上一些其餘的優化技巧，SigmaR=30（SigmaS對速度沒有影響）時，一副640*480的彩色圖像，在I3的筆記本上耗時約爲75ms(值域使用均值模糊）、125ms(值域使用高斯函數)。

     論文中提升的下采樣技術進行速度的提高，個人見解看狀況取捨。我本身也進行了編程，得出的結論是：

     一、下采樣的係數越小，結果和準確值誤差越大，而且此時由於下采樣形成高斯濾波或者均值濾波的加快已經在整個耗時裏佔用的比例不大了，此時主要的矛盾是最後的雙線性插值以及線性插值了，所以，整體時間上無明顯提高。所以，我建議採樣倍數不要超過3，即採樣圖的大小最小爲原圖的1/9。

     二、爲速度和效果綜合考慮，能夠採用下采樣係數爲2，這是雙線程插值實際上是求四個相鄰像素的平均值，所以能夠有較大的優化空間。

     一樣的640*480的圖像，使用2*2下采樣時約爲40ms(均值模糊）以及55ms(高斯模糊);

     在Real-Time O(1) Bilateral Filtering一文中有一下幾段話：

     As visible, our results are visually very similar to the exact even using very small number of PBFICs. To achieve acceptable PSNR value ( dB) for variance2
R ∈ , our method generally requires to PBFICs, and the running time is about 3.7ms to 15ms for MB image.

     For a typical 1MB image, Porikli’s method runs at about second. Our GPU implementation runs at about frames per second using 8 PBFICs (Computation complexity of Recursive Gaussian filtering is about twice the box filtering)......

     我對此速度表示嚴重懷疑，第一論文中說道他的算法佔用內存數是大概4倍圖像大小，那基本上就是採用上面代碼相似的流程，這個流程有個嚴重的後果就是兩個取樣點的計算必須按大小的順序進行，那這個並行就是個難題。另外，咱們知道，8個PBFICs的過程就包括16個均值模糊或高斯模糊的過程（1MB大小的圖像，就是1024*1024大小的灰度圖），就憑這個過程在3.5或者15ms能執行完畢的機器或許還不多見吧。GPU有着能耐？抑或是做者使用的是超級計算機，不知道各位大神贊成嗎？

    所以，論文的標題 Real - Time 是否是值得商榷呢？

    相關工程參考：http://files.cnblogs.com/Imageshop/FastBilateralFilterTest.rar