双边滤波器

   双边滤波器是一种具有保边去噪特性的非线性滤波器，它比一般的滤波器多了一个高斯方差，它是基于图像空间分布的高斯滤波函数，同时它还有一个基于图像像素差的高斯滤波函数，所以该滤波器不仅与图像灰度像素值有关，而且像素间的距离也会对滤波器的作用产生影响。

 

双边滤波器的公式如下： 

 

                                                                                                                               （3-5）

 

其中I为原图像，J为经双边滤波后的图像，p、q为图像中像素点的坐标，f、g分别为空间域影响函数和灰度域影响函数，一般使用高斯分布来描述，为归一化项，为设定的像素空间，|*|和||*||分别表示*的一阶范数和二阶范数。

 

其数学表达式为：

 

 

 

                   （3-6）         

 

其中d(p,q)和δ(Ip,Iq）分别为图像两像素点之间的欧几里得距离及像素的灰度差，和分别是空间域高斯函数的方差和灰度域高斯函数的方差，它们决定了双边滤波器的性能，它们通过指代像素之间的相对距离和亮度的变化范围来决定中心像素的数值。

 

由公式（3-2）知，双边滤波器对图像像素点进行了求平均，因此能够平滑图像，另外，由于f和g均为高斯分布，则它们的乘积会随着像素之间距离和像素值的变化迅速减小，即在一个像素空间范围内，对于与中间像素的亮度值相似的像素点的影响大，而与其亮度值差异大的像素点的影响小，同时只要σd的变化范围小于它的幅度增加空间σr，就对边缘没有影响，因此该滤波器可以保持图像的边缘。

 

图像的细节和边缘主要受σr的影响，随着参数σr的增加，双边滤波器逐渐接近于高斯模糊，因为这一系列的高斯函数比较平坦，在图像覆盖的强度间隔上几乎是个常量，而增大σd可以光滑更多的特征。

 

论文中使用的双边滤波器的大小固定为5*5，空间域标准差σd固定为2.0，灰度域标准差σr初始值设为0.5，每次迭代后将其乘以0.9，以便逐渐减弱双边滤波器的影响。

代码如下：

% Pre-process input and select appropriate filter.  
function B = bfilter2(A,w,sigma)   
% Verify that the input image exists and is valid.  
if ~exist('A','var') || isempty(A)  
   error('Input image A is undefined or invalid.');  
end  
if ~isfloat(A) || ~sum([1,3] == size(A,3)) || ...  
      min(A(:)) < 0 || max(A(:)) > 1  
   error(['Input image A must be a double precision ',...  
          'matrix of size NxMx1 or NxMx3 on the closed ',...  
          'interval [0,1].']);        
end  
  
% Verify bilateral filter window size.  
if ~exist('w','var') || isempty(w) || ...  
      numel(w) ~= 1 || w < 1  
   w = 5;  
end  
w = ceil(w);  
  
% Verify bilateral filter standard deviations.  
if ~exist('sigma','var') || isempty(sigma) || ...  
      numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0  
   sigma = [2 0.9];  
end  
  
% Apply either grayscale or color bilateral filtering.  
if size(A,3) == 1  
   B = bfltGray(A,w,sigma(1),sigma(2));  
else  
   B = bfltColor(A,w,sigma(1),sigma(2));  
end  
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% Implements bilateral filtering for grayscale images.  
function B = bfltGray(A,w,sigma_d,sigma_r)  
  
% Pre-compute Gaussian distance weights.  
[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));  
  
% Create waitbar.  
h = waitbar(0,'Applying bilateral filter...');  
set(h,'Name','Bilateral Filter Progress');  
  
% Apply bilateral filter.  
%计算值域核H 并与定义域核G 乘积得到双边权重函数F  
dim = size(A);  
B = zeros(dim);  
for i = 1:dim(1)  
   for j = 1:dim(2)  
        
         % Extract local region.  
         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  
        
         % Compute Gaussian intensity weights.  
         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);  
         B(i,j) = sum(F(:).*I(:))/sum(F(:));  
                 
   end  
   waitbar(i/dim(1));  
end  
  
% Close waitbar.  
close(h);  
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% Implements bilateral filter for color images.  
function B = bfltColor(A,w,sigma_d,sigma_r)  
  
% Convert input sRGB image to CIELab color space.  
if exist('applycform','file')  
   A = applycform(A,makecform('srgb2lab'));  
else  
   A = colorspace('Lab<-RGB',A);  
end  
  
% Pre-compute Gaussian domain weights.  
[X,Y] = meshgrid(-w:w,-w:w);  
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));  
  
% Rescale range variance (using maximum luminance).  
sigma_r = 100*sigma_r;  
  
% Create waitbar.  
h = waitbar(0,'Applying bilateral filter...');  
set(h,'Name','Bilateral Filter Progress');  
  
% Apply bilateral filter.  
dim = size(A);  
B = zeros(dim);  
for i = 1:dim(1)  
   for j = 1:dim(2)  
        
         % Extract local region.  
         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,:);  
        
         % Compute Gaussian range weights.  
         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));  
        
         % Calculate bilateral filter response.  
         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  
   waitbar(i/dim(1));  
end  
  
% Convert filtered image back to sRGB color space.  
if exist('applycform','file')  
   B = applycform(B,makecform('lab2srgb'));  
else    
   B = colorspace('RGB<-Lab',B);  
end  
  
% Close waitbar.