function [contour_list]=contour(img)

% TURTLE ALGORITHM

% 1. Initialization
direction=1; % on the right
nb_points=0;  % point number on the contour

% 2. Research of the first black pixel - we assume that there are 1!
s=size(img);
img_contour=img;
i=1;
j=1;
while (img(i,j)~=0)
   i=i+1;
   if i>s(1)
      i=1;
      j=j+1;
   end
end

% 3. Featured departure point
p(1)=i;
p(2)=j;
start_point(1)=i;
start_point(2)=j;

% 4. current direction = direction
current_direction = direction;

% 5. Select the neighbor q of p in the current_direction indicated by the Freeman code
% (it is assumed that there is)
loop_bool=0;
while (~loop_bool)
   neighbor=0;
   while (~neighbor)
      neighbor=1;
      q(1)=p(1);
      q(2)=p(2);
      switch current_direction
      case 0
         q(2)=q(2)+1;
      case 1
         q(1)=q(1)-1;
         q(2)=q(2)+1;
      case 2
         q(1)=q(1)-1;
      case 3
         q(1)=q(1)-1;
         q(2)=q(2)-1;
      case 4
         q(2)=q(2)-1;
      case 5
         q(1)=q(1)+1;
         q(2)=q(2)-1;
      case 6
         q(1)=q(1)+1;
      case 7
         q(1)=q(1)+1;
         q(2)=q(2)+1;
      end  
     
      % 6. if q is not a point of the region
      % OR is not neighbor of a white pixel (connectivity 4 V)
      % OR has already been identified as a contour point THEN...
      if (double(img(q(1)-1,q(2)))+double(img(q(1)+1,q(2)))+double(img(q(1),q(2)-1))+double(img(q(1),q(2)+1))==0)...
              || (img(q(1),q(2))~=0)...
              || (img_contour(q(1),q(2))==128)
        
         % If this is not the starting point then change direction and back to 5
         if ~((q(1)==start_point(1)) && (q(2)==start_point(2)))
            current_direction = mod(current_direction + 7 , 8); % It rotates in the clockwise direction
            neighbor=0;
         end
      end
   end
  
   % 7. Mark P as a contour point
   nb_points=nb_points+1;
   contour_list(nb_points,1)=p(1);
   contour_list(nb_points,2)=p(2);
  
   img_contour(contour_list(nb_points,1),contour_list(nb_points,2))=128;
  
   % 8. direction = current direction
   direction = current_direction;
   current_direction=mod(current_direction+3,8); % We look for the next neighbor: most on the "left"
  
   % 9. p = q
   p(1)=q(1);
   p(2)=q(2);
  
   % 10. While p is different from start_point then back to 5
   loop_bool = (p(1)==start_point(1))*(p(2)==start_point(2));
end
