function [rowS,colS,dl,dp]=trace_shape(region,S)

% shrink to include only region
% (later frame by 1 to avoid boundary problems)
Nrows=size(region,1);
[col,row]=meshgrid(1:size(region,2),1:size(region,1));
on=find(region);
minrow=min(row(on));
maxrow=max(row(on));
mincol=min(col(on));
maxcol=max(col(on));
region=region(minrow:maxrow,mincol:maxcol);

[maxrow2,maxcol2]=size(region);

region=frame(region,1);

row=floor(maxrow2/2);
col=1;
dir=0;
while region(row,col)==0, col=col+1; end

tmp=next(row,col,dir,region(row-1:row+1,col-1:col+1));
row(2)=tmp(1);
col(2)=tmp(2);
dir(2)=tmp(3);

while row(end)~=row(1) | col(end)~=col(1) | ...
       sum(next(row(end),col(end),dir(end),...
       region(row(end)-1:row(end)+1,col(end)-1:col(end)+1)) ~= ...
       [row(2),col(2),dir(2)])
   tmp=next(row(end),col(end),dir(end),...
       region(row(end)-1:row(end)+1,col(end)-1:col(end)+1));
   row(end+1)=tmp(1);
   col(end+1)=tmp(2);
   dir(end+1)=tmp(3);
end

row=row(1:end-1);
col=col(1:end-1);

% now (col,row) is a list of x,y pairs, tracing around the boundary

% take this list and turn it into a list of smoothed
% d(angle),d(length) pairs

if S>0
   smooth=gaussian(-3*S:3*S,0,S^2);
   colS=conv([col(end-3*S+1:end) col col(1:3*S)],smooth);
   colS=colS(2*3*S+1:end-2*3*S);
   rowS=conv([row(end-3*S+1:end) row row(1:3*S)],smooth);
   rowS=rowS(2*3*S+1:end-2*3*S);
else
   colS=col;
   rowS=row;
end

% now have a smoothed list of x,y pairs, tracing around the boundary
if 0 % disp on the region
   clf;imagen(region);hold on;plot(col,row,'*')
   plot(colS,rowS)
elseif 0 % show many scales
   Svect=[1 2 3 4 5 7 10 15 20 30];
   for i=1:length(Svect);subplot(4,3,i);S=Svect(i);
      smooth=gaussian(-3*S:3*S,0,S^2);
      colS=conv([col(end-3*S+1:end) col col(1:3*S)],smooth);
      colS=colS(2*3*S+1:end-2*3*S);
      rowS=conv([row(end-3*S+1:end) row row(1:3*S)],smooth);
      rowS=rowS(2*3*S+1:end-2*3*S);
      plot(colS,-rowS);title(int2str(S));axis('image')
   end
end


diff=[
   colS(2:end)-colS(1:end-1) colS(1)-colS(end)
   -1*(rowS(2:end)-rowS(1:end-1)) -1*[rowS(1)-rowS(end)]
   ];
dl=sqrt(sum(diff.^2,1));
phi=atan2(diff(2,:),diff(1,:));
dp=[0 phi(2:end)-phi(1:end-1)];

if 0 % reconstruct using dl,dp
   l=0;
   phi=0;
   xy=[0;0];
   lp=[];
   for i=1:length(dl)
      phi=phi+dp(i);
      l=l+dl(i);
      lp=[lp [l;phi]];
      xy(:,end+1)=xy(:,end)+[dl(i)*cos(phi) ; dl(i)*sin(phi)];
   end
end

colS=colS-1+(mincol-1); % 1 to undo framing, (mincol-1) to undo centering
rowS=rowS-1+(minrow-1); % 1 to undo framing, (minrow-1) to undo centering

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function out=next(row,col,dir,patch);

back=mod(dir+4,8); % direction we came from

lookup=[
   1 0
   1 1
   0 1
   -1 1
   -1 0
   -1 -1
   0 -1
   1 -1];
%% (x,y) step to go for each direction

patch=flipud(patch); % so y goes the right way

newdir=-1;
for i=[back-1:-1:0 7:-1:back]
   if patch(2+lookup(i+1,2),2+lookup(i+1,1))
      newdir=i;
      break
   end
end

if newdir<0, keyboard, end

row2=row-lookup(newdir+1,2);
col2=col+lookup(newdir+1,1);

out=[row2,col2,newdir];

return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1997-1999. The Regents of the University
% of California (Regents). All Rights Reserved.
% 
% Permission to use, copy, modify, and distribute this software and its
% documentation for educational, research, and not-for-profit purposes,
% without fee and without a signed licensing agreement, is hereby
% granted, provided that the above copyright notice, this paragraph and
% the following two paragraphs appear in all copies, modifications, and
% distributions. Contact The Office of Technology Licensing, UC Berkeley,
% 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, (510)
% 643-7201, for commercial licensing opportunities. Created by
% Chad Carson, Serge Belongie, and Jitendra Malik, Department of Electrical
% Engineering and Computer Sciences, University of California, Berkeley.
% For technical information, contact carson@eecs.berkeley.edu.
% 
% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF
% ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION
% TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

