要害的算法思想都来自于其它一篇随想710官方网站,       在对双边网格做总括此前

 一、引言    

  钻探两岸滤波有很长一段时间了,方今看了一篇Real-提姆e
O(1) Bilateral
Filtering的杂谈,标题很吸引人,就研读了一番,经过几天的求学,基本已了然其考虑,现将这一历程做一简短的总括。

1.前言

        近年来在看Deep Bilateral Learning for Real-Time Image
Enhancement
,是一篇用CNN完结实时图像增强的点子,可以在三哥大上高速完结HDR。文章中提到到两岸网格(Bilateral
Grid)的思维,查阅了汪洋的文献,对双方网格的论争做了一下个体的接头。

     双边滤波在图像处理领域中所有广阔的选取,比如去噪、去斯特拉斯堡克、光流估摸等等,近来,相比较盛行的Non-Local算法也足以作为是双方滤波的一种增添。自从汤姆asi
et
al等人提出该算法那一天起,怎样连忙的落到实处他,从来是人们议论和切磋的关键之一,在二零一一年及二〇一二年Kunal
N.
Chaudhury等人发布的相关杂文中,指出了依照三角函数关系的值域核算法,能立见成效而又准确的已毕急迅双边算法。本文首要对此杂谈提出的办法加以讲演。

     随想大于10MB,不可以上传至今日头条,可以在那个链接下载:http://www.cs.cityu.edu.hk/~qiyang/publications/cvpr-09-qingxiong-yang.pdf

2.三头滤波器的全速近似

       在对互相网格做统计此前,先介绍一下多头滤波器的马上近似方法(A
Fast Approximation of the Bilateral Filter using a Signal Processing
Approach
),它是彼此网格的雏形,是一个将二者滤波器伸张到高维空间拓展线性卷积的法子。而两端网格的撰稿人在这篇小说的根基上,将高维空间数据映射成3D
array,进而提出了双面网格。

    
双边滤波的边缘保持特色首假使由此在卷积的进度中结成空域函数和值域核函数来落实的,典型的核函数为高斯分布函数,如下所示:

     首先,先提交一个自己自己的下结论:那篇文章无什么新意,紧要的算法思想都来源于于其余一篇随想,Fast
Bilateral Filtering for the Display of High-Dynamic-Range
Images
,而且文中的局地实验结果自己认为存在较大的水分,不过,其中提到的算法如故相比较快的。

2.1.双面滤波器

       
首先介绍一下怎么是五头滤波器,那里引用一些客人的下结论博客:彼此滤波算法原理

       
不难地说,双边滤波器是一种高斯滤波器的伸张。传统的高斯滤波器有一个高斯核,通过对空中相邻的像素点以高斯函数位权值取均值,完毕对图像的平整。由于观念高斯滤波器只考虑了空荡荡的新闻,即使它能兑现对图像的可行平滑,但与此同时也搅乱了边缘新闻。双边滤波器的规律即在观念高斯滤波器的底子上添加了一个特性亮度差别的高斯核,即既考虑了空中的新闻,由考虑了值域的新闻。在灰度差距不大的限定,表征亮度高斯核的权值较大、接近于1,因而双方滤波退化为观念的高斯滤波器,达成对图像的坦荡;在边缘部分,即便相邻像素点的上空接近,空间距离小,但出于边缘部分,灰度差距较大,由此在值域上偏离较大,所以插手的新高斯核使得在该部分不执行平滑处理,保留图像的边缘。所以说,双边滤波是一种非线性(三个高斯核的乘积)的滤波方法,是整合图像的半空中邻近度和像素值相似度的一种折中处理,同时考虑空域新闻和灰度相似性,达到保边去噪的指标。

           
710官方网站 1                    
(1)

    
杂文中对双边模糊的优化思路大约是那般的:

2.2.连忙近似

710官方网站 2

双方网络接近例子

       
作者首先将二者滤波器的计算公式写成了齐次坐标的形式(关于齐次坐标的领悟,可以参考有关齐次坐标的接头(经典)):

