双边滤波模板主要有两个模板生成，第一个是高斯模板，第二个是以灰度级的差值作为函数系数生成的模板。然后这两个模板点乘就得到了最终的双边滤波模板。
第一个模板是全局模板，所以只需要生成一次。第二个模板需要对每个像素都计算一次，所以需要放到循环的里面来生成，这很像表面模糊啊。哦，表面模糊就是用了一个截尾滤波器。
这里的公式我参考了这里，不过她给的第二个好像不是截尾均值滤波器，而是以灰度差值为自变量的高斯滤波器。截尾均值滤波器这里有一些理论和实现，
代码如下：
1 clear all;
2 close all;
3 clc;
4
5 img=imread('lena.jpg');
6 img=mat2gray(img);
7 [m n]=size(img);
8 imshow(img);
9
10 r=10; %模板半径
11 imgn=zeros(m+2*r+1,n+2*r+1);
12 imgn(r+1:m+r,r+1:n+r)=img;
13 imgn(1:r,r+1:n+r)=img(1:r,1:n); %扩展上边界
14 imgn(1:m+r,n+r+1:n+2*r+1)=imgn(1:m+r,n:n+r); %扩展右边界
15 imgn(m+r+1:m+2*r+1,r+1:n+2*r+1)=imgn(m:m+r,r+1:n+2*r+1); %扩展下边界
16 imgn(1:m+2*r+1,1:r)=imgn(1:m+2*r+1,r+1:2*r); %扩展左边界
17
18 sigma_d=2;
19 sigma_r=0.1;
20 [x,y] = meshgrid(-r:r,-r:r);
21 w1=exp(-(x.^2+y.^2)/(2*sigma_d^2)); %以距离作为自变量高斯滤波器
22
23 h=waitbar(0,'wait...');
24 for i=r+1:m+r
25 for j=r+1:n+r
26 w2=exp(-(imgn(i-r:i+r,j-r:j+r)-imgn(i,j)).^2/(2*sigma_r^2)); %以周围和当前像素灰度差值作为自变量的高斯滤波器
27 w=w1.*w2;
28
29 s=imgn(i-r:i+r,j-r:j+r).*w;
30 imgn(i,j)=sum(sum(s))/sum(sum(w));
31
32 end
33 waitbar(i/m);
34 end
35 close(h)
36
37 figure;
38 imshow(mat2gray(imgn(r+1:m+r,r+1:n+r)));
效果：
原图
双边滤波后