研究雙邊濾波有很長一段時間了，最近看了一篇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)值，上式的分子就相當於對fR(I(x),I(y))*I(y)進行空域核卷積運算，分母則是對fR(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、第13、14行是獲取圖像的動態范圍，即具有最大亮度和最小亮度的像素值。
2、 第22行的m_nr_scale是計算在原圖中的取樣數。26-29行中的m_grayscale是用來記錄取樣點的值的，比如如果動態范圍是[0,255],取樣數,5，則m_grayscale的值分別為0、64、128、192、255， 即對式（1）中的I(x)先固定為這5個值，計算式（1）的結果。
3、第32到第40行直接的這些代碼其實是為了節省內存的，因為如果取樣點為5，那么就需要5*2倍原圖大小內存的空間來存儲取樣點的卷積值，但是如果我按取樣點的大小順序計算，那么每計算一個取樣點后（第一個除外，這就是70行的判斷），就可以把原圖中夾子於這個取樣點及這個取樣點之前那個取樣數據之間的像素的結果值通過兩者之間的線性插值獲取。這種方案就可以只需要2*2倍原圖大小的內存。但是這種方案就涉及到插值的順序，32到40就是處理這樣的過程的，實際的寫法你可以有很多種，上面的代碼寫的很爛的。
4、52到61之間的代碼是看空域你是用什么類型的卷積函數，這里可以使用任意的其他的卷積函數，而至於的卷積函數也可以時任意的，這個可以參考代碼中的color_weighted_table_update函數內的代碼。
5、第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(值域使用高斯函數)。
論文中提高的下采樣技術進行速度的提升，我的看法看情況取舍。我自己也進行了編程，得出的結論是：
1、下采樣的系數越小，結果和准確值偏差越大，並且此時因為下采樣造成高斯濾波或者均值濾波的加快已經在整個耗時里占用的比例不大了，此時主要的矛盾是最后的雙線性插值以及線性插值了，因此，總體時間上無明顯提升。因此，我建議采樣倍數不要超過3，即采樣圖的大小最小為原圖的1/9。
2、為速度和效果綜合考慮，可以采用下采樣系數為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
