双边滤波器—— Matlab实现

例:先用双边滤波器(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);   %生成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);  

 

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

双边滤波

高斯滤波

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

双边滤波

高斯滤波

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

下面给出一些测试结果

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

双边滤波,滤波后黑色领带与白色衬衫区域轮廓清晰。

高斯滤波,滤波后图像整体模糊。

 

输入图像为256×256 大沸沸彩色图像

双边滤波

高斯滤波

 

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

双边滤波

高斯滤波

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

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

 

参考资料:

1. http://blog.csdn.net/abcjennifer/article/details/7616663

2. http://blog.csdn.net/bugrunner/article/details/7170471

3. Bilateral Filtering for Gray and Color Images

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

  • 10
    点赞
  • 117
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
以下是一个基于局部自适应滤波双边滤波Matlab代码示例: ```matlab % 读取图像 image = imread('image.jpg'); % 双边滤波参数设置 sigma_spatial = 10; % 空间域标准差 sigma_range = 20; % 强度域标准差 % 双边滤波 filtered_image = bilateralFilter(image, sigma_spatial, sigma_range); % 显示原始图像和滤波后的图像 figure; subplot(1, 2, 1); imshow(image); title('原始图像'); subplot(1, 2, 2); imshow(filtered_image); title('滤波后的图像'); % 双边滤波函数 function filtered_image = bilateralFilter(image, sigma_spatial, sigma_range) [height, width, ~] = size(image); filtered_image = zeros(size(image)); % 对每个像素进行滤波 for i = 1:height for j = 1:width % 获取当前像素的空间域范围 spatial_range = getSpatialRange(image, i, j, sigma_spatial); % 获取当前像素的强度域范围 range_range = getRangeRange(image, i, j, sigma_range); % 计算加权平均值 filtered_image(i, j, :) = weightedAverage(image, i, j, spatial_range, range_range); end end end % 获取当前像素的空间域范围 function spatial_range = getSpatialRange(image, i, j, sigma_spatial) [height, width, ~] = size(image); spatial_range = zeros(height, width); % 计算每个像素与当前像素的空间距离 for x = 1:height for y = 1:width spatial_range(x, y) = sqrt((x - i)^2 + (y - j)^2); end end % 计算空间域权重 spatial_weight = exp(-(spatial_range.^2) / (2 * sigma_spatial^2)); % 归一化空间域权重 spatial_weight = spatial_weight / sum(spatial_weight(:)); % 返回空间域范围 spatial_range = spatial_weight; end % 获取当前像素的强度域范围 function range_range = getRangeRange(image, i, j, sigma_range) [height, width, ~] = size(image); range_range = zeros(height, width); % 计算每个像素与当前像素的强度差异 for x = 1:height for y = 1:width range_range(x, y) = norm(double(image(x, y, :)) - double(image(i, j, :))); end end % 计算强度域权重 range_weight = exp(-(range_range.^2) / (2 * sigma_range^2)); % 归一化强度域权重 range_weight = range_weight / sum(range_weight(:)); % 返回强度域范围 range_range = range_weight; end % 计算加权平均值 function weighted_average = weightedAverage(image, i, j, spatial_range, range_range) [height, width, channels] = size(image); weighted_average = zeros(1, 1, channels); % 对每个通道进行加权平均 for c = 1:channels weighted_average(c) = sum(sum(double(image(:, :, c)) .* spatial_range .* range_range)); end % 返回加权平均值 weighted_average = uint8(weighted_average); end ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值