对数字图像的线性空间域滤波、非线性空间域滤波、频率域滤波、小波滤波

学了一段时间冈萨雷斯的《数字图像处理》，结合自己查询，总结一些针对数字图像的空间域滤波与频率域滤波实现方法，代码中的解释得已经很详细了

主函数：

分别是线性空间域滤波、非线性空间域滤波、频率域滤波、小波滤波的实现，最后还计算了峰值信噪比PSNR

clc,clear,close all
f = imread('C:\Users\HS\Desktop\第一章\图像-灰度图.jpg');imshow(f);
title('原图');set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);

%% 空间域滤波

% 二维非线性空间滤波器--中值滤波（最有名的排序统计滤波器），
%调用工具箱medfilt2,默认在3*3领域上计算中值，尝试三个不同的边界填充选项
gm = medfilt2(f,[3 3],'zeros');figure,imshow(gm);hold on
title('中值滤波 zeros,PSNR=39.159');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);hold off

gms = medfilt2(f,[3 3],'symmetric');figure,imshow(gms);
title('中值滤波 symmetric,PSNR=40.038');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);

gmi = medfilt2(f,[3 3],'indexed');figure,imshow(gms);
title('中值滤波 indexed,PSNR=39.159');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);

psnr1 = getPSNR(f,gm);psnr2 = getPSNR(f,gms);psnr3 = getPSNR(f,gmi);

%线性空间域滤波
%先用函数fspecial生成滤波模板w，再由函数imfilter实现线性空间滤波，用replicate
w1 = fspecial('average',[3,3]);
g1 = imfilter(f,w1,'replicate');
figure;imshow(g1,[]);title('矩形平均滤波器，PSNR=39.602');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr4 = getPSNR(f,g1);

w2 = fspecial('disk',5);
g2 = imfilter(f,w2,'replicate');
figure;imshow(g2,[]);title('圆形平均滤波器，PSNR=32.463');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr5 = getPSNR(f,g2);

w3 = fspecial('gaussian',[3,3],0.5);
g3 = imfilter(f,w3,'replicate');
figure;imshow(g3,[]);title('高斯低通滤波器，PSNR=48.4455');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr6 = getPSNR(f,g3);

w4 = fspecial('laplacian',0.2);
g4 = imfilter(f,w4,'replicate');
figure;imshow(g4,[]);title('拉普拉斯滤波器，PSNR=6.614');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr7 = getPSNR(f,g4);

w5 = fspecial('log',[5,5],0.5);
g5 = imfilter(f,w5,'replicate');
figure;imshow(g5,[]);title('高斯-拉普拉斯滤波器，PSNR=6.766');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr8 = getPSNR(f,g5);

w6 = fspecial('motion',9,0);
g6 = imfilter(f,w6,'replicate');
figure;imshow(g6,[]);title('motion滤波器，PSNR=33.703');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr9 = getPSNR(f,g6);

w7 = fspecial('prewitt');
g7 = imfilter(f,w7,'replicate');
figure;imshow(g7,[]);title('prewitt滤波器，PSNR=6.891');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr10 = getPSNR(f,g7);

w8 = fspecial('sobel');
g8 = imfilter(f,w8,'replicate');
figure;imshow(g8,[]);title('sobel滤波器，PSNR=7.054');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr11 = getPSNR(f,g8);

w9 = fspecial('unsharp',0.2);
g9 = imfilter(f,w9,'replicate');
figure;imshow(g9,[]);title('钝化滤波器，PSNR=29.586');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr12 = getPSNR(f,g9);

%% 频率域滤波器
%dftfilt用于频率域滤波，
%tofloat把输入图像转换为浮点图像，
%paddedzsize获得填充参数P、Q，
%dftuv计算网格数组，
%lpfilter生成低通滤波器的传递函数，ideal-理想低通滤波器，btw-巴特沃兹，gaussian-高斯
%hpfilter生成高通滤波器，其实就是HP=1-LP
%clc,clear,close all
%f = imread('C:\Users\HS\Desktop\第一章\水雷图像-灰度图.jpg');

PQ = paddedsize(size(f));
D0 = 0.05*PQ(1);
h = lpfilter('gaussian',PQ(1),PQ(2),D0,2);%设计低通滤波器
HBW = hpfilter('gaussian',PQ(1),PQ(2),D0,2);%设计高通滤波器
H = 0.5+2*HBW;%补偿方法

