研究双边滤波有很长一段时间了，最近看了一篇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] ",(float)image_max,(float)image_min); 24 //printf("[sigma_range: %1.3f] ",m_sigma_range); 25 //printf("[nr_scale: %d] ",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
