% 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






