% MONOFILT - Apply monogenic filters to an image to obtain 2D analytic signal%% Implementation of Felsberg's monogenic filters%% Usage: [f, h1f, h2f, A, theta, psi] = ...%             monofilt(im, nscale, minWaveLength, mult, sigmaOnf, orientWrap)%                             3         4           2     0.65      1/0% 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          = 3;      Number of filter scales.%    minWaveLength   = 4;      Wavelength of smallest scale filter.%    mult            = 2;      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. %    orientWrap       1/0      Optional flag 1/0 to turn on/off%                              'wrapping' of orientation data from a%                              range of -pi .. pi to the range 0 .. pi.%                              This affects the interpretation of the%                              phase angle - see note below. Defaults to 0.% Returns:%  %        f  - cell array of bandpass filter responses with respect to scale.%      h1f  - cell array of bandpass h1 filter responses wrt scale.%      h2f  - cell array of bandpass h2 filter responses.%        A  - cell array of monogenic energy responses.%    theta  - cell array of phase orientation responses.%      psi  - cell array of phase angle responses.%% If orientWrap is 1 (on) theta will be returned in the range 0 .. pi and% psi (the phase angle) will be returned in the range -pi .. pi.  If% orientWrap is 0 theta will be returned in the range -pi .. pi and psi will% be returned in the range -pi/2 .. pi/2.  Try both options on an image of a% circle to see what this means!%% Experimentation with sigmaOnf can be useful depending on your application.% I have found values as low as 0.2 (a filter with a *very* large bandwidth)% to be useful on some occasions.%% See also: GABORCONVOLVE% References:% Michael Felsberg and Gerald Sommer. "A New Extension of Linear Signal% Processing for Estimating Local Properties and Detecting Features"% DAGM Symposium 2000, Kiel %% Michael Felsberg and Gerald Sommer. "The monogenic signal" IEEE% Transactions on Signal Processing, 49(12):3136-3144, December 2001% Copyright (c) 2004-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.%  October 2004 - Original version.%  May     2005 - Orientation wrapping and code cleaned up.%  August  2005 - Phase calculation improved.function [f, h1f, h2f, A, theta, psi] = ...	monofilt(im, nscale, minWaveLength, mult, sigmaOnf, orientWrap)    if nargin == 5	orientWrap = 0;  % Default is no orientation wrapping    end        if nargout > 4	thetaPhase = 1;   % Calculate orientation and phase    else	thetaPhase = 0;   % Only return filter outputs    end        [rows,cols] = size(im);        IM = fft2(double(im));        % Generate horizontal and vertical frequency grids that vary from    % -0.5 to 0.5     [u1, u2] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...			([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));    u1 = ifftshift(u1);   % Quadrant shift to put 0 frequency at the corners    u2 = ifftshift(u2);        radius = sqrt(u1.^2 + u2.^2);    % Matrix values contain frequency                                     % values as a radius from centre                                     % (but quadrant shifted)        % Get rid of the 0 radius value in the middle (at top left corner after    % fftshifting) so that taking the log of the radius, or dividing by the    % radius, will not cause trouble.    radius(1,1) = 1;        H1 = i*u1./radius;   % The two monogenic filters in the frequency domain    H2 = i*u2./radius;        % The two monogenic filters H1 and H2 are oriented in frequency space    % but are not selective in terms of the magnitudes of the    % frequencies.  The code below generates bandpass log-Gabor filters    % which are point-wise multiplied by H1 and H2 to produce different    % bandpass versions of H1 and H2        for s = 1:nscale	wavelength = minWaveLength*mult^(s-1);	fo = 1.0/wavelength;                  % Centre frequency of filter.	logGabor = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));  	logGabor(1,1) = 0;                    % undo the radius fudge.			% Generate bandpass versions of H1 and H2 at this scale	H1s = H1.*logGabor; 	H2s = H2.*logGabor; 		%  Apply filters to image in the frequency domain and get spatial        %  results 	f{s} = real(ifft2(IM.*logGabor));    	h1f{s} = real(ifft2(IM.*H1s));	h2f{s} = real(ifft2(IM.*H2s));		A{s} = sqrt(f{s}.^2 + h1f{s}.^2 + h2f{s}.^2);  % Magnitude of Energy.	% If requested calculate the orientation and phase angles	if thetaPhase   	    theta{s} = atan2(h2f{s}, h1f{s});              % Orientation.	    	    % Here phase is measured relative to the h1f-h2f plane as an	    % 'elevation' angle that ranges over +- pi/2	    psi{s} = atan2(f{s}, sqrt(h1f{s}.^2 + h2f{s}.^2));	    	    if orientWrap		% Wrap orientation values back into the range 0-pi		negind = find(theta{s}<0);		theta{s}(negind) = theta{s}(negind) + pi;				% Where orientation values have been wrapped we should                % adjust phase accordingly **check**		psi{s}(negind) = pi-psi{s}(negind);		morethanpi = find(psi{s}>pi);		psi{s}(morethanpi) = psi{s}(morethanpi)-2*pi;	    end	end	    end