1 | % FASTRADIAL - Loy and Zelinski's fast radial feature detector |
---|
2 | % |
---|
3 | % An implementation of Loy and Zelinski's fast radial feature detector |
---|
4 | % |
---|
5 | % Usage: S = fastradial(im, radii, alpha) |
---|
6 | % |
---|
7 | % Arguments: |
---|
8 | % im - image to be analysed |
---|
9 | % radii - array of integer radius values to be processed |
---|
10 | % suggested radii might be [1 3 5] |
---|
11 | % alpha - radial strictness parameter. |
---|
12 | % 1 - slack, accepts features with bilateral symmetry. |
---|
13 | % 2 - a reasonable compromise. |
---|
14 | % 3 - strict, only accepts radial symmetry. |
---|
15 | % ... and you can go higher |
---|
16 | % |
---|
17 | % Returns S - Symmetry map. Bright points with high symmetry are |
---|
18 | % marked with large positive values. Dark points of |
---|
19 | % high symmetry marked with large -ve values. |
---|
20 | % |
---|
21 | % To localize points use NONMAXSUPPTS on S, -S or abs(S) depending on |
---|
22 | % what you are seeking to find. |
---|
23 | |
---|
24 | % Reference: |
---|
25 | % Loy, G. Zelinsky, A. Fast radial symmetry for detecting points of |
---|
26 | % interest. IEEE PAMI, Vol. 25, No. 8, August 2003. pp 959-973. |
---|
27 | |
---|
28 | % Copyright (c) 2004-2005 Peter Kovesi |
---|
29 | % School of Computer Science & Software Engineering |
---|
30 | % The University of Western Australia |
---|
31 | % http://www.csse.uwa.edu.au/ |
---|
32 | % |
---|
33 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
34 | % of this software and associated documentation files (the "Software"), to deal |
---|
35 | % in the Software without restriction, subject to the following conditions: |
---|
36 | % |
---|
37 | % The above copyright notice and this permission notice shall be included in |
---|
38 | % all copies or substantial portions of the Software. |
---|
39 | % |
---|
40 | % The Software is provided "as is", without warranty of any kind. |
---|
41 | |
---|
42 | % November 2004 - original version |
---|
43 | % July 2005 - Bug corrected: magitude and orientation matrices were |
---|
44 | % not zeroed for each radius value used (Thanks to Ben |
---|
45 | % Jackson) |
---|
46 | |
---|
47 | function S = fastradial(im, radii, alpha) |
---|
48 | |
---|
49 | if any(radii ~= round(radii)) | any(radii < 1) |
---|
50 | error('radii must be integers and > 1') |
---|
51 | end |
---|
52 | |
---|
53 | [rows,cols]=size(im); |
---|
54 | |
---|
55 | % Use the Sobel masks to get gradients in x and y |
---|
56 | gx = [-1 0 1 |
---|
57 | -2 0 2 |
---|
58 | -1 0 1]; |
---|
59 | gy = gx'; |
---|
60 | |
---|
61 | imgx = filter2(gx,im); |
---|
62 | imgy = filter2(gy,im); |
---|
63 | mag = sqrt(imgx.^2 + imgy.^2)+eps; % (+eps to avoid division by 0) |
---|
64 | |
---|
65 | % Normalise gradient values so that [imgx imgy] form unit |
---|
66 | % direction vectors. |
---|
67 | imgx = imgx./mag; |
---|
68 | imgy = imgy./mag; |
---|
69 | |
---|
70 | S = zeros(rows,cols); % Symmetry matrix |
---|
71 | |
---|
72 | [x,y] = meshgrid(1:cols, 1:rows); |
---|
73 | |
---|
74 | for n = radii |
---|
75 | M = zeros(rows,cols); % Magnitude projection image |
---|
76 | O = zeros(rows,cols); % Orientation projection image |
---|
77 | |
---|
78 | % Coordinates of 'positively' and 'negatively' affected pixels |
---|
79 | posx = x + round(n*imgx); |
---|
80 | posy = y + round(n*imgy); |
---|
81 | |
---|
82 | negx = x - round(n*imgx); |
---|
83 | negy = y - round(n*imgy); |
---|
84 | |
---|
85 | % Clamp coordinate values to range [1 rows 1 cols] |
---|
86 | posx( find(posx<1) ) = 1; |
---|
87 | posx( find(posx>cols) ) = cols; |
---|
88 | posy( find(posy<1) ) = 1; |
---|
89 | posy( find(posy>rows) ) = rows; |
---|
90 | |
---|
91 | negx( find(negx<1) ) = 1; |
---|
92 | negx( find(negx>cols) ) = cols; |
---|
93 | negy( find(negy<1) ) = 1; |
---|
94 | negy( find(negy>rows) ) = rows; |
---|
95 | |
---|
96 | |
---|
97 | % Form the orientation and magnitude projection matrices |
---|
98 | for r = 1:rows |
---|
99 | for c = 1:cols |
---|
100 | O(posy(r,c),posx(r,c)) = O(posy(r,c),posx(r,c)) + 1; |
---|
101 | O(negy(r,c),negx(r,c)) = O(negy(r,c),negx(r,c)) - 1; |
---|
102 | |
---|
103 | M(posy(r,c),posx(r,c)) = M(posy(r,c),posx(r,c)) + mag(r,c); |
---|
104 | M(negy(r,c),negx(r,c)) = M(negy(r,c),negx(r,c)) - mag(r,c); |
---|
105 | end |
---|
106 | end |
---|
107 | |
---|
108 | % Clamp Orientation projection matrix values to a maximum of |
---|
109 | % +/-kappa, but first set the normalization parameter kappa to the |
---|
110 | % values suggested by Loy and Zelinski |
---|
111 | if n == 1, kappa = 8; else, kappa = 9.9; end |
---|
112 | |
---|
113 | O(find(O > kappa)) = kappa; |
---|
114 | O(find(O < -kappa)) = -kappa; |
---|
115 | |
---|
116 | % Unsmoothed symmetry measure at this radius value |
---|
117 | F = M./kappa .* (abs(O)/kappa).^alpha; |
---|
118 | |
---|
119 | % Generate a Gaussian of size proportional to n to smooth and spread |
---|
120 | % the symmetry measure. The Gaussian is also scaled in magnitude |
---|
121 | % by n so that large scales do not lose their relative weighting. |
---|
122 | A = fspecial('gaussian',[n n], 0.25*n) * n; |
---|
123 | |
---|
124 | S = S + filter2(A,F); |
---|
125 | |
---|
126 | end % for each radius |
---|
127 | |
---|
128 | S = S/length(radii); % Average |
---|