function [xc,yc]=contour(Ib)
%[xc,yc]=contour(Ib)
%
%Searches the boundary for an object in the binary image Ib.
%(where 1 is considered as object pixel), starting at the 
%first object pixel by scanning downwards each column from
%the left.
%
%xc and yc are the output boundary coordinates.
%

% Mats Kvarström 010110 modified 020123
% Modified 030121 by John Gustafsson

%%% Notation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (xc,yc): coordinates of boundary in image
% xx     : single index notation (i.e. xx = (xc-1)*Ny+yc)
% d      : denotes the position of a boundary pixel in relation to
%          the last boundary pixel in a counterclockwise manner,
%          such that
%                 d = | 0     | 1     | 2    | 3       
%          --------------------------------------------
%           came from | right | above | left | below
% i      : a counting index

[Ny,Nx] = size(Ib);

%%% Main program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First Point
tmp    = find(Ib==1); % look up the positions of the 1's
xx(1)  = tmp(1);      % take the first element in this vector as start
d(1)   = 1; % ...as if we came from above   

% Second point
[xxnew,dnew] = next_cand(xx(1),d(1),Ib,Ny);
xx(2) = xxnew;
d(2)  = dnew;

% Third point
[xxnew,dnew] = next_cand(xx(2),d(2),Ib,Ny);

% And iterate...
i = 2;
while xx(i)~=xx(1) | xxnew~=xx(2)
  i  = i+1;
  xx(i) = xxnew;
  d(i)  = dnew;
  % next candidate
  [xxnew,dnew]=next_cand(xx(i),d(i),Ib,Ny);
end

% Convert from single index to ordinary coordinates
yc = mod(xx-1,Ny)+1;
xc = (xx-yc)/Ny+1;


%%% Subfunctions %%%%%%%%%%%%%%%%%%%d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
function [xnew,dnew] = next_cand(x,d,I,Ny)
% find the next candidate for boundary pixel
dnew = mod(d+1,4);
xnew = move(x,dnew,Ny);
while I(xnew)~=1 
  % If xnew is not a part of the object then try 
  % the next neighbour:
  dnew = mod(dnew-1,4);
  xnew = move(x,dnew,Ny);
end

%%%
function xnew = move(x,d,Ny)
% Moves from x to xnew according to d
%
if d==0     % came from neigbour above: move to the right
  xnew = x+Ny;
elseif d==1 % came from neighbour to the right: move below
  xnew = x+1;
elseif d==2 % came from neighbour below: move to the left
  xnew = x-Ny;
else        % came form neighbour to the left: move up
  xnew = x-1;
end