-
Notifications
You must be signed in to change notification settings - Fork 6
/
phasesym.m
459 lines (399 loc) · 21.1 KB
/
phasesym.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
% PHASESYM - Function for computing phase symmetry on an image.
%
% This function calculates the phase symmetry of points in an image.
% This is a contrast invariant measure of symmetry. This function can be
% used as a line and blob detector. The greyscale 'polarity' of the lines
% that you want to find can be specified.
%
% There are potentially many arguments, here is the full usage:
%
% [phaseSym, orientation, totalEnergy, T] = ...
% phasesym(im, nscale, norient, minWaveLength, mult, ...
% sigmaOnf, k, polarity, noiseMethod)
%
% However, apart from the image, all parameters have defaults and the
% usage can be as simple as:
%
% phaseSym = phasesym(im);
%
% Arguments:
% Default values Description
%
% nscale 5 - Number of wavelet scales, try values 3-6
% norient 6 - Number of filter orientations.
% minWaveLength 3 - Wavelength of smallest scale filter.
% mult 2.1 - Scaling factor between successive filters.
% sigmaOnf 0.55 - 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.
% k 2.0 - No of standard deviations of the noise energy beyond
% the mean at which we set the noise threshold point.
% You may want to vary this up to a value of 10 or
% 20 for noisy images
% polarity 0 - Controls 'polarity' of symmetry features to find.
% 1 - just return 'bright' points
% -1 - just return 'dark' points
% 0 - return bright and dark points.
% noiseMethod -1 - Parameter specifies method used to determine
% noise statistics.
% -1 use median of smallest scale filter responses
% -2 use mode of smallest scale filter responses
% 0+ use noiseMethod value as the fixed noise threshold.
%
% Return values:
% phaseSym - Phase symmetry image (values between 0 and 1).
% orientation - Orientation image. Orientation in which local
% symmetry energy is a maximum, in degrees
% (0-180), angles positive anti-clockwise. Note
% the orientation info is quantized by the number
% of orientations
% totalEnergy - Un-normalised raw symmetry energy which may be
% more to your liking.
% T - Calculated noise threshold (can be useful for
% diagnosing noise characteristics of images). Once you know
% this you can then specify fixed thresholds and save some
% computation time.
%
% Notes on specifying parameters:
%
% The parameters can be specified as a full list eg.
% >> phaseSym = phasesym(im, 5, 6, 3, 2.5, 0.55, 2.0, 0);
%
% or as a partial list with unspecified parameters taking on default values
% >> phaseSym = phasesym(im, 5, 6, 3);
%
% or as a partial list of parameters followed by some parameters specified via a
% keyword-value pair, remaining parameters are set to defaults, for example:
% >> phaseSym = phasesym(im, 5, 6, 3, 'polarity',-1, 'k', 2.5);
%
% The convolutions are done via the FFT. Many of the parameters relate to the
% specification of the filters in the frequency plane. The values do not seem
% to be very critical and the defaults are usually fine. You may want to
% experiment with the values of 'nscales' and 'k', the noise compensation factor.
%
% Notes on filter settings to obtain even coverage of the spectrum
% sigmaOnf .85 mult 1.3
% sigmaOnf .75 mult 1.6 (filter bandwidth ~1 octave)
% sigmaOnf .65 mult 2.1
% sigmaOnf .55 mult 3 (filter bandwidth ~2 octaves)
%
% For maximum speed the input image should have dimensions that correspond to
% powers of 2, but the code will operate on images of arbitrary size.
%
% See Also: PHASECONG, PHASECONG2, GABORCONVOLVE, PLOTGABORFILTERS
% References:
% Peter Kovesi, "Symmetry and Asymmetry From Local Phase" AI'97, Tenth
% Australian Joint Conference on Artificial Intelligence. 2 - 4 December
% 1997. http://www.cs.uwa.edu.au/pub/robvis/papers/pk/ai97.ps.gz.
%
% Peter Kovesi, "Image Features From Phase Congruency". Videre: A
% Journal of Computer Vision Research. MIT Press. Volume 1, Number 3,
% Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html
% April 1996 Original Version written
% August 1998 Noise compensation corrected.
% October 1998 Noise compensation corrected. - Again!!!
% September 1999 Modified to operate on non-square images of arbitrary size.
% February 2001 Specialised from phasecong.m to calculate phase symmetry
% July 2005 Better argument handling + general cleanup and speed improvements
% August 2005 Made Octave compatible.
% January 2007 Small correction and cleanup of radius calculation for odd
% image sizes.
% May 2009 Noise compensation simplified reducing memory and
% computation overhead. Spread function changed to a cosine
% eliminating parameter dThetaOnSigma and ensuring even
% angular coverage.
% September 2017 Changed to use FILTERGRID
% Copyright (c) 1996-2017 Peter Kovesi
% http://www.peterkovesi.com
%
% 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.
function[phaseSym, orientation, totalEnergy, T] = phasesym(varargin)
% Get arguments and/or default values
[im, nscale, norient, minWaveLength, mult, sigmaOnf, k, ...
polarity, noiseMethod] = checkargs(varargin(:));
epsilon = 1e-4; % Used to prevent division by zero.
[rows,cols] = size(im);
imagefft = fft2(im); % Fourier transform of image
zero = zeros(rows,cols);
totalEnergy = zero; % Matrix for accumulating weighted phase
% congruency values (energy).
totalSumAn = zero; % Matrix for accumulating filter response
% amplitude values.
orientation = zero; % Matrix storing orientation with greatest
% energy for each pixel.
% Pre-compute data for filter construction
[radius, x, y] = FilterGrid(rows,cols);
theta = atan2(-y,x); % Matrix values contain polar angle.
% (note -ve y is used to give +ve
% anti-clockwise angles)
radius(1,1) = 1; % Get rid of the 0 radius value at the 0
% frequency point (now at top-left corner)
% so that taking the log of the radius will
% not cause trouble.
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta; % save a little memory
% Filters are constructed in terms of two components.
% 1) The radial component, which controls the frequency band that the filter
% responds to
% 2) The angular component, which controls the orientation that the filter
% responds to.
% The two components are multiplied together to construct the overall filter.
% Construct the radial filter components...
% First construct a low-pass filter that is as large as possible, yet falls
% away to zero at the boundaries. All log Gabor filters are multiplied by
% this to ensure no extra frequencies at the 'corners' of the FFT are
% incorporated as this seems to upset the normalisation process when
% calculating phase congrunecy.
lp = LowPassFilter([rows,cols],.4,10); % Radius .4, 'sharpness' 10
logGabor = cell(1,nscale);
for s = 1:nscale
wavelength = minWaveLength*mult^(s-1);
fo = 1.0/wavelength; % Centre frequency of filter.
logGabor{s} = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));
logGabor{s} = logGabor{s}.*lp; % Apply low-pass filter
logGabor{s}(1,1) = 0; % Set the value at the 0 frequency point of the filter
% back to zero (undo the radius fudge).
end
%% The main loop...
for o = 1:norient % For each orientation....
% Construct the angular filter spread function
angl = (o-1)*pi/norient; % Filter angle.
% For each point in the filter matrix calculate the angular distance from
% the specified filter orientation. To overcome the angular wrap-around
% problem sine difference and cosine difference values are first computed
% and then the atan2 function is used to determine angular distance.
ds = sintheta * cos(angl) - costheta * sin(angl); % Difference in sine.
dc = costheta * cos(angl) + sintheta * sin(angl); % Difference in cosine.
dtheta = abs(atan2(ds,dc)); % Absolute angular distance.
% Scale theta so that cosine spread function has the right wavelength
% and clamp to pi.
dtheta = min(dtheta*norient/2,pi);
% The spread function is cos(dtheta) between -pi and pi. We add 1,
% and then divide by 2 so that the value ranges 0-1
spread = (cos(dtheta)+1)/2;
sumAn_ThisOrient = zero;
Energy_ThisOrient = zero;
for s = 1:nscale % For each scale....
filter = logGabor{s} .* spread; % Multiply radial and angular
% components to get filter.
% Convolve image with even and odd filters returning the result in EO
EO = ifft2(imagefft .* filter);
An = abs(EO); % Amplitude of even & odd filter response.
sumAn_ThisOrient = sumAn_ThisOrient + An; % Sum of amplitude responses.
% At the smallest scale estimate noise characteristics from the
% distribution of the filter amplitude responses stored in sumAn.
% tau is the Rayleigh parameter that is used to describe the
% distribution.
if s == 1
if noiseMethod == -1 % Use median to estimate noise statistics
tau = median(sumAn_ThisOrient(:))/sqrt(log(4));
elseif noiseMethod == -2 % Use mode to estimate noise statistics
tau = rayleighmode(sumAn_ThisOrient(:));
end
end
% Now calculate the phase symmetry measure.
if polarity == 0 % look for 'white' and 'black' spots
Energy_ThisOrient = Energy_ThisOrient ...
+ abs(real(EO)) - abs(imag(EO));
elseif polarity == 1 % Just look for 'white' spots
Energy_ThisOrient = Energy_ThisOrient ...
+ real(EO) - abs(imag(EO));
elseif polarity == -1 % Just look for 'black' spots
Energy_ThisOrient = Energy_ThisOrient ...
- real(EO) - abs(imag(EO));
end
end % ... and process the next scale
%% Automatically determine noise threshold
%
% Assuming the noise is Gaussian the response of the filters to noise will
% form Rayleigh distribution. We use the filter responses at the smallest
% scale as a guide to the underlying noise level because the smallest scale
% filters spend most of their time responding to noise, and only
% occasionally responding to features. Either the median, or the mode, of
% the distribution of filter responses can be used as a robust statistic to
% estimate the distribution mean and standard deviation as these are related
% to the median or mode by fixed constants. The response of the larger
% scale filters to noise can then be estimated from the smallest scale
% filter response according to their relative bandwidths.
%
% This code assumes that the expected reponse to noise on the phase congruency
% calculation is simply the sum of the expected noise responses of each of
% the filters. This is a simplistic overestimate, however these two
% quantities should be related by some constant that will depend on the
% filter bank being used. Appropriate tuning of the parameter 'k' will
% allow you to produce the desired output.
if noiseMethod >= 0 % We are using a fixed noise threshold
T = noiseMethod; % use supplied noiseMethod value as the threshold
else
% Estimate the effect of noise on the sum of the filter responses as
% the sum of estimated individual responses (this is a simplistic
% overestimate). As the estimated noise response at succesive scales
% is scaled inversely proportional to bandwidth we have a simple
% geometric sum.
totalTau = tau * (1 - (1/mult)^nscale)/(1-(1/mult));
% Calculate mean and std dev from tau using fixed relationship
% between these parameters and tau. See
% http://mathworld.wolfram.com/RayleighDistribution.html
EstNoiseEnergyMean = totalTau*sqrt(pi/2); % Expected mean and std
EstNoiseEnergySigma = totalTau*sqrt((4-pi)/2); % values of noise energy
% Noise threshold, make sure it is not less than epsilon.
T = max(EstNoiseEnergyMean + k*EstNoiseEnergySigma, epsilon);
end
% Apply noise threshold, this is effectively wavelet denoising via
% soft thresholding. Note 'Energy_ThisOrient' will have -ve values.
% These will be floored out at the final normalization stage.
Energy_ThisOrient = Energy_ThisOrient - T;
% Update accumulator matrix for sumAn and totalEnergy
totalSumAn = totalSumAn + sumAn_ThisOrient;
totalEnergy = totalEnergy + Energy_ThisOrient;
% Update orientation matrix by finding image points where the energy in
% this orientation is greater than in any previous orientation (the
% change matrix) and then replacing these elements in the orientation
% matrix with the current orientation number.
if(o == 1),
maxEnergy = Energy_ThisOrient;
else
change = Energy_ThisOrient > maxEnergy;
orientation = (o - 1).*change + orientation.*(~change);
maxEnergy = max(maxEnergy, Energy_ThisOrient);
end
end % For each orientation
% Normalize totalEnergy by the totalSumAn to obtain phase symmetry
% totalEnergy is floored at 0 to eliminate -ve values
phaseSym = max(totalEnergy, 0) ./ (totalSumAn + epsilon);
% Convert orientation matrix values to degrees
orientation = fix(orientation * (180 / norient));
%------------------------------------------------------------------
% CHECKARGS
%
% Function to process the arguments that have been supplied, assign
% default values as needed and perform basic checks.
function [im, nscale, norient, minWaveLength, mult, sigmaOnf, ...
k, polarity, noiseMethod] = checkargs(arg)
nargs = length(arg);
if nargs < 1
error('No image supplied as an argument');
end
% Set up default values for all arguments and then overwrite them
% with with any new values that may be supplied
im = [];
nscale = 5; % Number of wavelet scales.
norient = 6; % Number of filter orientations.
minWaveLength = 3; % Wavelength of smallest scale filter.
mult = 2.1; % Scaling factor between successive filters.
sigmaOnf = 0.55; % 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.
k = 2.0; % No of standard deviations of the noise
% energy beyond the mean at which we set the
% noise threshold point.
polarity = 0; % Look for both black and white spots of symmetrry
noiseMethod = -1; % Use median response of smallest scale filter
% to estimate noise characteristics.
% Allowed argument reading states
allnumeric = 1; % Numeric argument values in predefined order
keywordvalue = 2; % Arguments in the form of string keyword
% followed by numeric value
readstate = allnumeric; % Start in the allnumeric state
if readstate == allnumeric
for n = 1:nargs
if isa(arg{n},'char')
readstate = keywordvalue;
break;
else
if n == 1, im = arg{n};
elseif n == 2, nscale = arg{n};
elseif n == 3, norient = arg{n};
elseif n == 4, minWaveLength = arg{n};
elseif n == 5, mult = arg{n};
elseif n == 6, sigmaOnf = arg{n};
elseif n == 7, k = arg{n};
elseif n == 8, polarity = arg{n};
elseif n == 9,noiseMethod = arg{n};
end
end
end
end
% Code to handle parameter name - value pairs
if readstate == keywordvalue
while n < nargs
if ~isa(arg{n},'char') || ~isa(arg{n+1}, 'double')
error('There should be a parameter name - value pair');
end
if strncmpi(arg{n},'im' ,2), im = arg{n+1};
elseif strncmpi(arg{n},'nscale' ,2), nscale = arg{n+1};
elseif strncmpi(arg{n},'norient' ,4), norient = arg{n+1};
elseif strncmpi(arg{n},'minWaveLength',2), minWaveLength = arg{n+1};
elseif strncmpi(arg{n},'mult' ,2), mult = arg{n+1};
elseif strncmpi(arg{n},'sigmaOnf',2), sigmaOnf = arg{n+1};
elseif strncmpi(arg{n},'k' ,1), k = arg{n+1};
elseif strncmpi(arg{n},'polarity',2), polarity = arg{n+1};
elseif strncmpi(arg{n},'noiseMethod',4), noiseMethod = arg{n+1};
else error('Unrecognised parameter name');
end
n = n+2;
if n == nargs
error('Unmatched parameter name - value pair');
end
end
end
if isempty(im)
error('No image argument supplied');
end
if ~isa(im, 'double')
im = double(im);
end
if nscale < 1
error('nscale must be an integer >= 1');
end
if norient < 1
error('norient must be an integer >= 1');
end
if minWaveLength < 2
error('It makes little sense to have a wavelength < 2');
end
if polarity ~= -1 && polarity ~= 0 && polarity ~= 1
error('Allowed polarity values are -1, 0 and 1')
end
%%-------------------------------------------------------------------------
% RAYLEIGHMODE
%
% Computes mode of a vector/matrix of data that is assumed to come from a
% Rayleigh distribution.
%
% Usage: rmode = rayleighmode(data, nbins)
%
% Arguments: data - data assumed to come from a Rayleigh distribution
% nbins - Optional number of bins to use when forming histogram
% of the data to determine the mode.
%
% Mode is computed by forming a histogram of the data over 50 bins and then
% finding the maximum value in the histogram. Mean and standard deviation
% can then be calculated from the mode as they are related by fixed
% constants.
%
% mean = mode * sqrt(pi/2)
% std dev = mode * sqrt((4-pi)/2)
%
% See
% http://mathworld.wolfram.com/RayleighDistribution.html
% http://en.wikipedia.org/wiki/Rayleigh_distribution
%
function rmode = rayleighmode(data, nbins)
if nargin == 1
nbins = 50; % Default number of histogram bins to use
end
mx = max(data(:));
edges = 0:mx/nbins:mx;
n = histc(data(:),edges);
[dum,ind] = max(n); % Find maximum and index of maximum in histogram
rmode = (edges(ind)+edges(ind+1))/2;