% This is vectorized version

% 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] = NMS2(Gx,Gy,t);

[M,N] = size(Gx);

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

[I,J] = find(m>0);
aGx = abs(Gx);
aGy = abs(Gy);

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

for i=1:size(I),
	if I(i)>1 & I(i)<M & J(i)>1 & J(i)<N
           if Gx(I(i),J(i))>0 
              dx = 1;
           else 
              dx = -1;
           end;
           if Gy(I(i),J(i))>0 
              dy = 1;
           else
	      dy = -1;
           end;

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

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