% PHASECONG - Computes phase congruency on an image.%% Usage: [pc or ft] = phasecong(im)%% This function calculates the PC_2 measure of phase congruency.  % For maximum speed the input image should be square and have a % size that is a power of 2, but the code will operate on images% of arbitrary size.  %%% Returned values:%                 pc - Phase congruency image (values between 0 and 1)   %                 or - Orientation image.  Provides orientation in which %                      local energy is a maximum in in degrees (0-180), %                      angles positive anti-clockwise.%                 ft - A complex valued image giving the weighted mean %                      phase angle at every point in the image for the %                      orientation having maximum energy.  Use the%                      function DISPFEAT to display this data.%% Parameters:  %                 im - A greyscale image to be processed.%                 % You can also specify numerous optional parameters.  See the code to find% out what they are.  The convolutions are done via the FFT.  Many of the% parameters relate to the specification of the filters in the frequency% plane.  Default values for parameters are set within the file rather than% being required as arguments because they rarely need to be changed - nor% are they very critical.  However, you may want to experiment with% specifying/editing the values of `nscales' and `noiseCompFactor'.%% Note this phase congruency code is very computationally expensive and uses% *lots* of memory.%%% Example MATLAB session:%% >> im = imread('picci.tif');   % >> image(im);                          % Display the image% >> [pc or ft] = phasecong(im);% >> imagesc(pc), colormap(gray);        % Display the phase congruency image%%% To convert the phase congruency image to an edge map (with my usual parameters):%% >> nm = nonmaxsup(pc, or, 1.5);        % Non-maxima suppression.  %    The parameter 1.5 can result in edges more than 1 pixel wide but helps%    in picking up `broad' maxima.% >> edgim = hysthresh(nm, 0.4, 0.2);    % Hysteresis thresholding.% >> edgeim = bwmorph(edgim,'skel',Inf); % Skeletonize the edgemap to fix%                                        % the non-maximal suppression.% >> imagesc(edgeim), colormap(gray);%%% To display the different feature types present in your image use:%% >> dispfeat(ft,edgim);%% With a small amount of editing the code can be modified to calculate% a dimensionless measure of local symmetry in the image.  The basis% of this is that one looks for points in the image where the local % phase is 90 or 270 degrees (the symmetric points in the cycle).% Editing instructions are within the code.%% Notes on filter settings to obtain even coverage of the spectrum% dthetaOnSigma 1.5% sigmaOnf  .85   mult 1.3% sigmaOnf  .75   mult 1.6     (bandwidth ~1 octave)% sigmaOnf  .65   mult 2.1  % sigmaOnf  .55   mult 3       (bandwidth ~2 octaves)%% References:%%     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% 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.% Original Version written April 1996     % Noise compensation corrected. August 1998% Noise compensation corrected. October 1998  - Again!!!% Modified to operate on non-square images of arbitrary size. September 1999% Modified to return feature type image. May 2001function[phaseCongruency,orientation, featType]=phasecong(im, nscale, norient, ...						minWaveLength, mult, ...						sigmaOnf, dThetaOnSigma, ...						k, cutOff)sze = size(im);if nargin < 2    nscale          = 3;     % Number of wavelet scales.endif nargin < 3    norient         = 6;     % Number of filter orientations.endif nargin < 4    minWaveLength   = 3;     % Wavelength of smallest scale filter.endif nargin < 5    mult            = 2;     % Scaling factor between successive filters.endif nargin < 6    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.endif nargin < 7    dThetaOnSigma   = 1.7;   % Ratio of angular interval between filter orientations			     % and the standard deviation of the angular Gaussian			     % function used to construct filters in the                             % freq. plane.endif nargin < 8    k               = 3.0;   % No of standard deviations of the noise energy beyond the			     % mean at which we set the noise threshold point.			     % standard deviation to its maximum effect                             % on Energy.endif nargin < 9    cutOff          = 0.4;   % The fractional measure of frequency spread                             % below which phase congruency values get penalized.end   g               = 10;    % Controls the sharpness of the transition in the sigmoid                         % function used to weight phase congruency for frequency                         % spread.epsilon         = .0001; % Used to prevent division by zero.thetaSigma = pi/norient/dThetaOnSigma;  % Calculate the standard deviation of the                                        % angular Gaussian function used to                                        % construct filters in the freq. plane.imagefft = fft2(im);                    % Fourier transform of imagesze = size(imagefft);rows = sze(1);cols = sze(2);zero = zeros(sze);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 = [];% Pre-compute some stuff to speed up filter constructionx = ones(rows,1) * (-cols/2 : (cols/2 - 1))/(cols/2);  y = (-rows/2 : (rows/2 - 1))' * ones(1,cols)/(rows/2);radius = sqrt(x.^2 + y.^2);       % Matrix values contain *normalised* radius from centre.radius(round(rows/2+1),round(cols/2+1)) = 1; % Get rid of the 0 radius value in the middle                                              % so that taking the log of the radius will                                              % not cause trouble.theta = atan2(-y,x);              % Matrix values contain polar angle.                                  % (note -ve y is used to give +ve                                  % anti-clockwise angles)sintheta = sin(theta);costheta = cos(theta);clear x; clear y; clear theta;      % save a little memory% The main loop...for o = 1:norient,                   % For each orientation.  disp(['Processing orientation ' num2str(o)]);  angl = (o-1)*pi/norient;           % Calculate filter angle.  wavelength = minWaveLength;        % Initialize filter wavelength.  sumE_ThisOrient   = zero;          % Initialize accumulator matrices.  sumO_ThisOrient   = zero;         sumAn_ThisOrient  = zero;        Energy_ThisOrient = zero;        EOArray = [];          % Array of complex convolution images - one for each scale.  ifftFilterArray = [];  % Array of inverse FFTs of filters  % 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.