gbw = dftfilt(f,HBW);%实现高频滤波
gbw = gscale(gbw);%标定图像灰度
figure,imshow(gbw);title('高频滤波,PSNR=15.717');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr13 = getPSNR(f,gbw);

ghf = dftfilt(f,H);%实现高频强调滤波
ghf = gscale(ghf);
figure,imshow(ghf);title('高频强调滤波,PSNR=18.524');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr14 = getPSNR(f,ghf);

g = dftfilt(f,h);%实现低频滤波
g = gscale(g);
figure,imshow(g);title('低频滤波,PSNR=17.905');
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr15 = getPSNR(f,g);

%% 小波去噪
%用小波函数sym4对x进行2层小波分解  
[c,s]=wavedec2(f,2,'sym4');  
%提取小波分解中第一层的低频图像，即实现了低通滤波去噪  
a1=wrcoef2('a',c,s,'sym4');  % a1为 double 型数据；
%画出去噪后的图像  
figure();imshow(uint8(a1)); 
% 注意 imshow()和image()显示图像有区别，imshow()不能显示 double 型数据，必须进行转换 uint8(a1)；
title('第一次去噪图像,PSNR=34.927');           % 并且image() 显示图像有坐标；
set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
%提取小波分解中第二层的低频图像，即实现了低通滤波去噪  
%相当于把第一层的低频图像经过再一次的低频滤波处理  
a2=wrcoef2('a',c,s,'sym4',2);  
%画出去噪后的图像   
figure();imshow(uint8(a2)); %image(a2);
title('第二次去噪图像,PSNR=34.927');set(gcf,'color','w');set(gca,'linewidth',1,'fontsize',16);
psnr16 = getPSNR(f,a1);psnr17 = getPSNR(f,a2);

下面是几个需要调用的函数：

hpfilter生成高通滤波器

function H = hpfilter(type, M, N, D0, n)
%HPFILTER Computes frequency domain highpass filters.
%   H = HPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a highpass filter, H, of the specified TYPE and size (M-by-N).
%   Valid values for TYPE, D0, and n are: 
%
%   'ideal'    Ideal highpass filter with cutoff frequency D0.  n
%              need not be supplied. D0 must be positive.
%
%   'btw'      Butterworth highpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0. D0 must be
%              positive.
%
%   'gaussian' Gaussian highpass filter with cutoff (standard
%              deviation) D0. n need not be supplied. D0 must be
%              positive.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/08/25 14:28:22 $

% The transfer function Hhp of a highpass filter is 1 - Hlp, 
% where Hlp is the transfer function of the corresponding lowpass 
% filter.  Thus, we can use function lpfilter to generate highpass 
% filters.

if nargin == 4
   n = 1; % Default value of n.
end

% Generate highpass filter.
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;

lpfilter生成低通滤波器的传递函数

function H = lpfilter(type, M, N, D0, n)
%LPFILTER Computes frequency domain lowpass filters.
%   H = LPFILTER(TYPE, M, N, D0, n) creates the transfer function of
%   a lowpass filter, H, of the specified TYPE and size (M-by-N). To
%   view the filter as an image or mesh plot, it should be centered
%   using H = fftshift(H). 
%
%   Valid values for TYPE, D0, and n are:
%
%   'ideal'    Ideal lowpass filter with cutoff frequency D0. n need
%              not be supplied.  D0 must be positive.
%
%   'btw'      Butterworth lowpass filter of order n, and cutoff
%              D0.  The default value for n is 1.0.  D0 must be
%              positive.
%
%   'gaussian' Gaussian lowpass filter with cutoff (standard
%              deviation) D0.  n need not be supplied.  D0 must be
%              positive. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.8 $  $Date: 2004/11/04 22:33:16 $

% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances. 
[U, V] = dftuv(M, N);

% Compute the distances D(U, V).
D = sqrt(U.^2 + V.^2);

% Begin filter computations.
switch type
case 'ideal'
   H = double(D <= D0);
case 'btw'
   if nargin == 4
      n = 1;	
   end
   H = 1./(1 + (D./D0).^(2*n));
case 'gaussian'
   H = exp(-(D.^2)./(2*(D0^2)));
otherwise
   error('Unknown filter type.')
end

paddedzsize获得填充参数

function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering. 
%   PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
%   computes the two-element size vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
%   vectors, computes the two-element size vector PQ.  The elements
%   of PQ are the smallest even integers greater than or equal to 
%   AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
%   PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]). 
    
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $

if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 & ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.
   
   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
