% 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 compatibilityfunction [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