% GABORCONVOLVE - function for convolving image with log-Gabor filters % % Usage: [EO, BP] = gaborconvolve(im, nscale, norient, minWaveLength, mult, ... % sigmaOnf, dThetaOnSigma, Lnorm, feedback) % % Arguments: % The convolutions are done via the FFT. Many of the parameters relate % to the specification of the filters in the frequency plane. % % Variable Suggested Description % name value % ---------------------------------------------------------- % im Image to be convolved. % 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. % Lnorm 0 Optional integer indicating what norm the % filters should be normalized to. A value of 1 % will produce filters with the same L1 norm, 2 % will produce filters with matching L2 % norm. the default value of 0 results in no % normalization (the filters have unit height % Gaussian transfer functions on a log frequency % scale) % feedback 0/1 Optional parameter. If set to 1 a message % indicating which orientation is being % processed is printed on the screen. % % Returns: % % EO - 2D cell array of complex valued convolution results % % EO{s,o} = convolution result for scale s and orientation o. % The real part is the result of convolving with the even % symmetric filter, the imaginary part is the result from % convolution with the odd symmetric filter. % % Hence: % abs(EO{s,o}) returns the magnitude of the convolution over the % image at scale s and orientation o. % angle(EO{s,o}) returns the phase angles. % % BP - Cell array of bandpass images corresponding to each scale s. % % % Notes on filter settings to obtain even coverage of the spectrum energy % dThetaOnSigma 1.2 - 1.3 % sigmaOnf .90 mult 1.15 % sigmaOnf .85 mult 1.2 % sigmaOnf .75 mult 1.4 (bandwidth ~1 octave) % sigmaOnf .65 mult 1.7 % sigmaOnf .55 mult 2.2 (bandwidth ~2 octaves) % % The determination of mult given sigmaOnf is entirely empirical. What I do is % plot out the sum of the squared filter amplitudes in the frequency domain and % see how even the coverage of the spectrum is. If there are concentric 'gaps' % in the spectrum one needs to reduce mult and/or reduce sigmaOnf (which % increases filter bandwidth) % % If there are 'gaps' radiating outwards then one needs to reduce dthetaOnSigma % (increasing angular bandwidth of the filters) % % For details of log-Gabor filters see: % D. J. Field, "Relations Between the Statistics of Natural Images and the % Response Properties of Cortical Cells", Journal of The Optical Society of % America A, Vol 4, No. 12, December 1987. pp 2379-2394 % Copyright (c) 2001-2010 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 % April 2010 - Reworked to tidy things up. Return of bandpass images added. % March 2013 - Restored use of dThetaOnSigma to control angular spread of filters. function [EO, BP] = gaborconvolve(im, nscale, norient, minWaveLength, mult, ... sigmaOnf, dThetaOnSigma, Lnorm, feedback) if ndims(im) == 3 warning('Colour image supplied: Converting to greyscale'); im = rgb2gray(im); end if ~exist('Lnorm','var'), Lnorm = 0; end if ~exist('feedback','var'), feedback = 0; end if ~isa(im,'double'), im = double(im); end [rows cols] = size(im); imagefft = fft2(im); % Fourier transform of image EO = cell(nscale, norient); % Pre-allocate cell array BP = cell(nscale,1); % 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. 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. This keeps the overall norm of each filter not too dissimilar. lp = lowpassfilter([rows,cols],.45,15); % Radius .45, 'sharpness' 15 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). % Compute bandpass image for each scale if Lnorm == 2 % Normalize filters to have same L2 norm L = sqrt(sum(logGabor{s}(:).^2)); elseif Lnorm == 1 % Normalize to have same L1 L = sum(sum(abs(real(ifft2(logGabor{s}))))); elseif Lnorm == 0 % No normalization L = 1; else error('Lnorm must be 0, 1 or 2'); end logGabor{s} = logGabor{s}./L; BP{s} = ifft2(imagefft .* logGabor{s}); end % The main loop... for o = 1:norient, % For each orientation. if feedback fprintf('Processing orientation %d \r', o); end angl = (o-1)*pi/norient; % Calculate filter angle. wavelength = minWaveLength; % Initialize filter wavelength. % Pre-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. % Calculate the standard deviation of the angular Gaussian % function used to construct filters in the freq. plane. thetaSigma = pi/norient/dThetaOnSigma; spread = exp((-dtheta.^2) / (2 * thetaSigma^2)); for s = 1:nscale, % For each scale. filter = logGabor{s} .* spread; % Multiply by the angular spread to get the filter if Lnorm == 2 % Normalize filters to have the same L2 norm ** why sqrt 2 **???? L = sqrt(sum(real(filter(:)).^2 + imag(filter(:)).^2 ))/sqrt(2); filter = filter./L; elseif Lnorm == 1 % Normalize to have same L1 L = sum(sum(abs(real(ifft2(filter))))); filter = filter./L; elseif Lnorm == 0 % No normalization ; end % Do the convolution, back transform, and save the result in EO EO{s,o} = ifft2(imagefft .* filter); wavelength = wavelength * mult; % Finally calculate Wavelength of next filter end % ... and process the next scale end % For each orientation if feedback, fprintf(' \r'); end