source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/kovesi/nonmaxsuppts.m @ 37

% 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