双边滤波和冲击滤波

http://blog.chinaaet.com/helimin/p/5100018184

http://blog.csdn.net/jfuck/article/details/8932978

function B = bfilter2(A,w,sigma)  
%A为给定图像，归一化到[0,1]的double矩阵  
%W为双边滤波器（核）的边长/2  
%定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2) 
%    This function implements 2-D bilateral filtering using
%    the method outlined in:
%
%       C. Tomasi and R. Manduchi. Bilateral Filtering for 
%       Gray and Color Images. In Proceedings of the IEEE 
%       International Conference on Computer Vision, 1998. 
%
%    B = bfilter2(A,W,SIGMA) performs 2-D bilateral filtering
%    for the grayscale or color image A. A should be a double
%    precision matrix of size NxMx1 or NxMx3 (i.e., grayscale
%    or color images, respectively) with normalized values in
%    the closed interval [0,1]. The half-size of the Gaussian
%    bilateral filter window is defined by W. The standard
%    deviations of the bilateral filter are given by SIGMA,
%    where the spatial-domain standard deviation is given by
%    SIGMA(1) and the intensity-domain standard deviation is
%    given by SIGMA(2).
%
% Douglas R. Lanman, Brown University, September 2006.
% dlanman@brown.edu, http://mesh.brown.edu/dlanman

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pre-process input and select appropriate filter.

% 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);    %大于w的最小整数

% Verify bilateral filter standard deviations.
if ~exist('sigma','var') || isempty(sigma) || ...
      numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0
   sigma = [3 0.1];
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.计算过程比较慢，创建waitbar可实时看到进度
h = waitbar(0,'Applying bilateral filter...');
set(h,'Name','Bilateral Filter Progress');

% Apply bilateral filter.
%计算值域核H 并与定义域核G 乘积得到双边权重函数F
dim = size(A);      %得到输入图像的width和height
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.
close(h);

clear all;
Image_pri = imread('1.png');
Image_normalized = im2double(Image_pri);
w = 5;      %窗口大小
sigma = [3 0.1];    %方差
Image_bf = bfilter2(Image_normalized,w,sigma);
Image_bfOut = uint8(Image_bf*255);
figure(1);
subplot(1,2,1);
imshow(Image_pri);
subplot(1,2,2);
imshow(Image_bfOut);

filter_gaussian = fspecial('gaussian',[5,5],3);   %生成gaussian空间波波器
gaussian_image = imfilter(Image_pri,filter_gaussian,'replicate');
figure(2);
subplot(1,2,1);
imshow(Image_pri);
subplot(1,2,2);
imshow(gaussian_image);

 双边滤波与高斯滤波器相比，对于图像的边缘信息能过更好的保存。其原理为一个与空间距离相关的高斯函数与一个灰度距离相关的高斯函数相乘。

         空间距离：指的是当前点与中心点的欧式距离。空间域高斯函数其数学形式为：

其中（xi,yi）为当前点位置，（xc,yc）为中心点的位置，sigma为空间域标准差。

         灰度距离：指的是当前点灰度与中心点灰度的差的绝对值。值域高斯函数其数学形式为：

其中gray(xi,yi)为当前点灰度值，gray(xc,yc)为中心点灰度值，sigma为值域标准差。

         对于高斯滤波，仅用空间距离的权值系数核与图像卷积后，确定中心点的灰度值。即认为离中心点越近的点，其权重系数越大。

         双边滤波中加入了对灰度信息的权重，即在邻域内，灰度值越接近中心点灰度值的点的权重更大，灰度值相差大的点权重越小。此权重大小，则由值域高斯函数确定。

         两者权重系数相乘，得到最终的卷积模板。由于双边滤波需要每个中心点邻域的灰度信息来确定其系数，所以其速度与比一般的滤波慢很多，而且计算量增长速度为核大小的平方。

空间域sigma选取：其中核大小通常为sigma的6*sigma + 1。因为离中心点3*sigma大小之外的系数与中点的系数只比非常小，可以认为此之外的点与中心点没有任何联系，及权重系数为0.OpenCV中默认的计算公式也是如此，OpenCV参考文档内容如下：“对应高斯参数的 Gaussian sigma (标准差). 如果为零，则标准差由下面的核尺寸计算： sigma = (n/2 - 1)*0.3 + 0.8, 其中 n=param1 对应水平核,n=param2对应垂直核.”

         值域sigma选取：另灰度差△g  =  abs（gray(xi,yi)- gray(xc,yc)），忽略常数的影响，因此其函数可以简化为：

由的图像可知：

已知 0≤△g≤255；

1）假设sigma = 255，当△g = 255时，系数为exp（-1） = 0.3679，当△g = 0时，系数为exp（-0）= 1.灰度最大点的系数与相差最小的灰度值系数之比为 0.3679.

2）假设sigma = 122.5，当△g = 255时，系数为exp（-4） = 0.0183，当△g = 0时，系数为exp（-0）= 1.灰度差最大点的系数与相差最小的灰度值系数之比为 0.0183.

结论：因为导数为，其增长速度为指数增长。

当simga较大时，灰度差最大值与最小值的系数在很小的一个范围之内，其比值较大。及灰度差较大的点，对于中心点也会有相应的较大的权值，此与双边滤波的保留边缘的初衷相违背。

当sigma较小时，灰度差最大值与最小值的系数在较大的一个范围之内，其比值很小，及灰度差较大的点，对应中心点仅有很小的权重。

综上分析可知：

Sigma越大，边缘越模糊，极限情况为simga无穷大，值域系数近似相等（忽略常数时，将近为exp（0）= 1），与高斯模板（空间域模板）相乘后可认为等效于高斯滤波。

Sigma越小，边缘越清晰，极限情况为simga无限接近0，值域系数近似相等（接近exp（-∞） =  0），与高斯模板（空间域模板）相乘后，可近似为系数皆相等，等效于源图像。