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