傅立叶变换滤波

接上文:频域滤波之傅立叶变换的一个例子

滤波步骤

1. 获得填充大小：PQ=paddedSize(size(f))
2. 获得填充后的傅立叶变换 F=fft2(f,PQ(1),PQ(2))
3. 生成一个滤波器，大小为PQ(1)*PQ(2)，需要使用fftshift把重心的波峰移到边角处
4. 频域点乘 G=H.*F
5. 获得逆傅立叶变换出的图像的实部g=real(ifft2(G))
6. 恢复图像原始大小 g=g(1:size(f,1),1:size(f,2))

第3步，分成两种情况先介绍第一种：从空域滤波生成频域滤波

从空域滤波生成频域滤波

下面测试几个基础函数的使用效果

% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %
% 测试几个频域常用的函数的用法
% % % % % % % % % % % % % % % %% % % % % % % % % % % % % % % %

clear all;
close all;
% 测试gscale 测试freqz2 测试fftshift与ifftshift
f=imread('E:\资料\onedrive\code\test\image\Fig0409(a)(bld).tif');
F=fft2(f);
S=fftshift(log(1+abs(F)));
%subplot(1,3,1);imshow(S);
%subplot(1,3,2);imshow(S,[]);
S=gscale(S);   % 测试得这个函数的作用是让频谱自适应地显示出来;
%subplot(1,3,3);imshow(S);
w=1;h=3;i=0;
% 由空域滤波器生成频域滤波器 对比fftshift与ifftshift
H=fspecial('sobel'); %梯度算子
H3=freqz2(H);         %由空域滤波器生成频域滤波器;
i=i+1;subplot(w,h,i);mesh(abs(H3));%生成函数图像？
% figure;
H1=ifftshift(H3);
i=i+1;subplot(w,h,i);mesh(abs(H1));

H2=fftshift(H3);
i=i+1;subplot(w,h,i);mesh(abs(H2));
%view(45,30); %设置看图像的视角;第一个是左右倾斜度，第二个是上下倾斜度
%我就感觉这个ifftshift和fftshift的效果是一样的

1. 测试后发现矩阵经过gscale后的效果与imshow(f,[])的效果一致。但是明显gscale适用范围更广，可以先把这个矩阵标准化，进行相应操作而不是简简单单的输出。
2. freqz2:frequcy表示生成频域滤波，用空域滤波生成频域滤波
3. fftshift:表示把位于中间的波峰移动边上，大概是这个意思？
4. mesh：二维显示频域函数
5. view：调整观测角度

paddedsize函数实现

计算傅立叶滤波需要的填充大小

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

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

dftfilt函数

函数的效果是合并了傅立叶变换滤波中的

2. 填充原始图像并进行傅立叶变换
3. 进行滤波
4. 傅立叶反变换并获得实部
5. 剪切成原始尺寸

所以剩下需要做的就是生成一个合适尺寸的频域滤波器并移动规范化

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 $

% 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));

空域生成频域的例子

clear all;
close all;

f=imread('E:\资料\onedrive\code\test\image\Fig0409(a)(bld).tif');
w=1;h=3;i=0;
h=fspecial('sobel'); %梯度算子
PQ=paddedsize(size(f));
H3=freqz2(h,PQ(1),PQ(2));         %由空域滤波器生成频域滤波器;

H1=fftshift(H3);
g1=dftfilt(f,H1);
g2=imfilter(double(f),h);
myImshow(g1);myImshow(g2);myImshow(0);

运行结果

untitled.png

从结果可以看到空域滤波和频域滤波效果差不多；
但是时间上频域滤波用时：0.044662 seconds.
空域滤波用时.046712 seconds.
这便是频域滤波的优势所在。