710官方网站 3

两岸滤波齐次表示

       
接着,小编引入δ函数,将公式从二维空间增加到三维空间(即一个2D的空间域和一个1D的值域/亮度域):

710官方网站 4

引入δ函数

        小编引入一些函数来代表上式中的一些操作:

710官方网站 5

高斯核

710官方网站 6

δ函数

        最终,将五头滤波器的公式改写为:

710官方网站 7

双面滤波的线性形式

710官方网站 8

改进后的两端滤波进度

       
小编提议,由于采样定理,可以通过对图像下采样做卷积处理,再上采样复苏原式分辨率,达到神速的目标。

710官方网站 9

一些标明的印证

    其中:

    
对于两岸模糊,离散化后的表明式大致如下所示:

2.4.伪代码

       
下图为相互近似的伪代码。1.初阶化所有的下采样齐次值,令所有下采样权值为0。2.找到全图的微小亮度点(这一步是为了持续将亮度域的限定平移到从0早先)。3.对于高分辨率图像的任一个像素点(X,Y)和亮度I(X,Y),(a)同样总括其齐次格局;(b)统计下采样后的坐标(那里运用了取整,因而,在原式高分辨率图像中的相邻空间,会被映射到低分辨率下的同一空间坐标,但在亮度域可能不在同一坐标上);(c)更新下采样空间的坐标值,即在低分辨率空间上丰裕高分辨率的齐次值(之所以用齐次值,是因为齐次坐标相加后,取出原始数据即成功了对相应数额的加权平均:(w1i1,w1)+(w2i2,w2)=(w1i1+w2i2,w1+w2),数据对应(w1i1+w2i2)/(w1+w2))。4.在采样域做高斯卷积,得到卷积后的齐次坐标。5.上采样回高分辨率图像,对于未知的数据点,使用三线性差值。

710官方网站 10

双面近似伪代码

           
710官方网站 11              (2)

      
710官方网站 12

3.三头网格

       
有大家计算了两边滤波器的近乎方法,使用一种3D数组的花样来代表那种在三维空间的三头滤波,建议了彼此网格的申辩(Real-time
edge-aware image processing with the bilateral
grid
)。

710官方网站 13

二者网格的相互滤波器完成原理

       
有了眼前双边滤波近似的理论,再来领悟两者网格就显得简单得多了。考虑上图所示的1D输入图像。双边网格是在空间域和亮度域举办采样,划分成网格。双边的边也就是从那里来的,空间(space)和亮度(range)。离散后,每个点的坐标和亮度音讯取整到相应的网格内。每个值域内的亮度值可由加权平均获得。通过在网格内展开滤波等拍卖,再由上采样方法插值出处理后的原始图像。

       
因而,构建一个双面网格的法子是:首先,依据空域和值域,划分出网格结构,并开始化所有网格的节点为0:

710官方网站 14

发轫化节点

       
对于高分辨率图像的坐标,分别除以空间采样间隔和亮度采样间隔并取整,再对应的节点处填入其原来的灰度值。考虑到取整函数会造成原图像的某一邻域内的像素点会下采样到同一个网格中,因而将灰度值写成齐次值(I(x,y),1)。并对在该网格内的具备齐次值叠加。那样,那些网格对应的尾声灰度最可以写成那一个灰度值的平均值:(I(x1,y1)+…+I(xn,yn))/n。

710官方网站 15

网格填充

       
那样,大家就将原来的图像数据转发为了一个3D的五头网格。任何图像操作函数都足以成效在那几个两岸网格上,即一定于将该双边网格左乘一个3D图像操作矩阵(例如,对于相互滤波器而言,那些函数是一个高斯卷积核,包含了半空中的方差和亮度的方差),得到一个在低分辨率意况下拍卖的结果。最终经过上采样,并开展插值,得到最终于高分辨率图像结果。上采样的规律是拔取一个参考图,对其擅自一个空间的像素举办空域和值域的采样,那里不举行取整操作,然后找到其在网格中的地点,由于没有了取整,必要动用三线性插值的方法,来兑现未知范围的亮度值的测算,那么些历程称作“slicing”。由插值上采样后,大家取得了最终的滤波结果。

    为归一化的成效。

    
