% Non-maximum suppression for a gradient image;
% magnitudes under t are ignored
% returns a binary image where maximal gradient values are non-zero
function [m1] = NMS(Gx,Gy,t);

m = sqrt(Gx.^2+Gy.^2);  % magnitude
I1 = find(m<t);
m(I1) = 0;    % zero values under t

aGx = abs(Gx);
aGy = abs(Gy);

[M,N] = size(m);
Z = zeros([M,N]);  % binary image

for i=2:M-1, 
    for j=2:N-1,
      % find direction
      if m(i,j)>0   % for all points for which m>t do
           if Gx(i,j)>0 
              dx = 1;
           else 
              dx = -1;
           end;
           if Gy(i,j)>0 
              dy = 1;
           else
	      dy = -1;
           end;

        if aGx(i,j)>aGy(i,j) % x is dominant direction
	   m01 = m(i,j+dx);
           m11 = m(i-dy,j+dx);
           m_01 = m(i,j-dx);
           m_11 = m(i+dy,j-dx);
           m1 = (aGx(i,j)-aGy(i,j))*m01 + aGy(i,j)*m11;
           m2 = (aGx(i,j)-aGy(i,j))*m_01 + aGy(i,j)*m_11;
           mag = m(i,j)*aGx(i,j);
           if (mag>m2) & (mag>=m1) 
              Z(i,j) = 1;
           end
	else % y is dominant direction
	   m01 = m(i-dy,j);
           m11 = m(i-dy,j+dx);
           m_01 = m(i+dy,j);
           m_11 = m(i+dy,j-dx);
           m1 = (aGy(i,j)-aGx(i,j))*m01 + aGx(i,j)*m11;
           m2 = (aGy(i,j)-aGx(i,j))*m_01 + aGx(i,j)*m_11;
           mag = m(i,j)*aGy(i,j);
           if (mag>m2) & (mag>=m1) 
              Z(i,j) = 1;
	   end
	end
      
	   end % if (m,j)>0
    end % for
end % for

I1 = find(Z==0);  % only keep points where Z=1
m(I1)=0;
m1 = uint8(m);