双边滤波器原理及其matlab实现

    之前做过图像细节增强方面的工作，处理的是红外灰度14bit图像，图像信号由14bit AD量化后，再经FPGA处理得到，使用非锐化掩模的方法，先用双边滤波器（BF）对原图像进行滤波得到低频部分，原图和低频作差后得到高频分量，高频分量和低频分量分别增强后再进行合成。



    双边滤波的特点是保边去噪，相较于高斯滤波，在平滑图像的同时，增加了对图像边缘的保护，其主要原因是由于该滤波器由两部分组成，一部分与像素空间距离相关，另一部分与像素点的像素差值相关。

    下面结合公式来说说为什么双边滤波在模糊图像的时候具有保边功能，双边滤波器公式为：



    其中，空间邻近度因子为



    亮度相似度因子为



双边滤波器的权重等于空间邻近度因子和亮度相似度因子的乘积



    空间邻近度因子为高斯滤波器系数，像素距离越远，权重越小，当

越大时，平滑效果越明显。

    亮度相似度因子与空间像素差值相关，像素差值越大，权重越小，这也是为什么双边滤波器能够保边去噪的原因。当

越大时，对同等灰度差的像素平滑作用越大，保边效果越差，论文中给出的参考是

一般取高斯噪声标准差的2倍。

 

下面列出matlab代码，代码下载地址，需要有Mathworks账号：

https://cn.mathworks.com/matlabcentral/fileexchange/12191-bilateral-filtering

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('academy.jpg');
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);   %生成sobel空间波波器
gaussian_image = imfilter(Image_pri,filter_gaussian,'replicate');
figure(2);
subplot(1,2,1);
imshow(Image_pri);
subplot(1,2,2);
imshow(gaussian_image);

对比双边滤波和高斯滤波结果：

输入图像为512×512的Lena灰度图像



双边滤波



高斯滤波

 

输入图像为512×512的Lena彩色图像



双边滤波



高斯滤波

注意看双边滤波图像，帽子边缘，模糊的效果很明显，但皮肤和背景处，仍然比较清晰，说明起到了保边平滑的作用。

 

下面给出一些测试结果

输入图像为256×256 Einstein灰度图像



双边滤波



高斯滤波

 

输入图像为256×256 大沸沸灰度图像



双边滤波



高斯滤波

 

输入图像为400×300的彩色图像



双边滤波



高斯滤波

 

可以看到，双边滤波相对于高斯滤波来说，主要的优势在于保边去噪，在平滑滤波的同时，保留了边缘信息。

此处代码运算速度很慢，2007年，麻省理工学院学者提出了快速双边滤波算法，并且给出了matlab代码，下载地址：

http://people.csail.mit.edu/jiawen/#code

 

参考资料：

http://blog.csdn.net/abcjennifer/article/details/7616663
http://blog.csdn.net/bugrunner/ 95c0
article/details/7170471

Bilateral Filtering for Gray and Color Images

Fast Bilateral Filtering for the Display of High-Dynamic-Range Images