f(s)是空荡荡核函数,f(r)是值域核函数,  难以直接优化上式的缘故是 f(r)的存在。  

4.两岸网格的上采样

        提议两岸网格的撰稿人接着又举办了投机的双方网格,公布了Bilateral
Guided
Upsampling
,介绍了什么运用双边网格完成部分图像操作算子。算法的主题绪想就是将一副高分辨率的图像通过下采样转换成一个两端网格,而在两者网格中,每个cell里提供一个图像变换算子,它的原理是在上空与值域相近的区域内,相似输入图像的亮度经算子变换后也应有是形似的,因此在种种cell里的操作算子能够当做是输入/输出的近乎曲线,也即一个仿射模型,而对于不一样的cell,小编通过给定输入和梦想输出去操练这几个双方网格完结其仿射模型的全局分段平滑。再经过上采样,得到高分辨率的拍卖后的图像。

710官方网站 16

Bilateral Guided Upsampling

   
σs为空白高斯函数的标准差,σr为值域高斯函数的标准差,Ω表示卷积的定义域。
可知,在图像的平缓区域,f(y)-f(x)的值变化很小,对应的值域权重接近于1,此时空域权重起第一成效,约等于直接对此区域进行高斯模糊,在边缘区域,f(y)-f(x)会有较大的出入,此时值域全面会骤降,从而致使那里整个核函数的遍布的下落,而保持了边缘的底细音信。

    
散文中提出一种思路,如果上式中固定I(X)的值,则对此每一个区其他I(y)值,上式的成员就相当于对fR(I(x),I(y))*I(y)拓展空域核卷积运算,分母则是对fR(I(x),I(y))进行空域核卷积元算,而那种卷积运算有着飞跃的算法。
这样,大家在图像的值域范围内选定若干有代表性的I(X)值,分别开展卷积,然后对于图像中的其余的像素值,进行线性插值得到。

    
直接的编码完成上述进程是分外耗时的,其时间复杂度为O(σs2) ,因此严重的限制住了该算法的加大和实际行使。不断有大家提议了然决的主意,其中Porikli基于有些假设对此进程进展了优化,比如自己就完成过里面一种:空域函数为均值函数,值域为其它其余函数,此时得以用直方图技术举行拍卖,可削减总括量,但自己的举行评释该算法那速度仍然慢,并且职能也不好。

    
算法的紧要性进献也就在此间,而以此想法是从Fast Bilateral Filtering for
the Display of High-Dynamic-Range
Images
一文中拿走的,并且在此文中还涉及了进行subsampleing举行进一步的优化,即那几个抽样卷积可以在原图的小图中开展,然后最后的结果在原图中通过双线性插值获取。

     在2011的杂谈《法斯特 O(1) bilateral
filtering using trigonometric range kernels》中,作者提议了用Raised
cosines函数来逼近高斯值域函数,并运用部分风味把值域函数分解为一些列函数的增大,从而达成函数的加快。下边大家紧要描述下该进度。

     
关于间接采样然后插值的算法源代码可以参照:http://files.cnblogs.com/Imageshop/qx_constant_time_bilateral_filter.rar

二、推导

      下边为其首要的兑现代码:

     1、一些基础理论和常识。

 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)
Cos函数在[-Pi/2,Pi/2]时期为非负、对称、在半周期内单调递增以及且有峰值的函数;

   我那边对中间的代码举行简易的讲述:

    (2) 欧拉公式:
 exp(ix)=cos(x)+isin(x);

     
1、第13、14行是获取图像的动态范围,即所有最大亮度和纤维亮度的像素值。

    (3) 分配律:
exp(a+b)=exp(a)*exp(b);

      2、
