图像处理（六）递归双边滤波磨皮

递归双边滤波是双边滤波的一种加速算法，加速比非常大，就像高斯模糊快速算法一样，加速起来，速度跟原算法相比，提高了十几倍。双边滤波的加速算法很多，文献都比较老旧，我这边主要讲一下比较新的算法：递归双边滤波，对应的paper为：《Recursive Bilateral Filtering》

这个算法比我另外一篇博文将的双指数滤波快一点，两篇文献的基本思想一样：

递归公式跟双指数的递归公式基本上一样。于是写代码就轻轻松松了，之前写过双指数的算法，接着只要在双指数的基础上，把代码改一改就OK了，这里就不详细讲解算法流程了，有兴趣的同学可以到paper的作者主页上，里面有代码，不过他的代码好像是要实现立体匹配的，把里面的相关的代码扣出来就可以了，下面是部分重要函数的代码。

[cpp] view plain copy

1. void CBRF::recursive_bilateral_filter(double *data,unsigned char**disparity,double***in,unsigned char***texture,double sigma_spatial,double sigma_range,int height,int width,int nr_channel,double***temp,double***temp_2w)  
2. {  
3.     double yp[100];  
4.     double alpha=exp(-sqrt(2.0)/(sigma_spatial*width));//filter kernel size  
5.     double inv_alpha=(1-alpha);  
6.     double range_table[256];  
7.     for(int i=0;i<=255;i++)   
8.         range_table[i]=exp(-double(i)/(sigma_range*255));  
9.   
10.   
11.     BYTE* p0 =m_pImage;  
12.     const int nChannel = 4;  
13.     int m_length =m_nWidth*m_nHeight;  
14. //水平迭代滤波  
15.     //double***in_=in;/*horizontal filtering*/  
16.     //double***out_=temp;  
17.     double ***out=new double**[m_nHeight];  
18.     for (int i=0;i<m_nHeight;i++)  
19.     {  
20.         out[i]=new double*[m_nWidth];  
21.     }  
22.     for (int i=0;i<m_nHeight;i++)  
23.     {  
24.         for (int j=0;j<m_nWidth;j++)  
25.         {  
26.             out[i][j]=new double [3];  
27.         }  
28.     }  
29.     double ***out=new double**[m_nHeight];  
30.   
31.     for(int i=0;i<m_nHeight;i++)  
32.     {  
33.         BYTE* p0 =m_pImage+i*m_nWidth;  
34.         memcpy(out[i][0],p0,sizeof(double)*m_nWidth*nChannel);  
35.         for(int j=1;j<width;j++)   
36.         {  
37.             double weight=range_table[euro_dist_rgb_max(out[i][j],out[i][j-1])];  
38.             double alpha_=weight*alpha;  
39.             for(int k=0;k<nChannel;k++)   
40.             {  
41.                 double ycc=inv_alpha*m_pImage[i*m_nWidth*nChannel+j*nChannel+k]+alpha_*yp[k];  
42.                 out[i][j][k]=ycc  
43.                 yp[k]=ycc;  
44.             }  
45.         }  
46.   
47.   
48.         int w1=width-1;  
49.         for(int c=0;c<nr_channel;c++)   
50.             out[y][w1][c]=0.5*(out[y][w1][c]+in[y][w1][c]);  
51.         memcpy(yp,out_[y][w1],sizeof(double)*nr_channel);  
52.         for(int x=width-2;x>=0;x--)   
53.         {  
54.             double weight=range_table[euro_dist_rgb_max(texture[y][x],texture[y][x+1])];  
55.             double alpha_=weight*alpha;  
56.             for(int c=0;c<nr_channel;c++)   
57.             {  
58.                 double ycc=inv_alpha*in_[y][x][c]+alpha_*yp[c];  
59.                 out_[y][x][c]=0.5*(out_[y][x][c]+ycc);  
60.                 yp[c]=ycc;  
61.             }  
62.         }  
63.     }  
64.   
65.   
66.   
67.   
68. //垂直迭代滤波  
69.     /*in_=temp; 
70.     alpha=exp(-sqrt(2.0)/(sigma_spatial*height));//filter kernel size 
71.     inv_alpha=(1-alpha); 
72.     double**ycy,**ypy,**xcy,**xpy; 
73.     unsigned char**tcy,**tpy; 
74.     memcpy(out[0][0],in_[0][0],sizeof(double)*width*nr_channel); 
75.     for(int y=1;y<height;y++) 
76.     { 
77.         tpy=texture[y-1]; 
78.         tcy=texture[y]; 
79.         xcy=in_[y]; 
80.         ypy=out[y-1]; 
81.         ycy=out[y]; 
82.         for(int x=0;x<width;x++) 
83.         { 
84.             double weight=range_table[euro_dist_rgb_max(tcy[x],tpy[x])]; 
85.             double alpha_=weight*alpha; 
86.             for(int c=0;c<nr_channel;c++)  
87.                 ycy[x][c]=inv_alpha*xcy[x][c]+alpha_*ypy[x][c]; 
88.         } 
89.     } 
90.     int h1=height-1; 
91.     ycy=temp_2w[0]; 
92.     ypy=temp_2w[1]; 
93.     memcpy(ypy[0],in_[h1][0],sizeof(double)*width*nr_channel); 
94.     for(int x=0;x<width;x++)  
95.     { 
96.         unsigned char disp=0;  
97.         double min_cost=0.5*(out[h1][x][0]+ypy[x][0]); 
98.         for(int c=1;c<nr_channel;c++) 
99.         { 
100.             double cost=0.5*(out[h1][x][c]+ypy[x][c]); 
101.             if(cost<min_cost) 
102.             { 
103.                 min_cost=cost; 
104.                 disp=c; 
105.             } 
106.         } 
107.         disparity[h1][x]=disp; 
108.     } 
109.     for(int y=h1-1;y>=0;y--) 
110.     { 
111.         tpy=texture[y+1]; 
112.         tcy=texture[y]; 
113.         xcy=in_[y]; 
114.         for(int x=0;x<width;x++) 
115.         { 
116.             double weight=range_table[euro_dist_rgb_max(tcy[x],tpy[x])]; 
117.             double alpha_=weight*alpha; 
118.             unsigned char disp=0; 
119.             double min_cost=255; 
120.             for(int c=0;c<nr_channel;c++)  
121.             { 
122.                 ycy[x][c]=inv_alpha*xcy[x][c]+alpha_*ypy[x][c]; 
123.                 double cost=0.5*(out[y][x][c]+ycy[x][c]); 
124.                 if(cost<min_cost) 
125.                 { 
126.                     min_cost=cost; 
127.                     disp=c; 
128.                 } 
129.             } 
130.             disparity[y][x]=disp; 
131.         } 
132.         memcpy(ypy[0],ycy[0],sizeof(double)*width*nr_channel); 
133.     }*/  
134. }  
135. unsigned char CBRF::euro_dist_rgb_max(unsigned char *a,unsigned char *b)  
136. {  
137.     unsigned char x,y,z;  
138.     x=abs(a[0]-b[0]);  
139.     y=abs(a[1]-b[1]);  
140.     z=abs(a[2]-b[2]);  
141.     return(max(max(x,y),z));  
142. }  

看一些我用这个算法写的demo的测试结果：

原图

美图秀秀

递归双边滤波

最后人一定要靠自己，想要深入理解一个算法，让记忆更深，一定要自己把代码敲过一遍，既然选择了IT，就要坚持写代码。本文地址：http://blog.csdn.net/hjimce/article/details/45421207    作者：hjimce     联系qq：1393852684    更多资源请关注我的博客：http://blog.csdn.net/hjimce  原创文章，转载请保留这几行信息