Bilateral Filter：保护边缘的平滑滤波器

转自：http://blog.163.com/yuyang_tech/blog/static/21605008320130242441886/

原理
滤波器是图像处理和计算机视觉中最基础的运算。而Bilateral Filter又是十分经典的一种滤波器，这主要得益于它的一个突出的特点，就是对图像进行平滑时，能进行边缘保护。
而Bilateral Flter的这个特性主要是因为他在平滑滤波时同时考虑了像素间的几何距离和色彩距离。
下面将详细具体讲述原因。
通俗的讲，对图像进行滤波就是一个加权平均的运算过程，滤波后图像中的每个像素点都是由其原图像中该点临域内多个像素点值的加权平均。不同的滤波器，最根本的差异就是权值不同。
Bilateral Filter具有两个权重，分别与几何距离和色彩距离相关，它的一般表达式为：

式中，权重和邻域像素点与中心像素点之间的几何距离相关，权重和邻域像素点与中心像素点之间的色彩距离相关。
也就是说，随着与之间几何距离和色彩距离的变化，像素点的权重和随之变化。这样的一个滤波器是一个智能的可自动调节权重的滤波器。至于如何自动调节，就取决于权重函数和的具体表达式了。
最简单和常见的一种Bilteral Filter是两个权重函数均为高斯函数，即：
 

其中，||x-y||求的是两个值x,y之间的欧式距离。
采用高斯函数之后，当邻域像素点与中心像素点的欧式距离比较大，或是两个像素值相差比较大时，该像素点的权重就比较小，导致该像素点对滤波后的像素点的影响较小。更进一步，每个滤波后像素点的值，受与他色彩相近并且距离较近的像素点的影响较大，而受其他临域像素点的影响较小，这样保护了图像中原像素点的色彩，滤掉了少数存在的“噪声点”，起到了平滑的作用，同时又保护了边缘。
举个例子，如下图所示：

（a）为图像中的一个patch，图中用高度代表该像素点的像素值。
（b）为求得（a）中右侧边缘上的一个像素点的权重（图中用高度代表权重的大小），可以看到，颜色越相近值越大，距离越近值越大，那么最终两个权重的积就如图中分布。
（c）对（a）图中的每一个像素点都进行Bilateral滤波，结果如图（c）所示，可以看到滤波的结果即保护了边缘，同时平滑了图像。
实验&Code
对于灰度图像进行bilateral滤波，matlab代码如下：

function out_img = bilateral_filter(in_img,w,sigma_c,sigma_s);
% matlab code for bilateral filter
% Pre-compute Gaussian distance weights.
[X,Y] = meshgrid(-w:w,-w:w);
C = exp(-(X.^2+Y.^2)/(2*sigma_c^2));
% Apply bilateral filter.
[m,n,c] = size(in_img);
out_img = zeros(m,n,c);
for i = 1:m
   for j = 1:n
       % Extract local region.
       I = in_img(max(i-w,1):min(i+w,m),max(j-w,1):min(j+w,n));
       % Compute Gaussian intensity weights.
       S = exp(-(I-in_img(i,j)).^2/(2*sigma_s^2));
       % Calculate bilateral filter response.
       F = S.*C((max(i-w,1):min(i+w,m))-i+w+1,(max(j-w,1):min(j+w,n))-j+w+1);
       out_img(i,j) = sum(F(:).*I(:))/sum(F(:));         
   end
end

对于彩色图像进行bilateral滤波，matlab代码如下：

function out_img = bilateral_filter(in_img,w,sigma_c,sigma_s);
% matlab code for bilateral filter]
% Convert input sRGB image to CIELab color space.
% in_img = colorspace('Lab<-RGB',in_img);
% Pre-compute Gaussian distance weights.
[X,Y] = meshgrid(-w:w,-w:w);
C = exp(-(X.^2+Y.^2)/(2*sigma_c^2));
% Rescale range variance (using maximum luminance).
sigma_s = 100*sigma_s;
% Apply bilateral filter.
[m,n,c] = size(in_img);
out_img = zeros(m,n,c);
for i = 1:m
   for j = 1:n
       % Extract local region.
       I = in_img(max(i-w,1):min(i+w,m),max(j-w,1):min(j+w,n),:);
       % Compute Gaussian intensity weights.
       dL = I(:,:,1)-in_img(i,j,1);
       da = I(:,:,2)-in_img(i,j,2);
       db = I(:,:,3)-in_img(i,j,3);
       S = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_s^2));
       
       % Calculate bilateral filter response.
       F = S.*C((max(i-w,1):min(i+w,m))-i+w+1,(max(j-w,1):min(j+w,n))-j+w+1);  
       out_img(i,j,1) = sum(sum(F.*I(:,:,1)))/sum(F(:));
       out_img(i,j,2) = sum(sum(F.*I(:,:,2)))/sum(F(:));
       out_img(i,j,3) = sum(sum(F.*I(:,:,3)))/sum(F(:));
   end
end
% Convert filtered image back to sRGB color space.
% out_img = colorspace('RGB<-Lab',out_img); 
两次运行的主程序：
clear all;
close all;
in_img = imread('lena.jpg');% or in_img = rgb2gray(imread('lena.jpg'));
in_img = rgb2gray(imread('1.jpg'));
in_img = double(in_img)/255;
in_img = in_img+0.05*randn(size(in_img));
in_img(in_img<0) = 0; in_img(in_img>1) = 1;
% Set bilateral filter parameters.
w     = 10;       % bilateral filter half-width
sigma_c = 3;% bilateral filter standard deviations
sigma_s = 0.1;
out_img = bilateral_filter(in_img,w,sigma_c,sigma_s);
figure,
subplot(1,2,1),imshow(in_img,[]),title('input image');
subplot(1,2,2),imshow(out_img,[]),title('output image'); 
实验结果：

参考文献：

C. Tomasi, R. Manduchi. “Bilateral Filtering for Gray and Color Images“, ICCV 1998.