//**********************************
//***      skeleton
// Written by Jarek Rossignac,  June 2006 
//**********************************

int[] skeleton = new int [maxnt];          // 0 = water, 1=bridge, 2=skeleton
int[] cluster = new int [maxnt];
boolean[] visited = new boolean [maxnt];
boolean[] visitedV = new boolean [maxnv];
int nclusters=1;
int rec=0;
int tic=0, maxTic=0;
int bestCluster=1;


void writeLake(int L) {
  boolean found=false;
  for (int i=0; i<nc; i++) {
   if ( !found  && (skeleton[t(i)]==0) && (cluster[t(i)]==L) && nb(i)  && (skeleton[t(o(i))]!=0) ) {c=i; found=true;};};
   int f=n(c);
   int cc=1;
   print("lake "+L+" = (");
   int fin=0;
   while (fin<1) {    
      cc=1;  while(nb(f)&&skeleton[t(o(f))]==0) {f=n(o(f)); cc++;};  print(","+cc);
      if (f==c) {fin++;};
      f=n(f);
      }; 
  println(")");
  }

boolean isSkelet (int i) {return(skeleton[t(i)]==2);}

void flattenLake() {
  maxTic++;
  for (int t=0; t<nt; t++) {visited[t]=false;};
  for (int v=0; v<nv; v++) {visitedV[v]=false;};
  c=0; while (!((skeleton[t(c)]==0)&&(cluster[t(c)]==lake))) {c++;};
  visited[t(c)]=true; visitedV[v(c)]=true; visitedV[v(n(c))]=true;  visitedV[v(p(c))]=true; 
  vec CB = g(p(c)).vecTo(g(n(c)));
  vec N=cross( CB , g(p(c)).vecTo(g(c)) );
  vec T=cross ( N , CB ); T.makeUnit();
  pt M=midPt(g(n(c)),g(p(c))); 
  M.addScaledVec( CB.norm()*sqrt(3)/2.0,T);
  g(c).setToPoint(M);
//  g(c).moveTowards(0.3,M);
  tic=0;
  flattenTri(c,1);
  tic=0;
  flattenTri(p(c),1); 
  tic=0;
  flattenTri(n(c),1);
  println("done flattening lake "+lake);
  }

void flattenTri(int c, int r) {  
  int o=o(c);
  if (nb(c)&&(!visited[t(o)])&&(skeleton[t(o)]==0)&&(cluster[t(o)]==lake)) {
     if (!visitedV[v(o)]) {
       pt T=g(p(c)).make(); T.addVec(g(c).vecTo(g(n(c)))); g(o).moveTowards(1,T);  
       visitedV[v(o)]=true;};
     visited[t(o)]=true;
     flattenTri(n(o),r+1);  flattenTri(p(o),r+1); 
    };
  };

void computeSkeleton() {
  int ts=0, tc=0;
  int r=1;
  int singV=0; for (int i=0; i<nv; i++) {if (Valence[i]!=6) {singV++;};};
  for(int i=0; i<nt; i++) {Mt[i]=0; skeleton[i]=0;};  
  for(int i=0; i<nc; i++) {if (Valence[v(i)]!=6) {Mt[t(i)]=1; skeleton[t(i)]=2;};};
  int singT=0; for (int i=0; i<nt; i++) {if (skeleton[i]==2) {singT++;};};

  for(int i=0; i<nc; i++) {
     if (!isSkelet(i)) {
         int nbrs=0; 
         if((o(i)!=-1)&&(isSkelet(o(i)))) {nbrs++;};
         if((r(i)!=-1)&&(isSkelet(r(i)))) {nbrs++;};
         if((l(i)!=-1)&&(isSkelet(l(i)))) {nbrs++;};
         if (nbrs>2) {Mt[t(i)]=1; skeleton[t(i)]=2; };
         };
     };
  ts=0; tc=0;   for(int i=0; i<nt; i++) {if(skeleton[i]==2) {ts++;} else {tc++;}; }; 
  println(nv+" vertices ("+singV+" irregular). "+nt+" triangles ("+singT+" incident on irregular vertices, "+ts+" after adding sinlge-triangle lakes)");
  for(int i=0; i<nt; i++) {cluster[i]=0;};  
  makeClusters();    
  singT=0; for (int i=0; i<nt; i++) {if (skeleton[i]==2) {singT++;};};
  println(nclusters+" islands with a total of "+singT+" triangles");
  for (int i=1; i<nclusters; i++) {int other=growCluster(bestCluster);  traceBridge(bestCluster); mergeCluster(other,bestCluster); };
  singT=0; for (int i=0; i<nt; i++) {if (skeleton[i]==2) {singT++;};};
  println("Skeleton (islands+bridges) has "+singT+" triangles");
  manifoldSkeleton();
  singT=0; for (int i=0; i<nt; i++) {if (skeleton[i]==2) {singT++;};};
  println("Manifold skeleton has a total of "+singT+" triangles");
  paintLakes();  
  }
   
