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.                        if s==1                EM_n = sum(sum(filter.^2)); % Record mean squared filter value at smallest            end                             % scale. This is used for noise estimation.        end                                 % ... and process the next scale        % Now calculate the phase symmetry measure.        if polarity == 0     % look for 'white' and 'black' spots            for s = 1:nscale,                                  Energy_ThisOrient = Energy_ThisOrient ...                    + abs(real(EO{s,o})) - abs(imag(EO{s,o}));            end                    elseif polarity == 1  % Just look for 'white' spots            for s = 1:nscale,                                  Energy_ThisOrient = Energy_ThisOrient ...                    + real(EO{s,o}) - abs(imag(EO{s,o}));            end                    elseif polarity == -1  % Just look for 'black' spots            for s = 1:nscale,                                  Energy_ThisOrient = Energy_ThisOrient ...                    - real(EO{s,o}) - abs(imag(EO{s,o}));            end                    end                % Compensate for noise        % We estimate the noise power from the energy squared response at the        % smallest scale.  If the noise is Gaussian the energy squared will        % have a Chi-squared 2DOF pdf.  We calculate the median energy squared        % response as this is a robust statistic.  From this we estimate the        % mean.  The estimate of noise power is obtained by dividing the mean        % squared energy value by the mean squared filter value                medianE2n = median(reshape(abs(EO{1,o}).^2,1,rows*cols));        meanE2n = -medianE2n/log(0.5);        estMeanE2n = [estMeanE2n meanE2n];                noisePower = meanE2n/EM_n;                       % Estimate of noise power.                % Now estimate the total energy^2 due to noise        % Estimate for sum(An^2) + sum(Ai.*Aj.*(cphi.*cphj + sphi.*sphj))                EstSumAn2 = zero;        for s = 1:nscale            EstSumAn2 = EstSumAn2+ifftFilterArray{s}.^2;        end                EstSumAiAj = zero;        for si = 1:(nscale-1)            for sj = (si+1):nscale                EstSumAiAj = EstSumAiAj + ifftFilterArray{si}.*ifftFilterArray{sj};            end        end                EstNoiseEnergy2 = 2*noisePower*sum(sum(EstSumAn2)) + 4*noisePower*sum(sum(EstSumAiAj));                tau = sqrt(EstNoiseEnergy2/2);                % Rayleigh parameter        EstNoiseEnergy = tau*sqrt(pi/2);              % Expected value of noise energy        EstNoiseEnergySigma = sqrt( (2-pi/2)*tau^2 );                T =  EstNoiseEnergy + k*EstNoiseEnergySigma;  % Noise threshold                % The estimated noise effect calculated above is only valid for the PC_1        % measure.  The PC_2 measure does not lend itself readily to the same        % analysis.  However empirically it seems that the noise effect is        % overestimated roughly by a factor of 1.7 for the filter parameters        % used here.        T = T/1.7;            % Apply noise threshold         Energy_ThisOrient = max(Energy_ThisOrient - T, zero);                                   % Update accumulator matrix for sumAn and totalEnergy        totalSumAn  = totalSumAn + sumAn_ThisOrient;        totalEnergy = totalEnergy + Energy_ThisOrient;                % Update orientation matrix by finding image points where the energy in        % this orientation is greater than in any previous orientation (the        % change matrix) and then replacing these elements in the orientation        % matrix with the current orientation number.                if(o == 1),            maxEnergy = Energy_ThisOrient;        else            change = Energy_ThisOrient > maxEnergy;            orientation = (o - 1).*change + orientation.*(~change);            maxEnergy = max(maxEnergy, Energy_ThisOrient);        end            end  % For each orientation    fprintf('                                   \r');        %    disp('Mean Energy squared values recorded with smallest scale filter at each orientation');%    disp(estMeanE2n);        % Normalize totalEnergy by the totalSumAn to obtain phase symmetry    phaseSym = totalEnergy ./ (totalSumAn + epsilon);        % Convert orientation matrix values to degrees    orientation = orientation * (180 / norient);        %------------------------------------------------------------------% CHECKARGS%% Function to process the arguments that have been supplied, assign% default values as needed and perform basic checks.    function [im, nscale, norient, minWaveLength, mult, sigmaOnf, ...          dThetaOnSigma,k, polarity] = checkargs(arg);     nargs = length(arg);        if nargs < 1        error('No image supplied as an argument');    end            % Set up default values for all arguments and then overwrite them    % with with any new values that may be supplied    im              = [];    nscale          = 5;     % Number of wavelet scales.        norient         = 6;     % Number of filter orientations.    minWaveLength   = 3;     % Wavelength of smallest scale filter.        mult            = 2.1;   % Scaling factor between successive filters.        sigmaOnf        = 0.55;  % 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.2;   % Ratio of angular interval between filter orientations                                 % and the standard deviation of the angular Gaussian                             % function used to construct filters in the                             % freq. plane.    k               = 2.0;   % No of standard deviations of the noise                             % energy beyond the mean at which we set the                             % noise threshold point.     polarity        = 0;     % Look for both black and white spots of symmetrry        % Allowed argument reading states    allnumeric   = 1;       % Numeric argument values in predefined order    keywordvalue = 2;       % Arguments in the form of string keyword                            % followed by numeric value    readstate = allnumeric; % Start in the allnumeric state        if readstate == allnumeric        for n = 1:nargs            if isa(arg{n},'char')                readstate = keywordvalue;                break;            else                if     n == 1, im            = arg{n};                 elseif n == 2, nscale        = arg{n};                              elseif n == 3, norient       = arg{n};                elseif n == 4, minWaveLength = arg{n};                elseif n == 5, mult          = arg{n};                elseif n == 6, sigmaOnf      = arg{n};                elseif n == 7, dThetaOnSigma = arg{n};                elseif n == 8, k             = arg{n};                              elseif n == 9, polarity      = arg{n};                                              end            end        end    end    % Code to handle parameter name - value pairs    if readstate == keywordvalue        while n < nargs                        if ~isa(arg{n},'char') | ~isa(arg{n+1}, 'double')                error('There should be a parameter name - value pair');            end                        if     strncmpi(arg{n},'im'      ,2), im =        arg{n+1};            elseif strncmpi(arg{n},'nscale'  ,2), nscale =    arg{n+1};            elseif strncmpi(arg{n},'norient' ,2), norient =   arg{n+1};            elseif strncmpi(arg{n},'minWaveLength',2), minWaveLength = arg{n+1};            elseif strncmpi(arg{n},'mult'    ,2), mult =      arg{n+1};            elseif strncmpi(arg{n},'sigmaOnf',2), sigmaOnf =  arg{n+1};            elseif strncmpi(arg{n},'dthetaOnSigma',2), dThetaOnSigma =  arg{n+1};            elseif strncmpi(arg{n},'k'       ,1), k =         arg{n+1};            elseif strncmpi(arg{n},'polarity',2), polarity =  arg{n+1};            else   error('Unrecognised parameter name');            end            n = n+2;            if n == nargs                error('Unmatched parameter name - value pair');            end                    end    end        if isempty(im)        error('No image argument supplied');    end    if ~isa(im, 'double')        im = double(im);    end        if nscale < 1        error('nscale must be an integer >= 1');    end        if norient < 1         error('norient must be an integer >= 1');    end        if minWaveLength < 2        error('It makes little sense to have a wavelength < 2');    end                      if polarity ~= -1 & polarity ~= 0 & polarity ~= 1        error('Allowed polarity values are -1, 0 and 1')    end