% NONMAXSUPPTS - Non-maximal suppression for features/corners
%
% Non maxima suppression and thresholding for points generated by a feature
% or corner detector.
%
% Usage: [r,c] = nonmaxsuppts(cim, radius, thresh, im)
% /
% optional
%
% [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)
%
% Arguments:
% cim - corner strength image.
% radius - radius of region considered in non-maximal
% suppression. Typical values to use might
% be 1-3 pixels.
% thresh - threshold.
% im - optional image data. If this is supplied the
% thresholded corners are overlayed on this
% image. This can be useful for parameter tuning.
% Returns:
% r - row coordinates of corner points (integer valued).
% c - column coordinates of corner points.
% rsubp - If four return values are requested sub-pixel
% csubp - localization of feature points is attempted and
% returned as an additional set of floating point
% coords. Note that you may still want to use the integer
% valued coords to specify centres of correlation windows
% for feature matching.
%
% Copyright (c) 2003-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% 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.
% September 2003 Original version
% August 2005 Subpixel localization and Octave compatibility
function [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)
v = version; Octave = v(1)<'5'; % Crude Octave test
subPixel = nargout == 4; % We want sub-pixel locations
[rows,cols] = size(cim);
% Extract local maxima by performing a grey scale morphological
% dilation and then finding points in the corner strength image that
% match the dilated image and are also greater than the threshold.
sze = 2*radius+1; % Size of dilation mask.
mx = ordfilt2(cim,sze^2,ones(sze)); % Grey-scale dilate.
% Make mask to exclude points within radius of the image boundary.
bordermask = zeros(size(cim));
bordermask(radius+1:end-radius, radius+1:end-radius) = 1;
% Find maxima, threshold, and apply bordermask
cimmx = (cim==mx) & (cim>thresh) & bordermask;
[r,c] = find(cimmx); % Find row,col coords.
if subPixel % Compute local maxima to sub pixel accuracy
if ~isempty(r) % ...if we have some ponts to work with
ind = sub2ind(size(cim),r,c); % 1D indices of feature points
w = 1; % Width that we look out on each side of the feature
% point to fit a local parabola
% Indices of points above, below, left and right of feature point
indrminus1 = max(ind-w,1);
indrplus1 = min(ind+w,rows*cols);
indcminus1 = max(ind-w*rows,1);
indcplus1 = min(ind+w*rows,rows*cols);
% Solve for quadratic down rows
cy = cim(ind);
ay = (cim(indrminus1) + cim(indrplus1))/2 - cy;
by = ay + cy - cim(indrminus1);
rowshift = -w*by./(2*ay); % Maxima of quadradic
% Solve for quadratic across columns
cx = cim(ind);
ax = (cim(indcminus1) + cim(indcplus1))/2 - cx;
bx = ax + cx - cim(indcminus1);
colshift = -w*bx./(2*ax); % Maxima of quadradic
rsubp = r+rowshift; % Add subpixel corrections to original row
csubp = c+colshift; % and column coords.
else
rsubp = []; csubp = [];
end
end
if nargin==4 & ~isempty(r) % Overlay corners on supplied image.
if Octave
warning('Only able to display points under Octave');
if subPixel
plot(csubp,rsubp,'r+'), title('corners detected');
else
plot(c,r,'r+'), title('corners detected');
end
[rows,cols] = size(cim);
axis([1 cols 1 rows]); axis('equal'); axis('ij')
else
figure(1), imshow(im,[]), hold on
if subPixel
plot(csubp,rsubp,'r+'), title('corners detected');
else
plot(c,r,'r+'), title('corners detected');
end
end
end
