双边滤波

1、原理介绍

双边滤波由C. Tomasi在1998年提出，是一种经典的非线性空间滤波方法。在滤波器稀疏的制定上，双边滤波同时考虑到了输出像素与邻域内其它像素的欧氏距离和取值的差异，即：同时考虑到了空间域和值域间的差别。如维纳滤波和高斯滤波等只考虑了空间域的滤波方法，在滤波后对边缘信息的保护效果不理想；如α-截尾均值滤波器等只考虑值域的滤波方法，在滤波后图像整体模糊，不能有效的保护细节信息。双边滤波器综合考量了空间域和值域对于滤波产生的影响，因而能达到保持边缘，降噪平滑的效果，是一种良好的边缘保持滤波器。
双边滤波通过基于高斯分布的加权平均方法实现。以图像中具体的像素值求解为里说明双边滤波的实现原理：图像在(i,j)处经过双边滤波后的输出像素值g依赖于邻域内像素值f的加权组合。
$$g(i,j)=\frac{{\sum }_{k,l}f(k,l)w(i,j,k,l)}{{\sum }_{k,l}w(i,j,k,l)}$$
其中，k,l表示邻域像素的位置，权重系数w(i,j,k,l)为空间域核d与值域核r的乘积。空间域核d是指基于高斯函数计算当前点与中心点的欧式距离。
$$d(i,j,k,l)=\exp \left(-\frac{(i-k{)}^2+(j-l{)}^2}{2{\sigma }_d^2}\right)$$
值域核r是指基于高斯函数计算当前点与中心点像素值的差的绝对值。
$$r(i,j,k,l)=\exp \left(-\frac{\Vert f(i,j)-f(k,l){\Vert }^2}{2{\sigma }_r^2}\right)$$
由空间域核与值域核的计算公式可得权重系数的计算公式为：
$$w(i,j,k,l)=\exp \left(-\frac{(i-k{)}^2+(j-l{)}^2}{2{\sigma }_d^2}-\frac{\Vert f(i,j)-f(k,l){\Vert }^2}{2{\sigma }_r^2}\right)$$

2、matlab程序实现

双边滤波实现函数：

%适用于单通道图像的双边滤波程序

function B = Bilater_Gray(A,w,sigma_d,sigma_r)
         %输出参数：
         % A为待滤波图像（double类型，取值在[0，1]）
         % w为滤波窗口的半径（e.g:3*3窗口的w值为1，w=3时的滤波效果较好）
         % sigma_d为定义域（空间域）核的方差，通常设置为3
         % sigma_r为值域核的方差，通常设置为0.1
         %输出参数：
         % B为滤波后的图像
% 预先计算高斯距离权重
[X,Y] = meshgrid(-w:w,-w:w);
%创建核距离矩阵，e.g.
%  [x,y]=meshgrid(-1:1,-1:1)
% 
% x =
% 
%     -1     0     1
%     -1     0     1
%     -1     0     1
% 
% 
% y =
% 
%     -1    -1    -1
%      0     0     0
%      1     1     1
%计算定义域核
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));


%计算值域核H 并与定义域核G 乘积得到双边权重函数F
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));
         %定义当前核所作用的区域为(iMin:iMax,jMin:jMax)
         I = A(iMin:iMax,jMin:jMax);%提取该区域的源图像值赋给I

         %计算值域核H.
         H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));

         % Calculate bilateral filter response.
         F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);
         %在计算边缘部分的点的时候H的大小会变化，例如在计算第一行第一列的点时，
         %H的大小为4*4，因为7*7的其余部分都在图像外面（没有对应的值）。
         %因此适当的对G进行裁切使得G的大小始终能和H对上
         B(i,j) = sum(F(:).*I(:))/sum(F(:));

   end
end

主程序中调用：

%双边滤波主程序
clear all;
close all;
clc;
% 输入图像
fname   = 'Images\lena.jpg';   #改成你要操作的图像
X = double(rgb2gray(imread(fname)))/255;
% 开始滤波
Y = Bilater_Gray(X,3,3,0.1);

%获取细节层，即：被过滤的部分
Z = X - Y;

%结果显示
figure; imshow(X,[]);
figure; imshow(Y,[]);
figure; imshow(Z,[]);
%因为Z中有负值，所以最终强制以图片的形式显示的时候为灰色（正常现象）

结果展示：