阅读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