[37] | 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 |
---|