之前做過圖像細節增強方面的工作，處理的是紅外灰階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.
% [email protected], 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/article/details/7170471
- Bilateral Filtering for Gray and Color Images
- Fast Bilateral Filtering for the Display of High-Dynamic-Range Images