plotgaborfilters.m

在MATLAB中调用gabor变换的代码

% PLOTGABORFILTERS - Plots log-Gabor filters
%
% The purpose of this code is to see what effect the various parameter
% settings have on the formation of a log-Gabor filter bank.
%
% Usage: [Ffilter, Efilter, Ofilter] = plotgaborfilters(sze,  nscale, norient,...
%                                       minWaveLength, mult, sigmaOnf, dThetaOnSigma)
%
% Arguments:
% Many of the parameters relate to the specification of the filters in the frequency plane.  
%
%   Variable       Suggested   Description
%   name           value
%  ----------------------------------------------------------
%    sze             = 200     Size of image grid on which the filters
%                              are calculated.  Note that the actual size
%                              of the filter is really specified by its
%                              wavelength. 
%    nscale          = 4;      Number of wavelet scales.
%    norient         = 6;      Number of filter orientations.
%    minWaveLength   = 3;      Wavelength of smallest scale filter.
%    mult            = 2;      Scaling factor between successive filters.
%    sigmaOnf        = 0.65;   Ratio of the standard deviation of the
%                              Gaussian describing the log Gabor filter's
%                              transfer function in the frequency domain
%                              to the filter center frequency. 
%    dThetaOnSigma   = 1.5;    Ratio of angular interval between filter
%                              orientations and the standard deviation of
%                              the angular Gaussian function used to
%                              construct filters in the freq. plane.
%
% Returns:
%    Ffilter - a 2D cell array of filters defined in the frequency domain.
%    Efilter - a 2D cell array of even filters defined in the spatial domain.
%    Ofilter - a 2D cell array of odd filters defined in the spatial domain.
%
%    Ffilter{s,o} = filter for scale s and orientation o.
%    The even and odd filters in the spatial domain for scale s,
%    orientation o, are obtained using.              
%
%    Efilter = ifftshift(real(ifft2(fftshift(filter{s,o}))));
%    Ofilter = ifftshift(imag(ifft2(fftshift(filter{s,o}))));
%
% Plots:
%    Figure 1         -  Sum of the filters in the frequency domain
%    Figure 2         -  Cross sections of Figure 1
%    Figures 3 and 4  -  Surface and intensity plots of filters in the
%                        spatial domain at the smallest and largest
%                        scales respectively.
%
% See also: GABORCONVOLVE, PHASECONG

% Copyright (c) 2001-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% May 2001      - Original version.
% February 2005 - Cleaned up.
% August   2005 - Ffilter,Efilter and Ofilter corrected to return with scale
%                 varying as the first index in the cell arrays.

function [Ffilter, Efilter, Ofilter] = ...
	plotgaborfilters(sze, nscale, norient, minWaveLength, mult, ...
				   sigmaOnf, dThetaOnSigma) 

    rows = sze; cols = sze;
thetaSigma = pi/norient/dThetaOnSigma;  % Calculate the standard deviation of the
                                        % angular Gaussian function used to
                                        % construct filters in the freq. plane.    
% Pre-compute some stuff to speed up filter construction

% Set up X and Y matrices with ranges normalised to +/- 0.5
% The following code adjusts things appropriately for odd and even values
% of rows and columns.
if mod(cols,2)
    xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1);
else
    xrange = [-cols/2:(cols/2-1)]/cols;	
end

if mod(rows,2)
    yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1);
else
    yrange = [-rows/2:(rows/2-1)]/rows;	
end

[x,y] = meshgrid(xrange, yrange);

radius = sqrt(x.^2 + y.^2);     % Normalised radius (frequency) values 0.0 - 0.5

% Get rid of the 0 radius value in the middle so that taking the log of
% the radius will not cause trouble.
radius(fix(rows/2+1),fix(cols/2+1)) = 1; 

theta = atan2(-y,x);              % Matrix values contain polar angle.
                                  % (note -ve y is used to give +ve
                                  % anti-clockwise angles)
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta;      % save a little memory