elseif nargin == 3
   m = max([AB CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P, P];
else 
   error('Wrong number of inputs.')
end

tofloat把输入图像转换为浮点图像

function [out, revertclass]=tofloat(in)
%tofloat convert image to floating point
 
identity=@(x) x;
tosingle=@im2single;
 
table={'uint8', tosingle, @im2uint8
    'uint16', tosingle, @im2uint16
    'int16', tosingle, @im2int16
    'logical', tosingle, @logical
    'double', identity, identity
    'single', identity, identity};
classIndex=find(strcmp(class(in),table(:,1)));
if isempty(classIndex)
error('Unsupported inut image class.')
end
 
out=table{classIndex,2}(in);
revertclass=table{classIndex,3};

dftfilt用于频率域滤波

function g = dftfilt(f, H)
%DFTFILT Performs frequency domain filtering.
%   G = DFTFILT(F, H) filters F in the frequency domain using the
%   filter transfer function H. The output, G, is the filtered
%   image, which has the same size as F.  DFTFILT automatically pads
%   F to be the same size as H.  Function PADDEDSIZE can be used to
%   determine an appropriate size for H.
%
%   DFTFILT assumes that F is real and that H is a real, uncentered
%   circularly-symmetric filter function. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/08/25 14:28:22 $
%[f.revertClass] = tofloat(f);
% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));

% Perform filtering. 
g = real(ifft2(H.*F));

% Crop to original size.
g = g(1:size(f, 1), 1:size(f, 2));

%convert the output to the same class as the inpt if so specified
%if nargin==2||strcmp(classout,'original')
%    g = revertClass(g);
%elseif strcmp(classout,'fltpoint')
%    return
%else
%    error('Undefined class for the outpt image.')
%end

gscale标定图像灰度

function g = gscale(f, varargin)
%GSCALE Scales the intensity of the input image.
%   G = GSCALE(F, 'full8') scales the intensities of F to the full
%   8-bit intensity range [0, 255].  This is the default if there is
%   only one input argument.
%
%   G = GSCALE(F, 'full16') scales the intensities of F to the full
%   16-bit intensity range [0, 65535]. 
%
%   G = GSCALE(F, 'minmax', LOW, HIGH) scales the intensities of F to
%   the range [LOW, HIGH]. These values must be provided, and they
%   must be in the range [0, 1], independently of the class of the
%   input. GSCALE performs any necessary scaling. If the input is of
%   class double, and its values are not in the range [0, 1], then
%   GSCALE scales it to this range before processing.
%
%   The class of the output is the same as the class of the input.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:36:09 $

if length(varargin) == 0 % If only one argument it must be f.
   method = 'full8';
else 
   method = varargin{1};
end

if strcmp(class(f), 'double') & (max(f(:)) > 1 | min(f(:)) < 0)
   f = mat2gray(f);
end

% Perform the specified scaling.
switch method
case 'full8'
   g = im2uint8(mat2gray(double(f)));
case 'full16'
   g = im2uint16(mat2gray(double(f)));
case 'minmax'
   low = varargin{2}; high = varargin{3};
   if low > 1 | low < 0 | high > 1 | high < 0
      error('Parameters low and high must be in the range [0, 1].')
   end
   if strcmp(class(f), 'double')
      low_in = min(f(:));
      high_in = max(f(:));
   elseif strcmp(class(f), 'uint8')
      low_in = double(min(f(:)))./255;
      high_in = double(max(f(:)))./255;
   elseif strcmp(class(f), 'uint16')
      low_in = double(min(f(:)))./65535;
      high_in = double(max(f(:)))./65535;    
   end
   % imadjust automatically matches the class of the input.
   g = imadjust(f, [low_in high_in], [low high]);   
otherwise
   error('Unknown method.')
end

getPSNR计算峰值信噪比

function [PSNR]=getPSNR(I,R)
% 本算法用于计算去噪的客观标准PSNR（越大越好）
%I---->原始无噪声图像
%R---->噪声恢复后的图像
%PSNR---->峰值信噪比，PSNR值越大，说明滤波效果越好
[rowNum,colNum]=size(I);
CNT=rowNum*colNum;
I=double(I);
R=double(R);
MSE=sum(sum((R-I).^2))/CNT; 
PSNR=10*log10(255^2/MSE);

还是建议大家去买实体书看看，虽然翻译得很糟糕呵呵哒~
全套代码有个大佬提供过了，可以自己下载下来搞一哈：
冈萨雷斯的《数字信号处理》免费资源