双边滤波

双边滤波(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