function E=TraceEdge(A,pos,SHOW)

% E=TraceEdge(A,y,x,jj)
%
% INPUT:
% A    [NxM]   a binary image.
% pos  [1x2]   the coordinates (x,y) of a point inside the current patch (a region 
%              of ones to be traced). 
%              pos(1) is the x (across) index and pos(2) the y (down) index.
% SHOW [1x1]   (Optional) If SHOW=1, the boundary will be plotted on the
%               current image. If SHOW=0, or omitted, nothing will be plotted.
%
% OUTPUT:
% E     [Kx2]   the list of edge coordinates of the current patch, obtained
%               counter clockwise from the first boundary point to the right of pos. 


if nargin==2, SHOW=0; end


si=size(A);
Ny=si(1); Nx=si(2);
if (pos(1)<0) | (pos(1)>Ny) | (pos(2)<0) | (pos(2)>Nx), E=[], return; end
if A(pos(2),pos(1))==0, E=[]; return; end

p=FirstPosition(A,[pos(2) pos(1)]');
p0=p;
 
E=p0([2 1])'; 
jj=2;
      
[p,jj]=NextPosition(A,p,jj);
E=[E; p([2 1])'];
if SHOW>0, 
    hold on
    plot(E(1,1),E(1,2),'bs','Markersize',5,'MarkerfaceColor','b');
    plot(E(end,1),E(end,2),'om','Markersize',5,'MarkerfaceColor','m'); 
end

while ~((p(1)==p0(1)) & (p(2)==p0(2)))
    [p,jj]=NextPosition(A,p,jj);
    E=[E; p([2 1])'];
    if SHOW>0, plot(E(end,1),E(end,2),'m.'); end
end 


E=E(1:end-1,:);
if SHOW>0, plot(E(:,1),E(:,2),'m'); end

% ================================================
function [p,jj]=NextPosition(S,p,jj)
si=size(S); Nx=si(2); Ny=si(1);
NN=[0 1 0 -1; 1 0 -1 0];
jj=mod(jj-1,4);
n=NN(:,jj+1);
q=p+n;  
for j=1:4
   if (q(1)>Ny) | (q(1)<1) | (q(2)>Nx) | (q(2)<1) | S(q(1),q(2))==0
       jj=mod(jj+1,4);
       n=NN(:,jj+1);
       q=p+n; 
    end    
end    
p=p+n;

% ====================================================

function p=FirstPosition(S,p)
Nx=size(S,2);
for j=p(1):size(S,2)
    if (p(2)+1>Nx), return;  end
    if S(p(1),p(2)+1)==0, return;  end
    p(2)=p(2)+1;
end