function [xs,BC,B]=I2BC(I)
%computes boundary B of closure of binary image I
%as represented with chain code BC, 8 neighbors
%Also computes 2-column coordinate matrix xs
I=bwmorph(I,'close');
I=bwmorph(I,'clean');
I=bwmorph(I,'spur');
%I=bwlabel(I,8);
B=bwperim(I);

B=bwmorph(B,'clean');
B=bwmorph(B,'spur');
B=double(B);
%length l of perimeter
l=sum(sum(B));
%find starting pixel
v=any(B);
[Y,colindex]=max(v);
[Z,rowindex]=max(B(:,colindex));
start=[rowindex,colindex];
%introduce window W of 8-neighborhood
W=[2 3 4;1 0 5; 8 7 6];
%and vector of offsets
off=[-1 0 1];
%also 2-column matrix with x and y increments for given code
incr=[0 -1;-1 -1;-1 0;-1 1;0 1;1 1;1 0; 1 -1];

%initializes coordinate matrix
xs=zeros(l,2);

%put frame of 0's around B
[l1,l2]=size(B);
Bframe=zeros(l1+2,l2+2);
Bframe(2:l1+1,2:l2+1)=B;

%remember frame!
start=start+[1 1];
Bframe(start(1),start(2))=0;

%initialize
Neigh=Bframe(start(1)+off,start(2)+off);
chaincode=Neigh.*W;
code=max(max(chaincode));
BC=code;
next=start+incr(code,:);
t=2;

%start iterating along boundary
while t<l-1
   Neigh=Bframe(next(1)+off,next(2)+off);
      chaincode=Neigh.*W;
      code=max(max(chaincode));
      BC=[BC;code];
      
      xs(t,:)=xs(t-1,:)+incr(code,:);
      Bframe(next(1),next(2))=0;
      next=next+incr(code,:);
      t=t+1;
      
      
end      
   Bframe(start(1),start(2))=1;
   Neigh=Bframe(next(1)+off,next(2)+off);
      chaincode=Neigh.*W;
      code=max(max(chaincode));
      BC=[BC;code];