ls Gradient.m
Gradient.m

ls Hysteresis.m
Hysteresis.m

ls NMS2.m
NMS2.m

type NMS2.m

% 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,m,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);

type Hysteresis.m

% Hysteresis thresholding: keep all points with m>t connected to 
% points with m>t2 (=2.5t)
function [m1] = Hysteresis(m, t)

[L, numc] = bwlabel(m);  
t2 = 2.5*t;

[I,J] = find((L>0) & (m>t2));  % find labeled points with mag>t2
comps = zeros(numc,1);  % initialize comps

for i=1:size(I),        % components with at least one point mag.>t2
  comps(L(I(i),J(i))) = 1;
end

[I,J] = find(L>0);         % find non-zero points
for i=1:size(I),       
  L(I(i),J(i)) = comps(L(I(i),J(i)));  %replace connected components labels
end;
   
m1 = zeros(size(m));     % initialize m
I2 = find(L>0);          % L can be treated as a vector
m1(I2) = m(I2);

% for i=1:numc,    % remove connected components that are too small
%     i, comps(i)
%     if comps(i)>0
%         ind = find(L == i);
%         i, length(ind)
%         if (length(ind) < MinEdgeLen) 
%             comps(i) = 0;
%         end;
%     end;
% end;
        



ls *.jpg
image1.jpg	image2.jpg	s_289.jpg	s_290.jpg

I = imread('s_290.jpg')



Ig = rgb2gray(I);
type Gradient.m

% Partial derivatives in X and Y directions
function [Gx,Gy] = Gradient(A, p);

X = [-1 0 1]/2; % filter
Y = [1 0 -1]'/2; % filter

Gx = filter2(X,A); % correlation
Gy = filter2(Y,A);
[M,N] = size(A);
Gx(1:p+1,:) = 0;
Gx(M-p-1:M,:) = 0;
Gx(:,1:p+1) = 0;
Gx(:,N-p-1:N) = 0;
Gy(1:p+1,:) = 0;
Gy(M-p-1:M,:) = 0;
Gy(:,1:p+1) = 0;
Gy(:,N-p-1:N) = 0;


[Gx,Gy] = Gradient(Ig, 4)
m_nms = NMS2(Gx,Gy,m,3);
m_h = Hysteresis(m_nms,3);
figure, imshow(m_h > 0);
[I,J] = find(m_h>0);
In = find(m_h > 0);
[M,N] = size(m_h)

M =

   480


N =

   640

quiver(J,M-I+1,Gx(In),Gy(In),1)
diary off;