% NONMAXSUP - Non-maxima suppression%% Usage:%          [im,location] = nonmaxsup(inimage, orient, radius);%% Function for performing non-maxima suppression on an image using an% orientation image.  It is assumed that the orientation image gives % feature normal orientation angles in degrees (0-180).%% Input:%   inimage - Image to be non-maxima suppressed.% %   orient  - Image containing feature normal orientation angles in degrees%             (0-180), angles positive anti-clockwise.% %   radius  - Distance in pixel units to be looked at on each side of each%             pixel when determining whether it is a local maxima or not.%             This value cannot be less than 1.%             (Suggested value about 1.2 - 1.5)%% Returns:%   im        - Non maximally suppressed image.%   location  - Complex valued image holding subpixel locations of edge%               points. For any pixel the real part holds the subpixel row%               coordinate of that edge point and the imaginary part holds%               the column coordinate.  (If a pixel value is 0+0i then it%               is not an edgepoint.)%               (Note that if this function is called without 'location'%               being specified as an output argument is not computed)%% Notes:%% The suggested radius value is 1.2 - 1.5 for the following reason. If the% radius parameter is set to 1 there is a chance that a maxima will not be% identified on a broad peak where adjacent pixels have the same value.  To% overcome this one typically uses a radius value of 1.2 to 1.5.  However% under these conditions there will be cases where two adjacent pixels will% both be marked as maxima.  Accordingly there is a final morphological% thinning step to correct this.%% This function is slow.  It uses bilinear interpolation to estimate% intensity values at ideal, real-valued pixel locations on each side of% pixels to determine if they are local maxima.% Copyright (c) 1996-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.% December  1996 - Original version% September 2004 - Subpixel localization added% August    2005 - Made Octave compatiblefunction [im, location] = nonmaxsup(inimage, orient, radius)if any(size(inimage) ~= size(orient))  error('image and orientation image are of different sizes');endif radius < 1  error('radius must be >= 1');endOctave = exist('OCTAVE_VERSION') ~= 0;  % Are we running under Octave?[rows,cols] = size(inimage);im = zeros(rows,cols);        % Preallocate memory for output imageif nargout == 2    location = zeros(rows,cols);endiradius = ceil(radius);% Precalculate x and y offsets relative to centre pixel for each orientation angle angle = [0:180].*pi/180;    % Array of angles in 1 degree increments (but in radians).xoff = radius*cos(angle);   % x and y offset of points at specified radius and angleyoff = radius*sin(angle);   % from each reference position.hfrac = xoff - floor(xoff); % Fractional offset of xoff relative to integer locationvfrac = yoff - floor(yoff); % Fractional offset of yoff relative to integer locationorient = fix(orient)+1;     % Orientations start at 0 degrees but arrays start                            % with index 1.% Now run through the image interpolating grey values on each side% of the centre pixel to be used for the non-maximal suppression.for row = (iradius+1):(rows - iradius)  for col = (iradius+1):(cols - iradius)     or = orient(row,col);   % Index into precomputed arrays    x = col + xoff(or);     % x, y location on one side of the point in question    y = row - yoff(or);    fx = floor(x);          % Get integer pixel locations that surround location x,y    cx = ceil(x);    fy = floor(y);    cy = ceil(y);    tl = inimage(fy,fx);    % Value at top left integer pixel location.    tr = inimage(fy,cx);    % top right    bl = inimage(cy,fx);    % bottom left    br = inimage(cy,cx);    % bottom right    upperavg = tl + hfrac(or) * (tr - tl);  % Now use bilinear interpolation to    loweravg = bl + hfrac(or) * (br - bl);  % estimate value at x,y    v1 = upperavg + vfrac(or) * (loweravg - upperavg);  if inimage(row, col) > v1 % We need to check the value on the other side...    x = col - xoff(or);     % x, y location on the `other side' of the point in question    y = row + yoff(or);    fx = floor(x);    cx = ceil(x);    fy = floor(y);    cy = ceil(y);    tl = inimage(fy,fx);    % Value at top left integer pixel location.    tr = inimage(fy,cx);    % top right    bl = inimage(cy,fx);    % bottom left    br = inimage(cy,cx);    % bottom right    upperavg = tl + hfrac(or) * (tr - tl);    loweravg = bl + hfrac(or) * (br - bl);    v2 = upperavg + vfrac(or) * (loweravg - upperavg);    if inimage(row,col) > v2            % This is a local maximum.      im(row, col) = inimage(row, col); % Record value in the output                                        % image.					      % Code for sub-pixel localization if it was requested					      if nargout == 2         % Solve for coefficients of parabola that passes through 	 % [-1, v1]  [0, inimage] and [1, v2]. 	 % v = a*r^2 + b*r + c	 c = inimage(row,col);	 a = (v1 + v2)/2 - c;	 b = a + c - v1;	 	 % location where maxima of fitted parabola occurs	 r = -b/(2*a);	 location(row,col) = complex(row + r*yoff(or), col - r*xoff(or));            end          end   end  endend% Finally thin the 'nonmaximally suppressed' image by pointwise% multiplying itself with a morphological skeletonization of itself.%% I know it is oxymoronic to thin a nonmaximally supressed image but % fixes the multiple adjacent peaks that can arise from using a radius% value > 1.if Octave    skel = bwmorph(im>0,'thin',Inf);   % Octave bwmorph only works on binary                                       % images and 'thin' seems to produce                                       % better results.else    skel = bwmorph(im,'skel',Inf);endim = im.*skel;if nargout == 2    location = location.*skel;end