数字图像处理（三）：高斯滤波和双边滤波

clear all;clc;close all
%% 彩色to灰度
img=imread(‘Fig0631(a)(strawberries_coffee_full_color).jpg‘);
gray=rgb2gray(img);                                     %把彩色图片转化成灰度图片
figure(1),imshow(gray),title(‘彩色原始图像转灰色图像）‘); %显示原始图像
gray=imnoise(gray,‘gaussian‘,0,0.001);                  %加入均值为0，方差为0.001的高斯噪声
figure(2),imshow(gray),title(‘加入高斯噪声之后的图像‘);   %显示加入高斯噪声之后的图像

%% 用matlab系统函数进行高斯滤波
sigma=0.5;                                                %滤波器的标准值，单位为像素
hsize=[3 3];                                              %模板尺寸
gsseq=fspecial(‘gaussian‘,hsize,sigma);                   %生成高斯序列
Y1=filter2(gsseq,gray)/255;                               %用生成的高斯序列进行滤波
figure(3),imshow(Y1),title(‘用Matlab自带函数进行高斯滤波‘); %显示滤波后的图像

%% 用重新编写的程序进行高斯滤波
gray=double(gray);                %将图像转为double型
gray=fft2(gray);                  %二维傅立叶变换
gray=fftshift(gray);              %频谱居中
[m,n]=size(gray);                 %计算图像大小
d0=80;                            %D0=sigma，也就是标准差
m1=fix(m/2);                      %计算图像中心
n1=fix(n/2);                      %计算图像中心
for i=1:m
    for j=1:n
        d=sqrt((i-m1)^2+(j-n1)^2);%计算像素点到图像中心的距离
        h(i,j)=exp(-d^2/2/d0^2);  %高斯滤波器
    end
end
g=gray.*h;                        %将图像进行高斯滤波，频域上表现为为两个函数相乘
g=ifftshift(g);                   %频域圆周移位
g=ifft2(g);                       %二维傅里叶反变换
g=mat2gray(real(g));              %归一化
figure(4),imshow(g),title(‘用重新编写的程序进行高斯滤波‘);%显示滤波后的图像

%双边滤波调用示例
I=imread(‘Fig0427(a)(woman).jpg‘); %读入图片
I=double(I)/255;  %转为double型并归一化
w     = 5;        % 双边滤波器半宽，w越大平滑作用越强
sigma = [3 0.1];  % 空间距离方差σd记为SIGMA(1),像素亮度方差σr记为SIGMA(2),即空间邻近度因子和亮度相似度因子的衰减程度
I1=bfilter2(I,w,sigma);                     %双边滤波器滤波
figure(1),imshow(I),title(‘原始图像‘);       %作出原始图像
figure(2),imshow(I1),title(‘双边滤波后的图像‘)%作出双边滤波后的图像

function B = bfltGray(A,w,sigma_d,sigma_r)
% 计算距离因子权重
[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));

%创建进度条
h = waitbar(0,‘Applying bilateral filter...‘);
set(h,‘Name‘,‘Bilateral Filter Progress‘);

% 应用双边滤波
%计算值域核H 并与定义域核G 乘积得到双边权重函数F
dim = size(A);
B = zeros(dim);
for i = 1:dim(1)
   for j = 1:dim(2)
         %边界限制
         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
      
         % 计算亮度因子权重
         H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));
      
         % 计算双边滤波结果
         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(h);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 对彩色图像进行双边滤波操作
function B = bfltColor(A,w,sigma_d,sigma_r)

% 将sRGB转换为CIELab色彩空间
if exist(‘applycform‘,‘file‘)
   A = applycform(A,makecform(‘srgb2lab‘));
else
   A = colorspace(‘Lab<-RGB‘,A);
end

% 计算空间距离因子权重
[X,Y] = meshgrid(-w:w,-w:w);
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));

% 调整亮度因子权重
sigma_r = 100*sigma_r;

% 创建进度条
h = waitbar(0,‘Applying bilateral filter...‘);
set(h,‘Name‘,‘Bilateral Filter Progress‘);

% 应用双边滤波
dim = size(A);
B = zeros(dim);
for i = 1:dim(1)
   for j = 1:dim(2)
         % 边界限制
         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,:);
      
         % 计算亮度因子权重
         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));
      
         % 计算双边滤波结果
         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

% 将图像转换为sRGB色彩空间.
if exist(‘applycform‘,‘file‘)
   B = applycform(B,makecform(‘lab2srgb‘));
else  
   B = colorspace(‘RGB<-Lab‘,B);
end

% 结束进度条
close(h);

%A为给定图像，归一化到[0,1]的矩阵
%w为双边滤波器（核）的边长/2
%定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function B = bfilter2(A,w,sigma)
% 检验给定图像是否存在并且有效
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

% 检验双边滤波器的半宽是否符合要求
if ~exist(‘w‘,‘var‘) || isempty(w) || ...
      numel(w) ~= 1 || w < 1
   w = 5;
end
w = ceil(w);

% 检验sigma参数是否符合要求
if ~exist(‘sigma‘,‘var‘) || isempty(sigma) || ...
      numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0
   sigma = [3 0.1];
end

%选择彩色模式或灰度模式
if size(A,3) == 1
   B = bfltGray(A,w,sigma(1),sigma(2));
else
   B = bfltColor(A,w,sigma(1),sigma(2));
end