% 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 = 1.7; 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.3; Ratio of angular interval between filter
% orientations and the standard deviation of
% the angular Gaussian function used to
% construct filters in the freq. plane.
%
% Note regarding the specification of norient: In the default case it is assumed
% that the angles of the filters are evenly spaced at intervals of pi/norient
% around the frequency plane. If you want to visualize filters at a specific
% set of orientations that are not necessarily evenly spaced you can set the
% orientations by passing a CELL array of orientations as the argument to
% norient. In this case the value supplied for dThetaOnSigma will be used as
% thetaSigma - the angular standard deviation of the filters. Yes, this is
% an ugly abuse of the argument list, but there it is!
% Example:
% View filters over 3 scales with orientations of -0.3 and +0.3 radians,
% minWaveLength of 6, mult of 2, sigmaOnf of 0.65 and thetaSigma of 0.4
% plotgaborfilters2(200, 3, {-.3 .3}, 6, 2, 0.65, 0.4);
%
% 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-2008 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.
% July 2008 - Allow specific filter orientations to be specified in
% norient via a cell array.
% March 2013 - Restored use of dThetaOnSigma to control angular spread of filters.
function [Ffilter, Efilter, Ofilter, filtersum] = ...
plotgaborfilters(sze, nscale, norient, minWaveLength, mult, ...
sigmaOnf, dThetaOnSigma)
if ~exist('dThetaOnSigma','var') || isempty(dThetaOnSigma), dThetaOnSigma = 1.3; end
if length(sze) == 1
rows = sze; cols = sze;
else
rows = sze(1); cols = sze(2);
end
if iscell(norient) % Filter orientations and spread have been specified
% explicitly
filterOrient = cell2mat(norient);
thetaSigma = dThetaOnSigma; % Use dThetaOnSigma directly as thetaSigma
norient = length(filterOrient);
else % Usual setup with filters evenly oriented
filterOrient = [0 : pi/norient : pi-pi/norient];
% Calculate the standard deviation of the angular Gaussian function
% used to construct filters in the frequency plane.
thetaSigma = pi/norient/dThetaOnSigma;
end
% Double up all the filter orientations by adding another set offset by
% pi. This allows us to see the overall orientation coverage of the
% filters a bit more easily.
filterOrient = [filterOrient filterOrient+pi];
% 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,10)); % Radius .4, 'sharpness' 10
% The main loop...
filtersum = zeros(rows,cols);
for o = 1:2*norient, % For each orientation.
angl = filterOrient(o);
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.
if ~dThetaOnSigma % If set to 0 use cosine angular weight spread
% function. (I do not think this is a good spread function)
