[37] | 1 | % NONMAXSUPPTS - Non-maximal suppression for features/corners |
---|
| 2 | % |
---|
| 3 | % Non maxima suppression and thresholding for points generated by a feature |
---|
| 4 | % or corner detector. |
---|
| 5 | % |
---|
| 6 | % Usage: [r,c] = nonmaxsuppts(cim, radius, thresh, im) |
---|
| 7 | % / |
---|
| 8 | % optional |
---|
| 9 | % |
---|
| 10 | % [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im) |
---|
| 11 | % |
---|
| 12 | % Arguments: |
---|
| 13 | % cim - corner strength image. |
---|
| 14 | % radius - radius of region considered in non-maximal |
---|
| 15 | % suppression. Typical values to use might |
---|
| 16 | % be 1-3 pixels. |
---|
| 17 | % thresh - threshold. |
---|
| 18 | % im - optional image data. If this is supplied the |
---|
| 19 | % thresholded corners are overlayed on this |
---|
| 20 | % image. This can be useful for parameter tuning. |
---|
| 21 | % Returns: |
---|
| 22 | % r - row coordinates of corner points (integer valued). |
---|
| 23 | % c - column coordinates of corner points. |
---|
| 24 | % rsubp - If four return values are requested sub-pixel |
---|
| 25 | % csubp - localization of feature points is attempted and |
---|
| 26 | % returned as an additional set of floating point |
---|
| 27 | % coords. Note that you may still want to use the integer |
---|
| 28 | % valued coords to specify centres of correlation windows |
---|
| 29 | % for feature matching. |
---|
| 30 | % |
---|
| 31 | |
---|
| 32 | % Copyright (c) 2003-2005 Peter Kovesi |
---|
| 33 | % School of Computer Science & Software Engineering |
---|
| 34 | % The University of Western Australia |
---|
| 35 | % http://www.csse.uwa.edu.au/ |
---|
| 36 | % |
---|
| 37 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
| 38 | % of this software and associated documentation files (the "Software"), to deal |
---|
| 39 | % in the Software without restriction, subject to the following conditions: |
---|
| 40 | % |
---|
| 41 | % The above copyright notice and this permission notice shall be included in all |
---|
| 42 | % copies or substantial portions of the Software. |
---|
| 43 | % |
---|
| 44 | % The Software is provided "as is", without warranty of any kind. |
---|
| 45 | |
---|
| 46 | % September 2003 Original version |
---|
| 47 | % August 2005 Subpixel localization and Octave compatibility |
---|
| 48 | |
---|
| 49 | function [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im) |
---|
| 50 | |
---|
| 51 | v = version; Octave = v(1)<'5'; % Crude Octave test |
---|
| 52 | subPixel = nargout == 4; % We want sub-pixel locations |
---|
| 53 | [rows,cols] = size(cim); |
---|
| 54 | |
---|
| 55 | % Extract local maxima by performing a grey scale morphological |
---|
| 56 | % dilation and then finding points in the corner strength image that |
---|
| 57 | % match the dilated image and are also greater than the threshold. |
---|
| 58 | |
---|
| 59 | sze = 2*radius+1; % Size of dilation mask. |
---|
| 60 | mx = ordfilt2(cim,sze^2,ones(sze)); % Grey-scale dilate. |
---|
| 61 | |
---|
| 62 | % Make mask to exclude points within radius of the image boundary. |
---|
| 63 | bordermask = zeros(size(cim)); |
---|
| 64 | bordermask(radius+1:end-radius, radius+1:end-radius) = 1; |
---|
| 65 | |
---|
| 66 | % Find maxima, threshold, and apply bordermask |
---|
| 67 | cimmx = (cim==mx) & (cim>thresh) & bordermask; |
---|
| 68 | |
---|
| 69 | [r,c] = find(cimmx); % Find row,col coords. |
---|
| 70 | |
---|
| 71 | |
---|
| 72 | if subPixel % Compute local maxima to sub pixel accuracy |
---|
| 73 | if ~isempty(r) % ...if we have some ponts to work with |
---|
| 74 | |
---|
| 75 | ind = sub2ind(size(cim),r,c); % 1D indices of feature points |
---|
| 76 | w = 1; % Width that we look out on each side of the feature |
---|
| 77 | % point to fit a local parabola |
---|
| 78 | |
---|
| 79 | % Indices of points above, below, left and right of feature point |
---|
| 80 | indrminus1 = max(ind-w,1); |
---|
| 81 | indrplus1 = min(ind+w,rows*cols); |
---|
| 82 | indcminus1 = max(ind-w*rows,1); |
---|
| 83 | indcplus1 = min(ind+w*rows,rows*cols); |
---|
| 84 | |
---|
| 85 | % Solve for quadratic down rows |
---|
| 86 | cy = cim(ind); |
---|
| 87 | ay = (cim(indrminus1) + cim(indrplus1))/2 - cy; |
---|
| 88 | by = ay + cy - cim(indrminus1); |
---|
| 89 | rowshift = -w*by./(2*ay); % Maxima of quadradic |
---|
| 90 | |
---|
| 91 | % Solve for quadratic across columns |
---|
| 92 | cx = cim(ind); |
---|
| 93 | ax = (cim(indcminus1) + cim(indcplus1))/2 - cx; |
---|
| 94 | bx = ax + cx - cim(indcminus1); |
---|
| 95 | colshift = -w*bx./(2*ax); % Maxima of quadradic |
---|
| 96 | |
---|
| 97 | rsubp = r+rowshift; % Add subpixel corrections to original row |
---|
| 98 | csubp = c+colshift; % and column coords. |
---|
| 99 | else |
---|
| 100 | rsubp = []; csubp = []; |
---|
| 101 | end |
---|
| 102 | end |
---|
| 103 | |
---|
| 104 | if nargin==4 & ~isempty(r) % Overlay corners on supplied image. |
---|
| 105 | if Octave |
---|
| 106 | warning('Only able to display points under Octave'); |
---|
| 107 | if subPixel |
---|
| 108 | plot(csubp,rsubp,'r+'), title('corners detected'); |
---|
| 109 | else |
---|
| 110 | plot(c,r,'r+'), title('corners detected'); |
---|
| 111 | end |
---|
| 112 | [rows,cols] = size(cim); |
---|
| 113 | axis([1 cols 1 rows]); axis('equal'); axis('ij') |
---|
| 114 | else |
---|
| 115 | figure(1), imshow(im,[]), hold on |
---|
| 116 | if subPixel |
---|
| 117 | plot(csubp,rsubp,'r+'), title('corners detected'); |
---|
| 118 | else |
---|
| 119 | plot(c,r,'r+'), title('corners detected'); |
---|
| 120 | end |
---|
| 121 | end |
---|
| 122 | end |
---|
| 123 | |
---|