[37] | 1 | % PHASECONG2 - Computes edge and corner phase congruency in an image. |
---|
| 2 | % |
---|
| 3 | % This function calculates the PC_2 measure of phase congruency. |
---|
| 4 | % This function supersedes PHASECONG |
---|
| 5 | % |
---|
| 6 | % There are potentially many arguments, here is the full usage: |
---|
| 7 | % |
---|
| 8 | % [M m or ft pc EO] = phasecong2(im, nscale, norient, minWaveLength, ... |
---|
| 9 | % mult, sigmaOnf, dThetaOnSigma, k, cutOff, g) |
---|
| 10 | % |
---|
| 11 | % However, apart from the image, all parameters have defaults and the |
---|
| 12 | % usage can be as simple as: |
---|
| 13 | % |
---|
| 14 | % M = phasecong2(im); |
---|
| 15 | % |
---|
| 16 | % Arguments: |
---|
| 17 | % Default values Description |
---|
| 18 | % |
---|
| 19 | % nscale 4 - Number of wavelet scales, try values 3-6 |
---|
| 20 | % norient 6 - Number of filter orientations. |
---|
| 21 | % minWaveLength 3 - Wavelength of smallest scale filter. |
---|
| 22 | % mult 2.1 - Scaling factor between successive filters. |
---|
| 23 | % sigmaOnf 0.55 - Ratio of the standard deviation of the Gaussian |
---|
| 24 | % describing the log Gabor filter's transfer function |
---|
| 25 | % in the frequency domain to the filter center frequency. |
---|
| 26 | % dThetaOnSigma 1.2 - Ratio of angular interval between filter orientations |
---|
| 27 | % and the standard deviation of the angular Gaussian |
---|
| 28 | % function used to construct filters in the |
---|
| 29 | % freq. plane. |
---|
| 30 | % k 2.0 - No of standard deviations of the noise energy beyond |
---|
| 31 | % the mean at which we set the noise threshold point. |
---|
| 32 | % You may want to vary this up to a value of 10 or |
---|
| 33 | % 20 for noisy images |
---|
| 34 | % cutOff 0.5 - The fractional measure of frequency spread |
---|
| 35 | % below which phase congruency values get penalized. |
---|
| 36 | % g 10 - Controls the sharpness of the transition in |
---|
| 37 | % the sigmoid function used to weight phase |
---|
| 38 | % congruency for frequency spread. |
---|
| 39 | % |
---|
| 40 | % Returned values: |
---|
| 41 | % M - Maximum moment of phase congruency covariance. |
---|
| 42 | % This is used as a indicator of edge strength. |
---|
| 43 | % m - Minimum moment of phase congruency covariance. |
---|
| 44 | % This is used as a indicator of corner strength. |
---|
| 45 | % or - Orientation image in integer degrees 0-180, |
---|
| 46 | % positive anticlockwise. |
---|
| 47 | % 0 corresponds to a vertical edge, 90 is horizontal. |
---|
| 48 | % ft - *Not correctly implemented at this stage* |
---|
| 49 | % A complex valued image giving the weighted mean |
---|
| 50 | % phase angle at every point in the image for each |
---|
| 51 | % orientation. |
---|
| 52 | % pc - Cell array of phase congruency images (values between 0 and 1) |
---|
| 53 | % for each orientation |
---|
| 54 | % EO - A 2D cell array of complex valued convolution results |
---|
| 55 | % |
---|
| 56 | % EO{s,o} = convolution result for scale s and orientation o. The real part |
---|
| 57 | % is the result of convolving with the even symmetric filter, the imaginary |
---|
| 58 | % part is the result from convolution with the odd symmetric filter. |
---|
| 59 | % |
---|
| 60 | % Hence: |
---|
| 61 | % abs(EO{s,o}) returns the magnitude of the convolution over the |
---|
| 62 | % image at scale s and orientation o. |
---|
| 63 | % angle(EO{s,o}) returns the phase angles. |
---|
| 64 | % |
---|
| 65 | % Notes on specifying parameters: |
---|
| 66 | % |
---|
| 67 | % The parameters can be specified as a full list eg. |
---|
| 68 | % >> [M m or ft pc EO] = phasecong2(im, 5, 6, 3, 2.5, 0.55, 1.2, 2.0, 0.4, 10); |
---|
| 69 | % |
---|
| 70 | % or as a partial list with unspecified parameters taking on default values |
---|
| 71 | % >> [M m or ft pc EO] = phasecong2(im, 5, 6, 3); |
---|
| 72 | % |
---|
| 73 | % or as a partial list of parameters followed by some parameters specified via a |
---|
| 74 | % keyword-value pair, remaining parameters are set to defaults, for example: |
---|
| 75 | % >> [M m or ft pc EO] = phasecong2(im, 5, 6, 3, 'cutOff', 0.3, 'k', 2.5); |
---|
| 76 | % |
---|
| 77 | % The convolutions are done via the FFT. Many of the parameters relate to the |
---|
| 78 | % specification of the filters in the frequency plane. The values do not seem |
---|
| 79 | % to be very critical and the defaults are usually fine. You may want to |
---|
| 80 | % experiment with the values of 'nscales' and 'k', the noise compensation factor. |
---|
| 81 | % |
---|
| 82 | % Notes on filter settings to obtain even coverage of the spectrum |
---|
| 83 | % dthetaOnSigma 1.2 norient 6 |
---|
| 84 | % sigmaOnf .85 mult 1.3 |
---|
| 85 | % sigmaOnf .75 mult 1.6 (filter bandwidth ~1 octave) |
---|
| 86 | % sigmaOnf .65 mult 2.1 |
---|
| 87 | % sigmaOnf .55 mult 3 (filter bandwidth ~2 octaves) |
---|
| 88 | % |
---|
| 89 | % For maximum speed the input image should have dimensions that correspond to |
---|
| 90 | % powers of 2, but the code will operate on images of arbitrary size. |
---|
| 91 | % |
---|
| 92 | % See Also: PHASECONG, PHASESYM, GABORCONVOLVE, PLOTGABORFILTERS |
---|
| 93 | |
---|
| 94 | % References: |
---|
| 95 | % |
---|
| 96 | % Peter Kovesi, "Image Features From Phase Congruency". Videre: A |
---|
| 97 | % Journal of Computer Vision Research. MIT Press. Volume 1, Number 3, |
---|
| 98 | % Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html |
---|
| 99 | % |
---|
| 100 | % Peter Kovesi, "Phase Congruency Detects Corners and |
---|
| 101 | % Edges". Proceedings DICTA 2003, Sydney Dec 10-12 |
---|
| 102 | |
---|
| 103 | % April 1996 Original Version written |
---|
| 104 | % August 1998 Noise compensation corrected. |
---|
| 105 | % October 1998 Noise compensation corrected. - Again!!! |
---|
| 106 | % September 1999 Modified to operate on non-square images of arbitrary size. |
---|
| 107 | % May 2001 Modified to return feature type image. |
---|
| 108 | % July 2003 Altered to calculate 'corner' points. |
---|
| 109 | % October 2003 Speed improvements and refinements. |
---|
| 110 | % July 2005 Better argument handling, changed order of return values |
---|
| 111 | % August 2005 Made Octave compatible |
---|
| 112 | |
---|
| 113 | % Copyright (c) 1996-2005 Peter Kovesi |
---|
| 114 | % School of Computer Science & Software Engineering |
---|
| 115 | % The University of Western Australia |
---|
| 116 | % http://www.csse.uwa.edu.au/ |
---|
| 117 | % |
---|
| 118 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
| 119 | % of this software and associated documentation files (the "Software"), to deal |
---|
| 120 | % in the Software without restriction, subject to the following conditions: |
---|
| 121 | % |
---|
| 122 | % The above copyright notice and this permission notice shall be included in all |
---|
| 123 | % copies or substantial portions of the Software. |
---|
| 124 | % |
---|
| 125 | % The software is provided "as is", without warranty of any kind. |
---|
| 126 | |
---|
| 127 | function [M, m, or, featType, PC, EO]=phasecong2(varargin) |
---|
| 128 | |
---|
| 129 | % Get arguments and/or default values |
---|
| 130 | [im, nscale, norient, minWaveLength, mult, sigmaOnf, ... |
---|
| 131 | dThetaOnSigma,k, cutOff, g] = checkargs(varargin(:)); |
---|
| 132 | |
---|
| 133 | v = version; Octave = v(1)<'5'; % Crude Octave test |
---|
| 134 | epsilon = .0001; % Used to prevent division by zero. |
---|
| 135 | |
---|
| 136 | thetaSigma = pi/norient/dThetaOnSigma; % Calculate the standard deviation of the |
---|
| 137 | % angular Gaussian function used to |
---|
| 138 | % construct filters in the freq. plane. |
---|
| 139 | |
---|
| 140 | [rows,cols] = size(im); |
---|
| 141 | imagefft = fft2(im); % Fourier transform of image |
---|
| 142 | |
---|
| 143 | zero = zeros(rows,cols); |
---|
| 144 | totalEnergy = zero; % Total weighted phase congruency values (energy). |
---|
| 145 | totalSumAn = zero; % Total filter response amplitude values. |
---|
| 146 | orientation = zero; % Matrix storing orientation with greatest |
---|
| 147 | % energy for each pixel. |
---|
| 148 | EO = cell(nscale, norient); % Array of convolution results. |
---|
| 149 | covx2 = zero; % Matrices for covariance data |
---|
| 150 | covy2 = zero; |
---|
| 151 | covxy = zero; |
---|
| 152 | |
---|
| 153 | estMeanE2n = []; |
---|
| 154 | ifftFilterArray = cell(1,nscale); % Array of inverse FFTs of filters |
---|
| 155 | |
---|
| 156 | % Pre-compute some stuff to speed up filter construction |
---|
| 157 | |
---|
| 158 | % Set up X and Y matrices with ranges normalised to +/- 0.5 |
---|
| 159 | % The following code adjusts things appropriately for odd and even values |
---|
| 160 | % of rows and columns. |
---|
| 161 | if mod(cols,2) |
---|
| 162 | xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1); |
---|
| 163 | else |
---|
| 164 | xrange = [-cols/2:(cols/2-1)]/cols; |
---|
| 165 | end |
---|
| 166 | |
---|
| 167 | if mod(rows,2) |
---|
| 168 | yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1); |
---|
| 169 | else |
---|
| 170 | yrange = [-rows/2:(rows/2-1)]/rows; |
---|
| 171 | end |
---|
| 172 | |
---|
| 173 | [x,y] = meshgrid(xrange, yrange); |
---|
| 174 | |
---|
| 175 | radius = sqrt(x.^2 + y.^2); % Matrix values contain *normalised* radius from centre. |
---|
| 176 | radius(rows/2+1, cols/2+1) = 1; % Get rid of the 0 radius value in the middle |
---|
| 177 | % so that taking the log of the radius will |
---|
| 178 | % not cause trouble. |
---|
| 179 | theta = atan2(-y,x); % Matrix values contain polar angle. |
---|
| 180 | % (note -ve y is used to give +ve |
---|
| 181 | % anti-clockwise angles) |
---|
| 182 | radius = ifftshift(radius); % Quadrant shift radius and theta so that filters |
---|
| 183 | theta = ifftshift(theta); % are constructed with 0 frequency at the corners. |
---|
| 184 | |
---|
| 185 | sintheta = sin(theta); |
---|
| 186 | costheta = cos(theta); |
---|
| 187 | clear x; clear y; clear theta; % save a little memory |
---|
| 188 | |
---|
| 189 | % Filters are constructed in terms of two components. |
---|
| 190 | % 1) The radial component, which controls the frequency band that the filter |
---|
| 191 | % responds to |
---|
| 192 | % 2) The angular component, which controls the orientation that the filter |
---|
| 193 | % responds to. |
---|
| 194 | % The two components are multiplied together to construct the overall filter. |
---|
| 195 | |
---|
| 196 | % Construct the radial filter components... |
---|
| 197 | |
---|
| 198 | % First construct a low-pass filter that is as large as possible, yet falls |
---|
| 199 | % away to zero at the boundaries. All log Gabor filters are multiplied by |
---|
| 200 | % this to ensure no extra frequencies at the 'corners' of the FFT are |
---|
| 201 | % incorporated as this seems to upset the normalisation process when |
---|
| 202 | % calculating phase congrunecy. |
---|
| 203 | lp = lowpassfilter([rows,cols],.45,15); % Radius .45, 'sharpness' 15 |
---|
| 204 | |
---|
| 205 | logGabor = cell(1,nscale); |
---|
| 206 | |
---|
| 207 | for s = 1:nscale |
---|
| 208 | wavelength = minWaveLength*mult^(s-1); |
---|
| 209 | fo = 1.0/wavelength; % Centre frequency of filter. |
---|
| 210 | logGabor{s} = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2)); |
---|
| 211 | logGabor{s} = logGabor{s}.*lp; % Apply low-pass filter |
---|
| 212 | logGabor{s}(1,1) = 0; % Set the value at the 0 frequency point of the filter |
---|
| 213 | % back to zero (undo the radius fudge). |
---|
| 214 | end |
---|
| 215 | |
---|
| 216 | % Then construct the angular filter components... |
---|
| 217 | |
---|
| 218 | spread = cell(1,norient); |
---|
| 219 | |
---|
| 220 | for o = 1:norient |
---|
| 221 | angl = (o-1)*pi/norient; % Filter angle. |
---|
| 222 | |
---|
| 223 | % For each point in the filter matrix calculate the angular distance from |
---|
| 224 | % the specified filter orientation. To overcome the angular wrap-around |
---|
| 225 | % problem sine difference and cosine difference values are first computed |
---|
| 226 | % and then the atan2 function is used to determine angular distance. |
---|
| 227 | |
---|
| 228 | ds = sintheta * cos(angl) - costheta * sin(angl); % Difference in sine. |
---|
| 229 | dc = costheta * cos(angl) + sintheta * sin(angl); % Difference in cosine. |
---|
| 230 | dtheta = abs(atan2(ds,dc)); % Absolute angular distance. |
---|
| 231 | spread{o} = exp((-dtheta.^2) / (2 * thetaSigma^2)); % Calculate the |
---|
| 232 | % angular filter component. |
---|
| 233 | end |
---|
| 234 | |
---|
| 235 | % The main loop... |
---|
| 236 | |
---|
| 237 | for o = 1:norient % For each orientation. |
---|
| 238 | fprintf('Processing orientation %d\r',o); |
---|
| 239 | if Octave fflush(1); end |
---|
| 240 | |
---|
| 241 | angl = (o-1)*pi/norient; % Filter angle. |
---|
| 242 | sumE_ThisOrient = zero; % Initialize accumulator matrices. |
---|
| 243 | sumO_ThisOrient = zero; |
---|
| 244 | sumAn_ThisOrient = zero; |
---|
| 245 | Energy = zero; |
---|
| 246 | |
---|
| 247 | for s = 1:nscale, % For each scale. |
---|
| 248 | filter = logGabor{s} .* spread{o}; % Multiply radial and angular |
---|
| 249 | % components to get the filter. |
---|
| 250 | |
---|
| 251 | % if o == 1 % accumulate filter info for noise compensation (nominally the same |
---|
| 252 | % for all orientations, hence it is only done once) |
---|
| 253 | ifftFilt = real(ifft2(filter))*sqrt(rows*cols); % Note rescaling to match power |
---|
| 254 | ifftFilterArray{s} = ifftFilt; % record ifft2 of filter |
---|
| 255 | % end |
---|
| 256 | |
---|
| 257 | % Convolve image with even and odd filters returning the result in EO |
---|
| 258 | EO{s,o} = ifft2(imagefft .* filter); |
---|
| 259 | |
---|
| 260 | An = abs(EO{s,o}); % Amplitude of even & odd filter response. |
---|
| 261 | sumAn_ThisOrient = sumAn_ThisOrient + An; % Sum of amplitude responses. |
---|
| 262 | sumE_ThisOrient = sumE_ThisOrient + real(EO{s,o}); % Sum of even filter convolution results. |
---|
| 263 | sumO_ThisOrient = sumO_ThisOrient + imag(EO{s,o}); % Sum of odd filter convolution results. |
---|
| 264 | |
---|
| 265 | if s==1 % Record mean squared filter value at smallest |
---|
| 266 | EM_n = sum(sum(filter.^2)); % scale. This is used for noise estimation. |
---|
| 267 | maxAn = An; % Record the maximum An over all scales. |
---|
| 268 | else |
---|
| 269 | maxAn = max(maxAn, An); |
---|
| 270 | end |
---|
| 271 | |
---|
| 272 | end % ... and process the next scale |
---|
| 273 | |
---|
| 274 | % Get weighted mean filter response vector, this gives the weighted mean |
---|
| 275 | % phase angle. |
---|
| 276 | |
---|
| 277 | XEnergy = sqrt(sumE_ThisOrient.^2 + sumO_ThisOrient.^2) + epsilon; |
---|
| 278 | MeanE = sumE_ThisOrient ./ XEnergy; |
---|
| 279 | MeanO = sumO_ThisOrient ./ XEnergy; |
---|
| 280 | |
---|
| 281 | % Now calculate An(cos(phase_deviation) - | sin(phase_deviation)) | by |
---|
| 282 | % using dot and cross products between the weighted mean filter response |
---|
| 283 | % vector and the individual filter response vectors at each scale. This |
---|
| 284 | % quantity is phase congruency multiplied by An, which we call energy. |
---|
| 285 | |
---|
| 286 | for s = 1:nscale, |
---|
| 287 | E = real(EO{s,o}); O = imag(EO{s,o}); % Extract even and odd |
---|
| 288 | % convolution results. |
---|
| 289 | Energy = Energy + E.*MeanE + O.*MeanO - abs(E.*MeanO - O.*MeanE); |
---|
| 290 | end |
---|
| 291 | |
---|
| 292 | % Compensate for noise |
---|
| 293 | % We estimate the noise power from the energy squared response at the |
---|
| 294 | % smallest scale. If the noise is Gaussian the energy squared will have a |
---|
| 295 | % Chi-squared 2DOF pdf. We calculate the median energy squared response |
---|
| 296 | % as this is a robust statistic. From this we estimate the mean. |
---|
| 297 | % The estimate of noise power is obtained by dividing the mean squared |
---|
| 298 | % energy value by the mean squared filter value |
---|
| 299 | |
---|
| 300 | medianE2n = median(reshape(abs(EO{1,o}).^2,1,rows*cols)); |
---|
| 301 | meanE2n = -medianE2n/log(0.5); |
---|
| 302 | estMeanE2n(o) = meanE2n; |
---|
| 303 | |
---|
| 304 | noisePower = meanE2n/EM_n; % Estimate of noise power. |
---|
| 305 | |
---|
| 306 | % if o == 1 |
---|
| 307 | % Now estimate the total energy^2 due to noise |
---|
| 308 | % Estimate for sum(An^2) + sum(Ai.*Aj.*(cphi.*cphj + sphi.*sphj)) |
---|
| 309 | |
---|
| 310 | EstSumAn2 = zero; |
---|
| 311 | for s = 1:nscale |
---|
| 312 | EstSumAn2 = EstSumAn2 + ifftFilterArray{s}.^2; |
---|
| 313 | end |
---|
| 314 | |
---|
| 315 | EstSumAiAj = zero; |
---|
| 316 | for si = 1:(nscale-1) |
---|
| 317 | for sj = (si+1):nscale |
---|
| 318 | EstSumAiAj = EstSumAiAj + ifftFilterArray{si}.*ifftFilterArray{sj}; |
---|
| 319 | end |
---|
| 320 | end |
---|
| 321 | sumEstSumAn2 = sum(sum(EstSumAn2)); |
---|
| 322 | sumEstSumAiAj = sum(sum(EstSumAiAj)); |
---|
| 323 | |
---|
| 324 | % end % if o == 1 |
---|
| 325 | |
---|
| 326 | EstNoiseEnergy2 = 2*noisePower*sumEstSumAn2 + 4*noisePower*sumEstSumAiAj; |
---|
| 327 | |
---|
| 328 | tau = sqrt(EstNoiseEnergy2/2); % Rayleigh parameter |
---|
| 329 | EstNoiseEnergy = tau*sqrt(pi/2); % Expected value of noise energy |
---|
| 330 | EstNoiseEnergySigma = sqrt( (2-pi/2)*tau^2 ); |
---|
| 331 | |
---|
| 332 | T = EstNoiseEnergy + k*EstNoiseEnergySigma; % Noise threshold |
---|
| 333 | |
---|
| 334 | % The estimated noise effect calculated above is only valid for the PC_1 measure. |
---|
| 335 | % The PC_2 measure does not lend itself readily to the same analysis. However |
---|
| 336 | % empirically it seems that the noise effect is overestimated roughly by a factor |
---|
| 337 | % of 1.7 for the filter parameters used here. |
---|
| 338 | |
---|
| 339 | T = T/1.7; % Empirical rescaling of the estimated noise effect to |
---|
| 340 | % suit the PC_2 phase congruency measure |
---|
| 341 | |
---|
| 342 | Energy = max(Energy - T, zero); % Apply noise threshold |
---|
| 343 | |
---|
| 344 | % Form weighting that penalizes frequency distributions that are |
---|
| 345 | % particularly narrow. Calculate fractional 'width' of the frequencies |
---|
| 346 | % present by taking the sum of the filter response amplitudes and dividing |
---|
| 347 | % by the maximum amplitude at each point on the image. |
---|
| 348 | |
---|
| 349 | width = sumAn_ThisOrient ./ (maxAn + epsilon) / nscale; |
---|
| 350 | |
---|
| 351 | % Now calculate the sigmoidal weighting function for this orientation. |
---|
| 352 | |
---|
| 353 | weight = 1.0 ./ (1 + exp( (cutOff - width)*g)); |
---|
| 354 | |
---|
| 355 | % Apply weighting to energy and then calculate phase congruency |
---|
| 356 | |
---|
| 357 | PC{o} = weight.*Energy./sumAn_ThisOrient; % Phase congruency for this orientation |
---|
| 358 | featType{o} = E+i*O; |
---|
| 359 | |
---|
| 360 | % Build up covariance data for every point |
---|
| 361 | covx = PC{o}*cos(angl); |
---|
| 362 | covy = PC{o}*sin(angl); |
---|
| 363 | covx2 = covx2 + covx.^2; |
---|
| 364 | covy2 = covy2 + covy.^2; |
---|
| 365 | covxy = covxy + covx.*covy; |
---|
| 366 | |
---|
| 367 | end % For each orientation |
---|
| 368 | |
---|
| 369 | fprintf(' \r'); |
---|
| 370 | |
---|
| 371 | % Edge and Corner calculations |
---|
| 372 | |
---|
| 373 | % The following is optimised code to calculate principal vector |
---|
| 374 | % of the phase congruency covariance data and to calculate |
---|
| 375 | % the minimumum and maximum moments - these correspond to |
---|
| 376 | % the singular values. |
---|
| 377 | |
---|
| 378 | % First normalise covariance values by the number of orientations/2 |
---|
| 379 | |
---|
| 380 | covx2 = covx2/(norient/2); |
---|
| 381 | covy2 = covy2/(norient/2); |
---|
| 382 | covxy = covxy/norient; % This gives us 2*covxy/(norient/2) |
---|
| 383 | |
---|
| 384 | denom = sqrt(covxy.^2 + (covx2-covy2).^2)+epsilon; |
---|
| 385 | sin2theta = covxy./denom; |
---|
| 386 | cos2theta = (covx2-covy2)./denom; |
---|
| 387 | or = atan2(sin2theta,cos2theta)/2; % Orientation perpendicular to edge. |
---|
| 388 | or = round(or*180/pi); % Return result rounded to integer |
---|
| 389 | % degrees. |
---|
| 390 | neg = or < 0; |
---|
| 391 | or = ~neg.*or + neg.*(or+180); % Adjust range from -90 to 90 |
---|
| 392 | % to 0 to 180. |
---|
| 393 | |
---|
| 394 | M = (covy2+covx2 + denom)/2; % Maximum moment |
---|
| 395 | m = (covy2+covx2 - denom)/2; % ... and minimum moment |
---|
| 396 | |
---|
| 397 | |
---|
| 398 | |
---|
| 399 | |
---|
| 400 | |
---|
| 401 | %------------------------------------------------------------------ |
---|
| 402 | % CHECKARGS |
---|
| 403 | % |
---|
| 404 | % Function to process the arguments that have been supplied, assign |
---|
| 405 | % default values as needed and perform basic checks. |
---|
| 406 | |
---|
| 407 | function [im, nscale, norient, minWaveLength, mult, sigmaOnf, ... |
---|
| 408 | dThetaOnSigma,k, cutOff, g] = checkargs(arg); |
---|
| 409 | |
---|
| 410 | nargs = length(arg); |
---|
| 411 | |
---|
| 412 | if nargs < 1 |
---|
| 413 | error('No image supplied as an argument'); |
---|
| 414 | end |
---|
| 415 | |
---|
| 416 | % Set up default values for all arguments and then overwrite them |
---|
| 417 | % with with any new values that may be supplied |
---|
| 418 | im = []; |
---|
| 419 | nscale = 4; % Number of wavelet scales. |
---|
| 420 | norient = 6; % Number of filter orientations. |
---|
| 421 | minWaveLength = 3; % Wavelength of smallest scale filter. |
---|
| 422 | mult = 2.1; % Scaling factor between successive filters. |
---|
| 423 | sigmaOnf = 0.55; % Ratio of the standard deviation of the |
---|
| 424 | % Gaussian describing the log Gabor filter's |
---|
| 425 | % transfer function in the frequency domain |
---|
| 426 | % to the filter center frequency. |
---|
| 427 | dThetaOnSigma = 1.2; % Ratio of angular interval between filter orientations |
---|
| 428 | % and the standard deviation of the angular Gaussian |
---|
| 429 | % function used to construct filters in the |
---|
| 430 | % freq. plane. |
---|
| 431 | k = 2.0; % No of standard deviations of the noise |
---|
| 432 | % energy beyond the mean at which we set the |
---|
| 433 | % noise threshold point. |
---|
| 434 | cutOff = 0.5; % The fractional measure of frequency spread |
---|
| 435 | % below which phase congruency values get penalized. |
---|
| 436 | g = 10; % Controls the sharpness of the transition in |
---|
| 437 | % the sigmoid function used to weight phase |
---|
| 438 | % congruency for frequency spread. |
---|
| 439 | |
---|
| 440 | % Allowed argument reading states |
---|
| 441 | allnumeric = 1; % Numeric argument values in predefined order |
---|
| 442 | keywordvalue = 2; % Arguments in the form of string keyword |
---|
| 443 | % followed by numeric value |
---|
| 444 | readstate = allnumeric; % Start in the allnumeric state |
---|
| 445 | |
---|
| 446 | if readstate == allnumeric |
---|
| 447 | for n = 1:nargs |
---|
| 448 | if isa(arg{n},'char') |
---|
| 449 | readstate = keywordvalue; |
---|
| 450 | break; |
---|
| 451 | else |
---|
| 452 | if n == 1, im = arg{n}; |
---|
| 453 | elseif n == 2, nscale = arg{n}; |
---|
| 454 | elseif n == 3, norient = arg{n}; |
---|
| 455 | elseif n == 4, minWaveLength = arg{n}; |
---|
| 456 | elseif n == 5, mult = arg{n}; |
---|
| 457 | elseif n == 6, sigmaOnf = arg{n}; |
---|
| 458 | elseif n == 7, dThetaOnSigma = arg{n}; |
---|
| 459 | elseif n == 8, k = arg{n}; |
---|
| 460 | elseif n == 9, cutOff = arg{n}; |
---|
| 461 | elseif n == 10,g = arg{n}; |
---|
| 462 | end |
---|
| 463 | end |
---|
| 464 | end |
---|
| 465 | end |
---|
| 466 | |
---|
| 467 | % Code to handle parameter name - value pairs |
---|
| 468 | if readstate == keywordvalue |
---|
| 469 | while n < nargs |
---|
| 470 | |
---|
| 471 | if ~isa(arg{n},'char') | ~isa(arg{n+1}, 'double') |
---|
| 472 | error('There should be a parameter name - value pair'); |
---|
| 473 | end |
---|
| 474 | |
---|
| 475 | if strncmpi(arg{n},'im' ,2), im = arg{n+1}; |
---|
| 476 | elseif strncmpi(arg{n},'nscale' ,2), nscale = arg{n+1}; |
---|
| 477 | elseif strncmpi(arg{n},'norient' ,2), norient = arg{n+1}; |
---|
| 478 | elseif strncmpi(arg{n},'minWaveLength',2), minWavelength = arg{n+1}; |
---|
| 479 | elseif strncmpi(arg{n},'mult' ,2), mult = arg{n+1}; |
---|
| 480 | elseif strncmpi(arg{n},'sigmaOnf',2), sigmaOnf = arg{n+1}; |
---|
| 481 | elseif strncmpi(arg{n},'dthetaOnSigma',2), dThetaOnSigma = arg{n+1}; |
---|
| 482 | elseif strncmpi(arg{n},'k' ,1), k = arg{n+1}; |
---|
| 483 | elseif strncmpi(arg{n},'cutOff' ,2), cutOff = arg{n+1}; |
---|
| 484 | elseif strncmpi(arg{n},'g' ,1), g = arg{n+1}; |
---|
| 485 | else error('Unrecognised parameter name'); |
---|
| 486 | end |
---|
| 487 | |
---|
| 488 | n = n+2; |
---|
| 489 | if n == nargs |
---|
| 490 | error('Unmatched parameter name - value pair'); |
---|
| 491 | end |
---|
| 492 | |
---|
| 493 | end |
---|
| 494 | end |
---|
| 495 | |
---|
| 496 | if isempty(im) |
---|
| 497 | error('No image argument supplied'); |
---|
| 498 | end |
---|
| 499 | |
---|
| 500 | if ~isa(im, 'double') |
---|
| 501 | im = double(im); |
---|
| 502 | end |
---|
| 503 | |
---|
| 504 | if nscale < 1 |
---|
| 505 | error('nscale must be an integer >= 1'); |
---|
| 506 | end |
---|
| 507 | |
---|
| 508 | if norient < 1 |
---|
| 509 | error('norient must be an integer >= 1'); |
---|
| 510 | end |
---|
| 511 | |
---|
| 512 | if minWaveLength < 2 |
---|
| 513 | error('It makes little sense to have a wavelength < 2'); |
---|
| 514 | end |
---|
| 515 | |
---|
| 516 | if cutOff < 0 | cutOff > 1 |
---|
| 517 | error('Cut off value must be between 0 and 1'); |
---|
| 518 | end |
---|
| 519 | |
---|
| 520 | |
---|
| 521 | |
---|
| 522 | |
---|
| 523 | |
---|
| 524 | |
---|