void manifoldSkeleton() {
  for(int c=0; c<nc; c++) {
    if ((skeleton[t(c)]==0)&&(skeleton[t(l(c))]==2)&&(skeleton[t(r(c))]==2)) { 
      int b=p(r(c));
      boolean gap=false;
      while (b!=c) {if (skeleton[t(b)]==0) {gap=true;}; b=p(r(b));};
      if (gap) {skeleton[t(c)]=2; };
      };  
    };
  }
   
void paintLakes() {
  nclusters=0;
  tic=0;
  for(int t=0; t<nt; t++) {if (skeleton[t]==0) {cluster[t]=0; Mt[t]=0;} else {tic++;}; };  
  println("The mesh has "+tic+" skeleton triangles and "+(nt-tic)+" water triangles separated into lakes:");
 
  for (int t=0; t<nt; t++) {
     int c=t*3; 
     if ( (skeleton[t]==0) && (cluster[t]==0) ) {
       tic=-2;
       nclusters++;   
       traceLake(t*3,nclusters); cluster[t(c)]=0;
       traceLake(t*3+1,nclusters); cluster[t(c)]=0;
       traceLake(t*3+2,nclusters); 
       println("  Lake "+nclusters+" has "+tic+" triangles");
       };
     };
   rings=nclusters+1;
   }  
   
void traceLake(int c, int n) {
  if ((skeleton[t(c)]==0) && (cluster[t(c)]==0)) {
     cluster[t(c)]=n; Mt[t(c)]=n; tic++;
     if (nb(n(c))) {traceLake(l(c),n);};  
     if (nb(p(c))) {traceLake(r(c),n);}; 
     };
  };   

void makeClusters() {
  nclusters=0;
  int bestTic=0;
  for (int t=0; t<nt; t++) {
     int c=t*3; 
     if ( (skeleton[t]==2) && (cluster[t]==0) ) {
       tic=0;
       nclusters++;   //   println("tracing cluster "+nclusters+" from c="+c);
       traceCluster(t*3,nclusters); cluster[t(c)]=0;
       traceCluster(t*3+1,nclusters); cluster[t(c)]=0;
       traceCluster(t*3+2,nclusters); 
       if (tic>bestTic) {bestTic=tic; bestCluster=nclusters;};
       };
     };
   }  
   
void mergeCluster(int other, int k) {
  for (int t=0; t<nt; t++) {
      if ( (skeleton[t]==2) && (cluster[t]==other) ) {cluster[t]=k;};
      if ( skeleton[t]==1 ) {cluster[t]=k; skeleton[t]=2;};
     };
   }    

void traceCluster(int c, int n) {
  if ((skeleton[t(c)]==2) && (cluster[t(c)]==0)) {
     cluster[t(c)]=n; tic++;
     if (nb(n(c))) {traceCluster(l(c),n);};  
     if (nb(p(c))) {traceCluster(r(c),n);}; 
      };
  };
  
int  growCluster(int k) {
   int tc=0; 
   for(int i=0; i<nc; i++) {P[i]=false;};
   for(int t=0; t<nt; t++) {
     Mt[t]=0; 
     if (skeleton[t]==0) { Mt[t]=0; tc++;};  
     if ( (skeleton[t]==2) && (cluster[t]==k) ) {Mt[t]=1;};
     }; 
   int r=1;
   int found=0;
   boolean reached=false;
   while ((tc>0) && (r<50)) {
     for(int i=0; i<nc; i++) {
       if ((skeleton[t(i)]==0)&&(Mt[t(i)]==0)) {       
           if(nb(i)&&(Mt[t(o(i))]==r)) {
              Mt[t(i)]=r+1; tc--; P[i]=true; 
              if (!reached) {
                 if ( (skeleton[t(l(i))]==2) && (cluster[t(l(i))]!=k) ) {c=i; reached=true; found=cluster[t(l(i))];};
                 if ( (skeleton[t(r(i))]==2) && (cluster[t(r(i))]!=k) ) {c=i; reached=true; found=cluster[t(r(i))];};
                 };
              };
           };
       };
     r++;
     };
   rings=r+1; 
   return(found);
   };
   
void  traceBridge(int k) {
   rings=1;
   skeleton[t(c)]=1;
   while (cluster[t(c)]!=k) {
     rings++;
     if (P[c]) {c=o(c);} else {if (P[p(c)]) {c=r(c);} else {c=l(c);};}; 
     Mt[t(c)]=rings; skeleton[t(c)]=1;
     };
  };