% 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