[37] | 1 | % MONOFILT - Apply monogenic filters to an image to obtain 2D analytic signal |
---|
| 2 | % |
---|
| 3 | % Implementation of Felsberg's monogenic filters |
---|
| 4 | % |
---|
| 5 | % Usage: [f, h1f, h2f, A, theta, psi] = ... |
---|
| 6 | % monofilt(im, nscale, minWaveLength, mult, sigmaOnf, orientWrap) |
---|
| 7 | % 3 4 2 0.65 1/0 |
---|
| 8 | % Arguments: |
---|
| 9 | % The convolutions are done via the FFT. Many of the parameters relate |
---|
| 10 | % to the specification of the filters in the frequency plane. |
---|
| 11 | % |
---|
| 12 | % Variable Suggested Description |
---|
| 13 | % name value |
---|
| 14 | % ---------------------------------------------------------- |
---|
| 15 | % im Image to be convolved. |
---|
| 16 | % nscale = 3; Number of filter scales. |
---|
| 17 | % minWaveLength = 4; Wavelength of smallest scale filter. |
---|
| 18 | % mult = 2; Scaling factor between successive filters. |
---|
| 19 | % sigmaOnf = 0.65; Ratio of the standard deviation of the |
---|
| 20 | % Gaussian describing the log Gabor filter's |
---|
| 21 | % transfer function in the frequency domain |
---|
| 22 | % to the filter center frequency. |
---|
| 23 | % orientWrap 1/0 Optional flag 1/0 to turn on/off |
---|
| 24 | % 'wrapping' of orientation data from a |
---|
| 25 | % range of -pi .. pi to the range 0 .. pi. |
---|
| 26 | % This affects the interpretation of the |
---|
| 27 | % phase angle - see note below. Defaults to 0. |
---|
| 28 | % Returns: |
---|
| 29 | % |
---|
| 30 | % f - cell array of bandpass filter responses with respect to scale. |
---|
| 31 | % h1f - cell array of bandpass h1 filter responses wrt scale. |
---|
| 32 | % h2f - cell array of bandpass h2 filter responses. |
---|
| 33 | % A - cell array of monogenic energy responses. |
---|
| 34 | % theta - cell array of phase orientation responses. |
---|
| 35 | % psi - cell array of phase angle responses. |
---|
| 36 | % |
---|
| 37 | % If orientWrap is 1 (on) theta will be returned in the range 0 .. pi and |
---|
| 38 | % psi (the phase angle) will be returned in the range -pi .. pi. If |
---|
| 39 | % orientWrap is 0 theta will be returned in the range -pi .. pi and psi will |
---|
| 40 | % be returned in the range -pi/2 .. pi/2. Try both options on an image of a |
---|
| 41 | % circle to see what this means! |
---|
| 42 | % |
---|
| 43 | % Experimentation with sigmaOnf can be useful depending on your application. |
---|
| 44 | % I have found values as low as 0.2 (a filter with a *very* large bandwidth) |
---|
| 45 | % to be useful on some occasions. |
---|
| 46 | % |
---|
| 47 | % See also: GABORCONVOLVE |
---|
| 48 | |
---|
| 49 | % References: |
---|
| 50 | % Michael Felsberg and Gerald Sommer. "A New Extension of Linear Signal |
---|
| 51 | % Processing for Estimating Local Properties and Detecting Features" |
---|
| 52 | % DAGM Symposium 2000, Kiel |
---|
| 53 | % |
---|
| 54 | % Michael Felsberg and Gerald Sommer. "The monogenic signal" IEEE |
---|
| 55 | % Transactions on Signal Processing, 49(12):3136-3144, December 2001 |
---|
| 56 | |
---|
| 57 | % Copyright (c) 2004-2005 Peter Kovesi |
---|
| 58 | % School of Computer Science & Software Engineering |
---|
| 59 | % The University of Western Australia |
---|
| 60 | % http://www.csse.uwa.edu.au/ |
---|
| 61 | % |
---|
| 62 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
| 63 | % of this software and associated documentation files (the "Software"), to deal |
---|
| 64 | % in the Software without restriction, subject to the following conditions: |
---|
| 65 | % |
---|
| 66 | % The above copyright notice and this permission notice shall be included in |
---|
| 67 | % all copies or substantial portions of the Software. |
---|
| 68 | % |
---|
| 69 | % The Software is provided "as is", without warranty of any kind. |
---|
| 70 | |
---|
| 71 | % October 2004 - Original version. |
---|
| 72 | % May 2005 - Orientation wrapping and code cleaned up. |
---|
| 73 | % August 2005 - Phase calculation improved. |
---|
| 74 | |
---|
| 75 | function [f, h1f, h2f, A, theta, psi] = ... |
---|
| 76 | monofilt(im, nscale, minWaveLength, mult, sigmaOnf, orientWrap) |
---|
| 77 | |
---|
| 78 | if nargin == 5 |
---|
| 79 | orientWrap = 0; % Default is no orientation wrapping |
---|
| 80 | end |
---|
| 81 | |
---|
| 82 | if nargout > 4 |
---|
| 83 | thetaPhase = 1; % Calculate orientation and phase |
---|
| 84 | else |
---|
| 85 | thetaPhase = 0; % Only return filter outputs |
---|
| 86 | end |
---|
| 87 | |
---|
| 88 | [rows,cols] = size(im); |
---|
| 89 | IM = fft2(double(im)); |
---|
| 90 | |
---|
| 91 | % Generate horizontal and vertical frequency grids that vary from |
---|
| 92 | % -0.5 to 0.5 |
---|
| 93 | [u1, u2] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ... |
---|
| 94 | ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2))); |
---|
| 95 | |
---|
| 96 | u1 = ifftshift(u1); % Quadrant shift to put 0 frequency at the corners |
---|
| 97 | u2 = ifftshift(u2); |
---|
| 98 | |
---|
| 99 | radius = sqrt(u1.^2 + u2.^2); % Matrix values contain frequency |
---|
| 100 | % values as a radius from centre |
---|
| 101 | % (but quadrant shifted) |
---|
| 102 | |
---|
| 103 | % Get rid of the 0 radius value in the middle (at top left corner after |
---|
| 104 | % fftshifting) so that taking the log of the radius, or dividing by the |
---|
| 105 | % radius, will not cause trouble. |
---|
| 106 | radius(1,1) = 1; |
---|
| 107 | |
---|
| 108 | H1 = i*u1./radius; % The two monogenic filters in the frequency domain |
---|
| 109 | H2 = i*u2./radius; |
---|
| 110 | |
---|
| 111 | % The two monogenic filters H1 and H2 are oriented in frequency space |
---|
| 112 | % but are not selective in terms of the magnitudes of the |
---|
| 113 | % frequencies. The code below generates bandpass log-Gabor filters |
---|
| 114 | % which are point-wise multiplied by H1 and H2 to produce different |
---|
| 115 | % bandpass versions of H1 and H2 |
---|
| 116 | |
---|
| 117 | for s = 1:nscale |
---|
| 118 | wavelength = minWaveLength*mult^(s-1); |
---|
| 119 | fo = 1.0/wavelength; % Centre frequency of filter. |
---|
| 120 | logGabor = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2)); |
---|
| 121 | logGabor(1,1) = 0; % undo the radius fudge. |
---|
| 122 | |
---|
| 123 | % Generate bandpass versions of H1 and H2 at this scale |
---|
| 124 | H1s = H1.*logGabor; |
---|
| 125 | H2s = H2.*logGabor; |
---|
| 126 | |
---|
| 127 | % Apply filters to image in the frequency domain and get spatial |
---|
| 128 | % results |
---|
| 129 | f{s} = real(ifft2(IM.*logGabor)); |
---|
| 130 | h1f{s} = real(ifft2(IM.*H1s)); |
---|
| 131 | h2f{s} = real(ifft2(IM.*H2s)); |
---|
| 132 | |
---|
| 133 | A{s} = sqrt(f{s}.^2 + h1f{s}.^2 + h2f{s}.^2); % Magnitude of Energy. |
---|
| 134 | |
---|
| 135 | % If requested calculate the orientation and phase angles |
---|
| 136 | if thetaPhase |
---|
| 137 | |
---|
| 138 | theta{s} = atan2(h2f{s}, h1f{s}); % Orientation. |
---|
| 139 | |
---|
| 140 | % Here phase is measured relative to the h1f-h2f plane as an |
---|
| 141 | % 'elevation' angle that ranges over +- pi/2 |
---|
| 142 | psi{s} = atan2(f{s}, sqrt(h1f{s}.^2 + h2f{s}.^2)); |
---|
| 143 | |
---|
| 144 | if orientWrap |
---|
| 145 | % Wrap orientation values back into the range 0-pi |
---|
| 146 | negind = find(theta{s}<0); |
---|
| 147 | theta{s}(negind) = theta{s}(negind) + pi; |
---|
| 148 | |
---|
| 149 | % Where orientation values have been wrapped we should |
---|
| 150 | % adjust phase accordingly **check** |
---|
| 151 | psi{s}(negind) = pi-psi{s}(negind); |
---|
| 152 | morethanpi = find(psi{s}>pi); |
---|
| 153 | psi{s}(morethanpi) = psi{s}(morethanpi)-2*pi; |
---|
| 154 | end |
---|
| 155 | end |
---|
| 156 | |
---|
| 157 | end |
---|
| 158 | |
---|
| 159 | |
---|