第22行的m_nr_scale是测算在原图中的取样数。26-29行中的m_grayscale是用来记录取样点的值的,比如要是动态范围是[0,255],取样数,5,则m_grayscale的值分别为0、64、128、192、255,
即对式(1)中的I(x)先固定为这5个值,计算式(1)的结果。

    (4)
图像的动态范围:[0,T],比如对于灰度图像即为[0,255];

     
3、第32到第40行直接的那么些代码其实是为了节约内存的,因为要是取样点为5,那么就须要5*2倍原图大小内存的长空来囤积取样点的卷积值,但是只要自己按取样点的分寸顺序统计,那么每计算一个取样点后(第二个除外,这就是70行的判定),就可以把原图中夹子于这些取样点及那么些取样点以前尤其取样数量里面的像素的结果值通过两者之间的线性插值获取。那种方案就足以只须求2*2倍原图大小的内存。可是那种方案就关系到插值的相继,32到40就是拍卖那样的经过的,实际的写法你可以有好多种,上边的代码写的很烂的。

     2、一些立竿见影的论据

     
4、52到61里面的代码是看空域你是用哪些类型的卷积函数,那里能够动用任意的其余的卷积函数,而有关的卷积函数也得以时任意的,那么些能够参考代码中的color_weighted_table_update函数内的代码。

     (1) 对于式子:

      
5、第72到87行的代码就是对其飞取样点的数据开展插值的经过,注意一些边缘的处理进程。

            
  710官方网站 17                                                                      
(3)

   
用插值+Sub萨姆pleing的代码可以从这边下载:http://files.cnblogs.com/Imageshop/qx_constant_time_bilateral_filter%28%E5%A2%9E%E5%BC%BA%E7%89%88%29.rar

        
其中s是自变量,取值范围[-T,T],令γ= Pi /
2T,则γs的值在[-Pi/2,Pi/2]内。此时,可以证实:

    试验结果(SigmaS=10,SigmaR=30,使用高斯卷积核函数):

           
710官方网站 18               
(4)

710官方网站 19  710官方网站 20 
710官方网站 21710官方网站 22 
710官方网站 23 
710官方网站 24

       (2)
当N丰裕大时,有下式创制:

          原图                       
 上述算法的结果                     
标准的结果

          
710官方网站 25                (5)

  上述的抽样数是按照第22行的m_nr_scale设置的,可知,视觉上如同两者之间没有怎么差别。

    即使令ρ=γσ,则上式就成为:

     按照m_nr_scale的计算办法,假若SigmaR比较小,m_nr_scale值也会比较大,对于一些工程应用,往往SigmaR就是要取相比较小的值才能保险住边缘。因而,我们品尝修改m_nr_scale的值,实际的测试注明,将m_nr_scale的值再该小一半,也能获得很为卓越的功用,而速度确可以拉长一倍。

           710官方网站 26       (6)

    
别的,上述代码还应对m_nr_scale的矮小值做个限制,m_nr_scale必须至少当先等于2的,否则不可以插值的。

      
同样,上面创造的尺码也必须有:

     在进度上,使用那种格局足够有的别样的优化技术,SigmaR=30(SigmaS对速度没有影响)时,一副640*480的彩色图像,在I3的记录簿上耗时约为75ms(值域使用均值模糊)、125ms(值域使用高斯函数)。

          
710官方网站 27

    
随笔中增强的下采样技术拓展速度的升官,我的见地看状态接纳。我自己也开展了编程,得出的定论是:

       当γs的值在[-Pi/2,Pi/2]时,因而只须求710官方网站 28
即可,此时需要
710官方网站 29

    
1、下采样的全面越小,结果和准确值偏差越大,并且此时因为下采样造成高斯滤波或者均值滤波的增速已经在总体耗时里占有的比重不大了,此时主要的龃龉是最后的双线性插值以及线性插值了,由此,总体时间上无显明提高。因此,我提出采样倍数不要跨越3,即采样图的深浅最小为原图的1/9。

         式6中,最右面部分即为高斯函数,此时验证,可以用 Raised
