% PHASESYM - Function for computing phase symmetry on an image.%% This function calculates the phase symmetry of points in an image.% This is a contrast invariant measure of symmetry.  This function can be% used as a line and blob detector.  The greyscale 'polarity' of the lines% that you want to find can be specified.%% There are potentially many arguments, here is the full usage:%%   [phaseSym, orientation, totalEnergy] = ...%                phasesym(im, nscale, norient, minWaveLength, mult, ...%                         sigmaOnf, dThetaOnSigma, k, polarity)%% However, apart from the image, all parameters have defaults and the% usage can be as simple as:%%    phaseSym = phasesym(im);% % Arguments:%              Default values      Description%%    nscale           5    - Number of wavelet scales, try values 3-6%    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.%                            You may want to vary this up to a value of 10 or%                            20 for noisy images %    polarity         0    - Controls 'polarity' of symmetry features to find.%                             1 - just return 'bright' points%                            -1 - just return 'dark' points%                             0 - return bright and dark points.%% Return values:%    phaseSym              - Phase symmetry image (values between 0 and 1).%    orientation           - Orientation image. Orientation in which local%                            symmetry energy is a maximum, in degrees%                            (0-180), angles positive anti-clockwise.%    totalEnergy           - Un-normalised raw symmetry energy which may be%                            more to your liking.%%% Notes on specifying parameters:  %% The parameters can be specified as a full list eg.%  >> phaseSym = phasesym(im, 5, 6, 3, 2.5, 0.55, 1.2, 2.0, 0);%% or as a partial list with unspecified parameters taking on default values%  >> phaseSym = phasesym(im, 5, 6, 3);%% or as a partial list of parameters followed by some parameters specified via a% keyword-value pair, remaining parameters are set to defaults, for example:%  >> phaseSym = phasesym(im, 5, 6, 3, 'polarity',-1, 'k', 2.5);% % 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.%% Notes on filter settings to obtain even coverage of the spectrum% dthetaOnSigma 1.2    norient 6% sigmaOnf       .85   mult 1.3% sigmaOnf       .75   mult 1.6     (filter bandwidth ~1 octave)% sigmaOnf       .65   mult 2.1  % sigmaOnf       .55   mult 3       (filter bandwidth ~2 octaves)%% 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.%% See Also:  PHASECONG, PHASECONG2, GABORCONVOLVE, PLOTGABORFILTERS% References:%     Peter Kovesi, "Symmetry and Asymmetry From Local Phase" AI'97, Tenth%     Australian Joint Conference on Artificial Intelligence. 2 - 4 December%     1997. http://www.cs.uwa.edu.au/pub/robvis/papers/pk/ai97.ps.gz.%%     Peter Kovesi, "Image Features From Phase Congruency". Videre: A%     Journal of Computer Vision Research. MIT Press. Volume 1, Number 3,%     Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html% April 1996     Original Version written % August 1998    Noise compensation corrected. % October 1998   Noise compensation corrected.   - Again!!!% September 1999 Modified to operate on non-square images of arbitrary size. % February 2001  Specialised from phasecong.m to calculate phase symmetry % July 2005      Better argument handling + general cleanup and speed improvements% August 2005    Made Octave compatible.% January 2007   Small correction and cleanup of radius calculation for odd%                image sizes.% Copyright (c) 1996-2005 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.function[phaseSym, orientation, totalEnergy] = phasesym(varargin)    % Get arguments and/or default values        [im, nscale, norient, minWaveLength, mult, sigmaOnf,  dThetaOnSigma,k, ...     polarity] = checkargs(varargin(:));      Octave = exist('OCTAVE_VERSION') ~= 0;  % Are we running under Octave?        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    % The main loop...    for o = 1:norient,                   % For each orientation.        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.