双边滤波 MATLAB代码实现

双边滤波与一般的高斯滤波的不同就是：双边滤波既利用了位置信息<or 几何信息——高斯滤波只用了位置信息>又利用了像素信息来定义滤波窗口的权重。

像素值越接近，权重越大。双边滤波会去除图像的细节信息，又能保持边界。

对于彩色图像，像素值的接近与否不能使用RGB空间值，双边滤波的原始文献建议使用CIE颜色空间。

代码如下：

   
   

     
     

      
      

     
     
     
     

      
      

       
       function resultI = BilateralFilt2(I,d,sigma)
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %%%
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %
       
       Author
       
       ：
       
       LiFeiteng
       
       
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %
       
       Version
       
       ：1.0——灰色图像  
       
       Time
       
       ：2013/05/01
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %
       
       Version
       
       ：1.1——灰色/彩色图像  
       
       Time
       
       ：2013/05/02 2013/05/05
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %
       
       d
       
        半窗口宽度
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       I
       
        = 
       
       double
       
       (I)
       
       ;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       if size(I,
       
       3)==
       
       1
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           resultI = BilateralFiltGray(I,d,sigma);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       elseif size(I,
       
       3)==
       
       3
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           resultI = BilateralFiltColor(I,d,sigma);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       else 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           error(
       
       ’Incorrect image size’)    
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       function resultI = BilateralFiltGray(I,d,sigma)
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       [
       
       m
       
        
       
       n
       
       ] = 
       
       size
       
       (I)
       
       ;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI = ReflectEdge(I,d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       resultI = zeros(m,n);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       width = 
       
       2*d+
       
       1;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %Distance
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       D = fspecial(
       
       ’gaussian’,[width,width],sigma(
       
       1));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       S = zeros(width,width);%pix Similarity
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       h = waitbar(
       
       0,
       
       ’Applying bilateral filter…’);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       set(h,
       
       ’Name’,
       
       ’Bilateral Filter Progress’);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       for i=
       
       1+d:m+d
      
      
     
     

     
     

      
      

     
     
     
     

      
      
    
       
       for j=
       
       1+d:n+d
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               pixValue = newI(i-d:i+d,j-d:j+d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               subValue = pixValue-newI(i,j);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               S = exp(-subValue.^
       
       2/(
       
       2*sigma(
       
       2)^
       
       2));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               H = S.*D;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               resultI(i-d,j-d) = sum(pixValue(:).*H(:))/sum(H(:)); 
      
      
     
     

     
     

      
      

     
     
     
     

      
      
    
       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           waitbar(i/m);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       close(h);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       function resultI = BilateralFiltColor(I,d,sigma)
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       I
       
        = 
       
       applycform
       
       (I,makecform(’srgb2lab’)
       
       );
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       [m n ~] = size(I);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI = ReflectEdge(I,d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       resultI = zeros(m,n,
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       width = 
       
       2*d+
       
       1;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %Distance
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       D = fspecial(
       
       ’gaussian’,[width,width],sigma(
       
       1));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       % [X,Y] = meshgrid(-d:d,-d:d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       % D = exp(-(X.^
       
       2+Y.^
       
       2)/(
       
       2*sigma(
       
       1)^
       
       2));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       S = zeros(width,width);%pix Similarity
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       h = waitbar(
       
       0,
       
       ’Applying bilateral filter…’);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       set(h,
       
       ’Name’,
       
       ’Bilateral Filter Progress’);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       sigma_r = 
       
       100*sigma(
       
       2);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       for i=
       
       1+d:m+d
      
      
     
     

     
     

      
      

     
     
     
     

      
      
    
       
       for j=
       
       1+d:n+d
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               pixValue = newI(i-d:i+d,j-d:j+d,
       
       1:
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               %subValue = pixValue-repmat(newI(i,j,
       
       1:
       
       3),width,width);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               dL = pixValue(:,:,
       
       1)-newI(i,j,
       
       1);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               da = pixValue(:,:,
       
       2)-newI(i,j,
       
       2);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               db = pixValue(:,:,
       
       3)-newI(i,j,
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               S = exp(-(dL.^
       
       2+da.^
       
       2+db.^
       
       2)/(
       
       2*sigma_r^
       
       2));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               H = S.*D;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               H = H./sum(H(:));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               resultI(i-d,j-d,
       
       1) = sum(sum(pixValue(:,:,
       
       1).*H)); 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               resultI(i-d,j-d,
       
       2) = sum(sum(pixValue(:,:,
       
       2).*H));    
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
               resultI(i-d,j-d,
       
       3) = sum(sum(pixValue(:,:,
       
       3).*H));    
      
      
     
     

     
     

      
      

     
     
     
     

      
      
    
       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           waitbar(i/m);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       close(h);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       resultI = applycform(resultI,makecform(
       
       ’lab2srgb’));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     
   
   

其中newI = ReflectEdge(I,d); %对称地扩展边界，在原始图像I的边界处镜像映射像素值 

   
   

     
     

      
      

     
     
     
     

      
      

       
       function newI = ReflectEdge(I,d)
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %Version：
       
       1.0——灰色图像  Time：
       
       2013/
       
       05/
       
       01
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %Version：
       
       1.1——灰色/彩色图像  Time：
       
       2013/
       
       05/
       
       02
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %考虑到实用性，决定不添加更多的边界处理选择，统一使用：reflect across edge
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       if size(I,
       
       3)==
       
       1
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           newI = ReflectEdgeGray(I,d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       elseif size(I,
       
       3)==
       
       3
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           newI = ReflectEdgeColor(I,d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       else 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
           error(
       
       'Incorrect image size')    
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       function newI = ReflectEdgeGray(I,d)
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       [m n] = size(I);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI = zeros(m+
       
       2*d,n+
       
       2*d);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %中间部分
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(d+
       
       1
       
       :d+m,d+
       
       1
       
       :d+n) = I;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %上
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       1
       
       :d,d+
       
       1
       
       :d+n) = I(
       
       d:-
       
       1
       
       :
       
       1,
       
       :);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %下
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       end-
       
       d:
       
       end,d+
       
       1
       
       :d+n) = I(
       
       end:-
       
       1
       
       :end-d,
       
       :);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %左
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       :,
       
       1
       
       :d) = newI(
       
       :,
       
       2*
       
       d:-
       
       1
       
       :d+
       
       1);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %右
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       :,m+d+
       
       1
       
       :m+
       
       2*d) = newI(
       
       :,m+
       
       d:-
       
       1
       
       :m+
       
       1);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       function newI = ReflectEdgeColor(I,d)
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %扩展图像边界
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       [m n ~] = size(I);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI = zeros(m+
       
       2*d,n+
       
       2*d,
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %中间部分
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(d+
       
       1
       
       :d+m,d+
       
       1
       
       :d+n,
       
       1
       
       :
       
       3) = I;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %上
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       1
       
       :d,d+
       
       1
       
       :d+n,
       
       1
       
       :
       
       3) = I(
       
       d:-
       
       1
       
       :
       
       1,
       
       :,
       
       1
       
       :
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %下
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       end-
       
       d:
       
       end,d+
       
       1
       
       :d+n,
       
       1
       
       :
       
       3) = I(
       
       end:-
       
       1
       
       :end-d,
       
       :,
       
       1
       
       :
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %左
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       :,
       
       1
       
       :d,
       
       1
       
       :
       
       3) = newI(
       
       :,
       
       2*
       
       d:-
       
       1
       
       :d+
       
       1,
       
       1
       
       :
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %右
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       newI(
       
       :,m+d+
       
       1
       
       :m+
       
       2*d,
       
       1
       
       :
       
       3) = newI(
       
       :,m+
       
       d:-
       
       1
       
       :m+
       
       1,
       
       1
       
       :
       
       3);
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       end
      
      
     
     
   
   

测试用例：

   
   

     
     

      
      

     
     
     
     

      
      

       
       img = imread(
       
       '.\lena.tif');
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %%img = imread(
       
       '.\images\lena_gray.tif');
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       img = 
       
       double(img)/
       
       255;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       img = img+
       
       0.05*randn(size(img));
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       img(img<
       
       0) = 
       
       0; img(img>
       
       1) = 
       
       1;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       %img = imnoise(img,
       
       'gaussian');
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       figure, imshow(img,[])
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       title(
       
       '原始图像')
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       d = 
       
       6;
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       sigma = [
       
       3 
       
       0.1];
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       resultI = BilateralFilt2(
       
       double(img), d, sigma);
      
      
     
     

     
     

      
      

     
     
     
     

      
      
 
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       figure, imshow(resultI,[])
      
      
     
     

     
     

      
      

     
     
     
     

      
      

       
       title(
       
       '双边滤波后的图像')
      
      
     
     
   
   

结果：

Reference：

1.C Tomasi, R Manduchi.Bilateral Filtering for Gray and Color Images, - Computer Vision, 1998.

            </div>