双边滤波(bilateral filter)是一种非线性滤波器，该算法结合空间信息和亮度相似性对图像进行滤波处理，在平滑滤波的同时能大量保留图像的边缘和细节特征。 定义如下：
其中
为输出图像，
为输入图像，
是以像素点
为中心的邻域窗口，
为滤波核。双边滤波的滤波核由两部分乘积组成：空域核
与值域核
。两个滤波核通常都采用高斯函数形式。
其中
为空域高斯函数的标准差，
为值域高斯函数的标准差。空域滤波系数由像素间的空间距离决定，距离越小，系数越大。值域滤波系数由像素间的相似度决定，像素值越接近，系数越大。
在灰度变化平缓区域，值域滤波系数接近1，此时空域滤波起主要作用，双边滤波器退化为传统的高斯低通滤波器，对图像进行平滑操作。而在图像变化剧烈的部分（即图像边缘），像素间差异较大，值域滤波起主要作用，因而能保持边缘信息。
下图是双边滤波原来示意图：
双边滤波原理示意
对彩色图像：[2]中将其转换至Lab空间，值域核由L、a、b三分量的差值决定：
但实验结果表明它与直接分别对R、G、B通道差别并不大。
程序：
clear all;
close all;
I=double(imread('mandrill.jpg'));
I=I/255;
w = 5;
sigma=[3 0.1];
O=bfilter2(I,w,sigma);
figure;
imshow([I,O],[]);
function B = bfilter2(A,w,sigma)
if size(A,3) == 1
B = bfltGray(A,w,sigma(1),sigma(2));
else
B = bfltColor(A,w,sigma(1),sigma(2));
end
function [B,K] = bfltGray(A,w,sigma_d,sigma_r)
%%
%A为给定图像，归一化到[0,1]的矩阵
%W为双边滤波器（核）的边长/2
%定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2)
%%
%计算空域核
[X,Y] = meshgrid(-w:w,-w:w);
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));
%计算值域核H并与定义域核G乘积得到双边权重函数F
[m,n] = size(A);
B = zeros(m,n);%结果图像
K = zeros(m,n);%双边滤波器
for i = 1:m
for j = 1:n
%确定邻域边界
iMin = max(i-w,1);
iMax = min(i+w,m);
jMin = max(j-w,1);
jMax = min(j+w,n);
%定义当前核所作用的区域为(iMin:iMax,jMin:jMax)
I = A(iMin:iMax,jMin:jMax);%提取该区域的源图像值赋给I
H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));
F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);
K(i,j) = sum(F(:));
B(i,j) = sum(F(:).*I(:))/K(i,j);
end
end
function B = bfltColor(A,w,sigma_d,sigma_r)
%%
%A为给定图像，归一化到[0,1]的矩阵
%W为双边滤波器（核）的边长/2
%定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2)
%%
%颜色空间转换：RGB-->Lab
if exist('applycform','file')
A = applycform(A,makecform('srgb2lab'));
else
A = colorspace('Lab<-RGB',A);
end
%计算空域核
[X,Y] = meshgrid(-w:w,-w:w);
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));
%改变值域参数以适应Lab空间的取值范围
sigma_r = 100*sigma_r;
dim = size(A);
B = zeros(dim);
for i = 1:dim(1)
for j = 1:dim(2)
%确定邻域边界
iMin = max(i-w,1);
iMax = min(i+w,dim(1));
jMin = max(j-w,1);
jMax = min(j+w,dim(2));
I = A(iMin:iMax,jMin:jMax,:);
% 计算值域核
dL = I(:,:,1)-A(i,j,1);
da = I(:,:,2)-A(i,j,2);
db = I(:,:,3)-A(i,j,3);
H = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_r^2));
F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);%双边滤波核
norm_F = sum(F(:));
B(i,j,1) = sum(sum(F.*I(:,:,1)))/norm_F;
B(i,j,2) = sum(sum(F.*I(:,:,2)))/norm_F;
B(i,j,3) = sum(sum(F.*I(:,:,3)))/norm_F;
end
end
%转回RGB空间
if exist('applycform','file')
B = applycform(B,makecform('lab2srgb'));
else
B = colorspace('RGB<-Lab',B);
end
参考：
[1 ]http://blog.csdn.net/abcjennifer/article/details/7616663
[2] TomasiC, Manduchi R. Bilateral filtering for gray and color images[C]//ComputerVision, 1998. Sixth International Conference on. IEEE, 1998: 839-846.
[3] http://www.mathworks.com/matlabcentral/fileexchange/12191-bilateral-filtering