cosines函数来就如的模仿高斯函数,我们用一段matlab函数来证实该结果:

    
2、为速度和机能综合考虑,可以选用下采样全面为2,这是双线程插值其实是求三个相邻像素的平均值,因而可以有较大的优化空间。

 

     同样的640*480的图像,使用2*2下采样时约为40ms(均值模糊)以及55ms(高斯模糊);

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

    
在Real-提姆e O(1) Bilateral Filtering一文中有眨眼间间几段话:

    
710官方网站 30

     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)……

     
从左图的曲线分布可知,N=2(肉色),N=3(红色)两条曲线一览无遗不合乎大家的定义域的要 求,分别出现了非单调递增和负值得场地。之后乘机N的持续增大,曲线越来越接近高斯分布曲线(青色曲线)。 那从实际的角度表达了公式6的不易。

     我对此速度表示沉痛困惑,第一舆论中协商他的算法占用内存数是几乎4倍图像大小,那大多就是使用地点代码类似的流水线,那个流程有个严重的结果就是四个取样点的猜想必须按大小的顺序举办,那那么些互动就是个难题。其它,大家知晓,8个PBFICs的进度就包蕴16个均值模糊或高斯模糊的经过(1MB大大小小的图像,就是1024*1024大小的灰度图),就凭这几个进度在3.5如故15ms能实施完结的机器或许还很少见吗。GPU有着能耐?抑或是小编利用的是极品统计机,不精晓诸位大神同意吗?

  3、推导

    由此,随笔的题目 Real – Time
是或不是值得商榷呢?

     
将公式(4)带入公式(6)中,获得:

   
相关工程参考:http://files.cnblogs.com/Imageshop/FastBilateralFilterTest.rar

    
     710官方网站 31      (7)

     710官方网站 32

      
将公式(7)带入原始的两头滤波的公式(1)中,有:

 

   710官方网站 33 (8)

 

其中:

 

            
710官方网站 34

 

    
式(8)的第三行使用了眼前基础理论中的第三条:分配率。注意式中的积分是对y积分,因而可以把f(x)相关部分提取到积分符号的外面。

   

     同样,对于η,可以推导得到:

 

       710官方网站 35        
(9)

 

    
注意(8)和(9)两式中后的后半有些,比如710官方网站 36 那些,可以看到那一个实在就是对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)
去除噪音

   
710官方网站 37  
710官方网站 38

                        
 lena图增添标准差为20的高斯噪音                                  
 使用σs=20、σr=40三头滤波后的结果

    
(2)普通图像的边缘保持结果

  
710官方网站 39 
710官方网站 40

            
带有噪音的玉女图1                                 使用σs=10、σr=35双面滤波后的结果   

 

710官方网站 41 710官方网站 42 710官方网站 43

  
  带有噪音的美丽的女子图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正如小时,大家观望其分布曲线,如下图:

  710官方网站 44

      
              σ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左右,由此那几个算法说实在的功效如故要命。

     
那两篇作品分别是二零一一年及二〇一二年登载的,应该是代表着脚下可比升高的技巧,我在网上日常看看有人说双边滤波不过实时,实在是不晓得那么些高人用的是何等理论,抑或是何许一级机器。

     
同自己里面的一篇博文中双指数边缘平滑滤波器用于磨皮算法的尝尝 关联的Beeps边缘保留算法比较,那里的快慢就要慢很多了,而两者的效能比较几乎大约,所以其实很纠结。

      希望有更好的法门用于该算法。

     
我用C编写了一个标准DLL,其C#调用示例可知:
http://files.cnblogs.com/Imageshop/BilateralFilterTest.rar

  710官方网站 45

 

 710官方网站 46

 

*********************************作者:
laviewpbt   时间: 2013.11.4   联系QQ:  33184777
 转发请保留本行新闻************************

 

相关文章