% gaborWavelet - Function for computing gabor features of a gray-scale image
%
% This function calculates gabor features. Mean-squared energy & meanAmplitude
% for each scale % and orientation is returned.
%
% There are potentially many arguments, here is the full usage:
%
% [gaborSquareEnergy, gaborMeanAmplitude] = ...
% gaborWavelet(im, nscale, norient )
%NOTE: nscale & norient are optional arguments
%
% However, apart from the image, all parameters have defaults and the
% usage can be as simple as:
%
% [gaborSquareEnergy, gaborMeanAmplitude ]= gaborWavelet(im);
%
% Arguments:
% Default values Description
%
% nscale 5 - Number of wavelet scales, try values 3-6
% norient 6 - Number of filter orientations.
%
% Return values:
% msEnergy - Mean square energy
% orientation - Mean amplitude
%
% The convolutions are done via the FFT. Many of the parameters relate to the
% specification of the filters in the frequency plane. The values do not seem
% to be very critical and the defaults are usually fine. You may want to
% experiment with the values of 'nscales' and 'k', the noise compensation factor.
%
% For maximum speed the input image should have dimensions that correspond to
% powers of 2, but the code will operate on images of arbitrary size.
function[gaborSquareEnergy, gaborMeanAmplitude] = gaborWavelet(varargin)
% Get arguments and/or default values
[im, nscale, norient, minWaveLength, mult, sigmaOnf, dThetaOnSigma,k, ...
polarity] = checkargs(varargin(:));
v = version; Octave = v(1) < '5'; % Crude Octave test
epsilon = .0001; % Used to prevent division by zero.
% Calculate the standard deviation of the angular Gaussian function
% used to construct filters in the frequency plane.
thetaSigma = pi/norient/dThetaOnSigma;
[rows,cols] = size(im);
imagefft = fft2(im); % Fourier transform of image
zero = zeros(rows,cols);
totalEnergy = zero; % Matrix for accumulating weighted phase
% congruency values (energy).
totalSumAn = zero; % Matrix for accumulating filter response
% amplitude values.
orientation = zero; % Matrix storing orientation with greatest
% energy for each pixel.
estMeanE2n = [];
EO = cell(nscale, norient); % Cell array of convolution results
ifftFilterArray = cell(1, nscale); % Cell array of inverse FFTs of filters
% 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); % Matrix values contain *normalised* radius from centre.
theta = atan2(-y,x); % Matrix values contain polar angle.
% (note -ve y is used to give +ve
% anti-clockwise angles)
radius = ifftshift(radius); % Quadrant shift radius and theta so that filters
theta = ifftshift(theta); % are constructed with 0 frequency at the corners.
radius(1,1) = 1; % Get rid of the 0 radius value at the 0
% frequency point (now at top-left corner)
% so that taking the log of the radius will
% not cause trouble.
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta; % save a little memory
% Filters are constructed in terms of two components.
% 1) The radial component, which controls the frequency band that the filter
% responds to
% 2) The angular component, which controls the orientation that the filter
% responds to.
% The two components are multiplied together to construct the overall filter.
% Construct the radial filter components...
% First construct 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 no extra frequencies at the 'corners' of the FFT are
% incorporated as this seems to upset the normalisation process when
% calculating phase congrunecy.
lp = lowpassfilter([rows,cols], .4, 10); % Radius .4, 'sharpness' 10
logGabor = cell(1, nscale);
for s = 1:nscale
wavelength = minWaveLength*mult^(s-1);
fo = 1.0/wavelength; % Centre frequency of filter.
logGabor{s} = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));
logGabor{s} = logGabor{s}.*lp; % Apply low-pass filter
logGabor{s}(1, 1) = 0; % Set the value at the 0 frequency point of the filter
% back to zero (undo the radius fudge).
end
% Then construct the angular filter components...
spread = cell(1, norient);
for o = 1:norient
angl = (o-1)*pi/norient; % Filter angle.
% 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{o} = exp((-dtheta.^2) / (2 * thetaSigma^2)); % Calculate the
% angular filter component.
end
count = 1;
gaborSquareEnergy = [];
gaborMeanAmplitude = [];
% The main loop...
for o = 1:norient, % For each orientation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fprintf('Processing orientation %d \r', o);
if Octave fflush(1); end
sumAn_ThisOrient = zero;
Energy_ThisOrient = zero;
for s = 1:nscale, % For each scale.
filter = logGabor{s} .* spread{o}; % Multiply radial and angular
% components to get filter.
ifftFilt = real(ifft2(filter))*sqrt(rows*cols); % Note rescaling to match power
ifftFilterArray{s} = ifftFilt; % record ifft2 of filter
% Convolve image with even and odd filters returning the result in EO
EO{s, o} = ifft2(imagefft .* filter);
An = abs(EO{s,o}); % Amplitude of even & odd filter response.
sumAn_ThisOrient = sumAn_ThisOrient + An; % Sum of amplitude responses.
% % Display of individual components
% if o == 6
% figure, imagesc( An ), colormap(gray), title