% Define a low-pass filter that is as large as possible, yet falls away to zero
% at the boundaries.  All log Gabor filters are multiplied by this to ensure
% that filters are as similar as possible across orientations (Eliminate the
% extra frequencies at the 'corners' of the FFT)
lp = fftshift(lowpassfilter([rows,cols],.45,15));   % Radius .4, 'sharpness' 10

% The main loop...

filtersum = zeros(rows,cols);

for o = 1:2*norient,                 % For each orientation.
  angl = (o-1)*pi/norient;           % Calculate filter angle.
  wavelength = minWaveLength;        % Initialize filter wavelength.

  % Compute filter data specific to this orientation
  % For each point in the filter matrix calculate the angular distance from the
  % specified filter orientation.  To overcome the angular wrap-around problem
  % sine difference and cosine difference values are first computed and then
  % the atan2 function is used to determine angular distance.

  ds = sintheta * cos(angl) - costheta * sin(angl); % Difference in sine.
  dc = costheta * cos(angl) + sintheta * sin(angl); % Difference in cosine.
  dtheta = abs(atan2(ds,dc));                       % Absolute angular distance.
  spread = exp((-dtheta.^2) / (2 * thetaSigma^2));  % The angular filter component.

  for s = 1:nscale,                  % For each scale.

    % Construct the filter - first calculate the radial filter component.
    fo = 1.0/wavelength;                  % Centre frequency of filter.

    logGabor = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));  
    logGabor(round(rows/2+1),round(cols/2+1)) = 0; % Set value at center of the filter
                                                   % back to zero (undo the radius fudge).

    logGabor = logGabor.*lp;             % Apply low-pass filter 
    Ffilter{s,o} = logGabor .* spread;   % Multiply by the angular
                                         % spread to get the filter.
					       
    filtersum = filtersum + Ffilter{s,o};
        
    Efilter{s,o} = ifftshift(real(ifft2(fftshift(Ffilter{s,o}))));
    Ofilter{s,o} = ifftshift(imag(ifft2(fftshift(Ffilter{s,o}))));
    					  
    wavelength = wavelength*mult;
  end
end

% Plot sum of filters and slices radially and tangentially
figure(1), clf, show(filtersum,1), title('sum of filters');

figure(2), clf
subplot(2,1,1), plot(filtersum(round(rows/2+1),:))
title('radial slice through sum of filters');

ang = [0:pi/32:2*pi];
r = rows/4;
tslice = improfile(filtersum,r*cos(ang)+cols/2,r*sin(ang)+rows/2);
subplot(2,1,2), plot(tslice), axis([0 length(tslice) 0 1.1*max(tslice)]);
title('tangential slice through sum of filters at f = 0.25');	   

% Plot Even and Odd filters at the largest and smallest scales
h = figure(3); clf
set(h,'name',sprintf('Filters: Wavelenth = %.2f',minWaveLength));
subplot(3,2,1), surfl(Efilter{1,1}), shading interp, colormap(gray), 
title('Even Filter');
subplot(3,2,2), surfl(Ofilter{1,1}), shading interp, colormap(gray)
title('Odd Filter');
subplot(3,2,3),imagesc(Efilter{1,1}), axis image, colormap(gray)
subplot(3,2,4),imagesc(Ofilter{1,1}), axis image, colormap(gray)
subplot(3,2,5),imagesc(Ffilter{1,1}), axis image, colormap(gray)
title('Frequency Domain');

h = figure(4); clf
set(h,'name',sprintf('Filters: Wavelenth = %.2f',minWaveLength*mult^(nscale-1)));
subplot(3,2,1), surfl(Efilter{nscale,1}), shading interp, colormap(gray)
title('Even Filter');
subplot(3,2,2), surfl(Ofilter{nscale,1}), shading interp, colormap(gray)
title('Odd Filter');
subplot(3,2,3),imagesc(Efilter{nscale,1}), axis image, colormap(gray)
subplot(3,2,4),imagesc(Ofilter{nscale,1}), axis image, colormap(gray)
subplot(3,2,5),imagesc(Ffilter{nscale,1}), axis image, colormap(gray)
title('Frequency Domain');