% find boundary edges of a triangulation

function B=find_boundary(tri,G)

   N=size(tri,1); % number of triangles
   
   adj=sum(G');
   count=0;
   for k=1:N
       
     if adj(k)==2
	    count=count+1;
	    row=find(G(k,:)>0);
	    a=tri(k,:);
	    b=tri(row(1),:);
	    c=tri(row(2),:);
	    d=setdiff(a,intersect(b,c));
	    B(count,:)=d;
     end % if 
      
     if adj(k)==1
	    count=count+1;
	    row=find(G(k,:)>0);
	    a=tri(k,:);
	    b=tri(row(1),:);
	    d=setdiff(a,b);
	    e=intersect(b,a);
	    B(count,:)=[d,e(1)];
	    count=count+1;
	    B(count,:)=[d,e(2)];
      end % if 

   end % for k 

 %  for k=1:size(B,1)
 %     x1=x(B(k,1));
 %     y1=y(B(k,1));
 %     x2=x(B(k,2));
 %     y2=y(B(k,2));
 %     plot([x1;x2],[y1;y2],'r');
 %  end 

return 
