phasecong.m

hopfield neural network for binary image recognition

% 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 2001



function[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.
end
if nargin < 3
    norient         = 6;     % Number of filter orientations.
end
if nargin < 4
    minWaveLength   = 3;     % Wavelength of smallest scale filter.
end
if nargin < 5
    mult            = 2;     % Scaling factor between successive filters.
end
if 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.
end
if 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.
end
if 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.
end
if 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 image
sze = 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 construction

x = 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.