int[] Mv = new int[maxnv];                  // vertex markers
int rings=1;                               // number of rings for colorcoding
int[] Tt = new int[maxnt];                 // triangle type
int[] Tc = new int[3*maxnt];               // corner type
int[] Mt = new int[maxnt];                 // triangle markers for distance and other things   

/****/
int [] Valence = new int [maxnv];          // vertex valence (count of incident triangles)
boolean [] Border = new boolean [maxnv];   // vertex is border
boolean [] VisitedV = new boolean [maxnv];  // vertex visited
boolean [] VisitedT = new boolean [maxnv];  // vertex visited
boolean[] Land = new boolean[maxnt];       // triangle type land (and bridge) or water
int[] Lake = new int[maxnt];               // Which lake , for water triangles   
int lake=1;                                // currently selected lake >1;
int nLakes=1;                              // number of lakes
int lakeSize;                              // triangle count in current lake
int[] Island = new int[maxnt];             // triangle markers for distance and other things   
int nIslands=1;                            // number of islands
int islandSize;                            // triangle count in current island
int largestIsland;
int sizeOfLargestIsland;
int[] Distance = new int[maxnt];           // triangle markers for distance fields 
int[] Tv = new int[maxnv];                 // vertex type: 0=water, 1=border, 2=land
boolean[] P = new boolean [3*maxnt];       // marker of corners in a path to parent triangle
int c = 0;                                 // current corner shown in image and manipulated with keys: n, p, o, l, r
int sc=0;
int maxValence=0; 
int longestShoreLake=0;
int longestShore=0;

int step=1;
void findLakes() {
if (step==1) {print("Initialize valences for "+nv+" vertices");
 for (int i=0; i<nv; i++) {Valence[i]=0;};  printVertices();
 for (int c=0; c<nc; c++) {Valence[v(c)]++; } setVertexTypesToValence(); printVertices();
 maxValence=0;  for (int i=0; i<nv; i++) {if (Valence[i]>maxValence) { maxValence=Valence[i];};}; println("  Maximum valences = "+maxValence);
};if (step==2) { println("Initialize Land and Distance for "+nt+" triangles");
for(int t=0; t<nt; t++) {Distance[t]=0; Land[t]=false;};  
};if (step==3) { println("Make islands around singular vertices");
for(int c=0; c<nc; c++) {if (Valence[v(c)]!=6) {Distance[t(c)]=1; Land[t(c)]=true;};}; 
};if (step==4) { println("Dry coves (water triangles with 2 or 3 land neighbors)");
dryCoves(); 
};if (step==5) { print("Conquer Islands: ");
for(int t=0; t<nt; t++) {Island[t]=0;}; 
conquerIslands();
if(nIslands==0) {print("                  Found no island! "); nIslands=1;}
else {println(" Found "+nIslands+" islands. Island "+largestIsland+" has "+sizeOfLargestIsland+" triangles"); };
printSkeleton(); setVertexTypes(); printVertexTypes();
};if (step==6) { println("Merge islands "); 
  for (int i=1; i<nIslands; i++) {
    int other=expandLargestIsland(); c=sc; 
     backTraceBridge(sc, other); 
    invadeIsland(other); 
    };
printSkeleton(); setVertexTypes(); printVertexTypes();
};if (step==7) { println("Dry and merge coves"); 
dryMergeCoves(); printSkeleton();
};if (step==8) {println("Dry and merge more coves"); 
dryMergeCoves();dryMergeCoves();dryMergeCoves();dryMergeCoves(); 
 setVertexTypes(); printVertexTypes(); printSkeleton();
};if (step==9) { println("Conquer lakes "); 
conquerLakes(); println("                 Found "+nLakes+" Lakes"); 
};if (step==10) { println("Measure lakes "); 
setVertexTypes(); printVertexTypes(); measureLakes(); 
};if (step==11) { println("Fill Small Lakes ");
fillSmallLakes();  setVertexTypes(); printVertexTypes(); printSkeleton(); measureLakes();
};if (step==12) { println("fan lakes");
fanLakes(); 
};if (step==13) { println("Collapse all lakes"); 
excludeInvisibleTriangles(); fanHoles(); compactVO(); compactV(); setVertexTypes(); printVertexTypes();
//  for (int i=1; i<=nLakes; i++) {println("                ...  collapsing lake "+i+" out of "+nLakes); collapseLake(i);}; 
printSkeleton();
};if (step==14) { println("Measure lakes ");
measureLakes(); c=sc; printSkeleton();
println("Set c on Lake "+longestShoreLake+" with shore length = "+longestShore); 
};if (step==15) { println("DONE ");

}; 
  }
  
void dryCoves() {for(int c=0; c<nc; c++) {if (Land[t(l(c))] && Land[t(r(c))]) { Land[t(c)]=true;}; }; } 
void dryMergeCoves() {for(int c=0; c<nc; c++) {if (Land[t(l(c))] && Land[t(r(c))]) { Land[t(c)]=true; Island[t(c)]=largestIsland;}; }; };

void collapseLake(int L) {
  println("       > collapsing lake "+L);
   pt LakeCenter= new pt(0,0,0);
   for (int v=0; v<nv; v++) {VisitedV[v]=false;};
   int interior=0, border=0;
   for (int c=0; c<nc; c++) {if ( !Land[t(c)] && (Lake[t(c)]==L) && !VisitedV[v(c)] ) {
     if (Tv[v(c)]==0) {interior++; LakeCenter.addPt(g(c)); } 
     else {border++;}; VisitedV[v(c)]=true;}; };
   LakeCenter.mul(1.0/interior);
   print("      collapsing "+interior+" interior points to "); LakeCenter.write();
   for (int c=0; c<nc; c++) {if ( !Land[t(c)] && (Lake[t(c)]==L) ) {if (Tv[v(c)]==0) { g(c).setToPoint(LakeCenter); /***/ }; }; };
  G[nv++].setToPoint(LakeCenter);
  boolean found=false;
  int cc=1;
  for (cc=0; cc<nc; cc++) {
   if ( !found  && (Land[t(cc)]) && nb(cc) && (Lake[t(o(cc))]==L)   && (!Land[t(o(cc))]) ) {found=true; c=cc;};};
   int f=n(c); Island[t(f)]=0;
   boolean fin=false;
   int lastf=c;
   println("      found starting corner "+f);
   while (!fin) {    
      while(nb(f)&&Land[t(o(f))]) {f=n(o(f)); }; 
      Island[t(f)]=0;  V[o(f)]=nv-1;        // *****
      O[n(o(f))]=p(o(lastf)); O[p(o(lastf))]=n(o(f));
      lastf=f;
      if (f==c) {fin=true;};
      f=n(f);
      }; 
    f=o(c); 
    vec bulge=triNormalFromPts( g(f) , g(n(f)) , g(p(f)) );
    f=n(l(f));
 //   print("adding bulge ... "); sc=c;
 //   while (o(f)!=c) {bulge.add(triNormalFromPts( g(f) , g(n(f)) , g(p(f)) )); f=n(l(f));};
 //   println(" done"); 
    for (int c=0; c<nc; c++) {if ( !Land[t(c)] && (Lake[t(c)]==L) ) {if (Tv[v(c)]==0) {g(c).setToPoint(G[nv-1]); }; }; };
    if (L==longestShoreLake) {sc=p(o(c)); println("setting sc to "+sc+" for longest shore lake "+L); };
    };  
  
void fillSmallLakes() {
  int interior, border;
  int remainingLakes=0;
  for (int v=0; v<nv; v++) {VisitedV[v]=false;};
  println(nLakes+" lakes:");
  for (int i=1; i<=nLakes; i++) {
    interior=0; border=0; 
    for (int c=0; c<nc; c++) {if ( !Land[t(c)] && (Lake[t(c)]==i) && !VisitedV[v(c)] ) {if (Tv[v(c)]==0) {interior++;} else {border++;}; VisitedV[v(c)]=true;}; };
    if (border<=maxValence) {fillLake(i);  println("          > filled lake "+i+" with "+interior+" interior & "+border+" border vertices"); }
    else {remainingLakes++;    for (int t=0; t<nt; t++) {if (!Land[t] && (Lake[t]==i)) { Lake[t]=remainingLakes;}; };};    
    };   
  nLakes=remainingLakes; 
  };

void fillLake(int L) {
  for (int t=0; t<nt; t++) {if (!Land[t] && (Lake[t]==L)) {Land[t]=true; Island[t]=largestIsland;} ;};  
  };
  
void measureLakes() {
  int interior, border;
  for (int v=0; v<nv; v++) {VisitedV[v]=false;};
  println(nLakes+" lakes:");
  longestShore=0;longestShoreLake=0;
  for (int i=1; i<=nLakes; i++) {
    interior=0; border=0; 
    for (int c=0; c<nc; c++) {if ( !Land[t(c)] && (Lake[t(c)]==i) && !VisitedV[v(c)] ) {if (Tv[v(c)]==0) {interior++;} else {border++;}; VisitedV[v(c)]=true;}; };
     println("          > "+interior+" interior & "+border+" border vertices in lake "+i);   
     if (border>longestShore) {longestShore=border; longestShoreLake=i;};
     }; 
  println("          > Lake "+longestShoreLake+" has the longest shore of "+longestShore); 
  for (int c=0; c<nc; c++) {if (!Land[t(c)] && (Lake[t(c)]==longestShoreLake) && (Lake[t(r(c))]==longestShoreLake) && (Lake[t(l(c))]!=longestShoreLake) ) {sc=c;};};
  };
  
void traceLakeBorder(int L) {
  boolean found=false;
  int cc=1;
  for (cc=0; cc<nc; cc++) {
   if ( !found  && (!Land[t(cc)]) && (Lake[t(cc)]==L) && nb(cc)  && (Land[t(o(cc))]) ) {found=true; c=cc;};};
   int f=n(c);
   cc=1;
   print("       * lake "+L+" = (");
   boolean fin=false;
   while (!fin) {    
      cc=1;  while(nb(f)&&!Land[t(o(f))]) {f=n(o(f)); cc++;};  print(","+cc);
      if (f==c) {fin=true;};
      f=n(f);
      }; /**/
  println(")");
  }

/**** STACK OVERLFOW ****/
void conquerLakes() {
  nLakes=0;
  lakeSize=0;
//  for(int t=0; t<nt; t++) {if (!Land[t]) {Lake[t]=0;} else {lakeSize++;}; };  
  for(int t=0; t<nt; t++) {if (!Land[t]) {Lake[t]=0;lakeSize++;}; };  
  println("The mesh has "+lakeSize+" water triangles separated into lakes:"); 
  for (int t=0; t<nt; t++) {
     int c=t*3; 
     if ( !Land[t] && (Lake[t]==0) ) {
       lakeSize=1;  nLakes++;   
       if (nb(c)) {conquerLake(o(c));};  if (nb(n(c))) {conquerLake(l(c));};  if (nb(p(c))) {conquerLake(r(c));}; 
       println("  Lake "+nclusters+" has "+lakeSize+" triangles");
       };
     };
   }  
   
void conquerLake(int c) {
  if (!Land[t(c)] && (Lake[t(c)]==0)) {
     Lake[t(c)]=nLakes; lakeSize++;
     if (nb(n(c))) {conquerLake(l(c));};  if (nb(p(c))) {conquerLake(r(c));}; 
     };
  };   
  
void conquerIslands() {
  nIslands=0;
  sizeOfLargestIsland=0;
  largestIsland=0;
  for (int t=0; t<nt; t++) {
    if (Land[t]&&Island[t]==0) {c=3*t;
       nIslands++; Island[t(c)]=nIslands; islandSize=1;
       if (nb(c)) {conquerIsland(o(c));};  if (nb(n(c))) {conquerIsland(l(c));};  if (nb(p(c))) {conquerIsland(r(c));}; 
       if (islandSize>sizeOfLargestIsland) {largestIsland=nIslands; sizeOfLargestIsland=islandSize;};
       };
     };
  }


void conquerIsland (int c) {
  if ((Land[t(c)]) && (Island[t(c)]==0)) {
     Island[t(c)]=nIslands; islandSize++;  
     if (nb(n(c))) {conquerIsland(l(c));};   if (nb(p(c))) {conquerIsland(r(c));}; 
     };
  };

int expandLargestIsland() {
   int tc=0;                               // counts water triangles
   for(int i=0; i<nc; i++) {P[i]=false;}; // clear path to parent
   for(int t=0; t<nt; t++) {
     Distance[t]=0;
     if (!Land[t]) {  tc++;};  
     if ( (Land[t]) && (Island[t]==largestIsland) ) {Distance[t]=1;};
     }; 
   int r=1;
   int found=0;
   boolean reached=false;
   while ((tc>0) && (r<15) && (!reached)) {
     for(int c=0; c<nc; c++) {
       if ((!Land[t(c)])&&(Distance[t(c)]==0)) {       
           if(nb(c)&&(Distance[t(o(c))]==r)) {
              Distance[t(c)]=r+1; tc--; P[c]=true; 
              if (!reached) {
        //         print("Testing "+c); println(" Land="+Land[t(l(c))]+", Island="+Island[t(l(c))]+", Distance="+Distance[t(c)]); 
                 if ( Land[t(l(c))] && (Island[t(l(c))]!=largestIsland) ) {sc=c; reached=true; found=Island[t(l(c))];  };
                 if ( Land[t(r(c))] && (Island[t(r(c))]!=largestIsland) ) {sc=c; reached=true; found=Island[t(r(c))]; };
                 };
              };
           };
       };
     r++;
     }; /**/
   rings=r+1; 
   return(found);
   };
 
void  backTraceBridge(int c, int other) {
   while (Island[t(c)]!=largestIsland) {
      Land[t(c)]=true; Island[t(c)]=other;
      if (P[c]) {c=o(c);} else {if (P[p(c)]) {c=r(c);} else {c=l(c);};}; 
      };
  };
      
void invadeIsland(int other) {for (int t=0; t<nt; t++) {if (Land[t]&&(Island[t]==other)) {Island[t]=largestIsland;};  }; };

void printSkeleton() {int ground=0, water=0; 
  for (int t=0; t<nt; t++) {if (Land[t]) {ground++;} else {water++;};}; 
  println("            > "+ground+" ground & "+water+" water triangles. Total="+nt);};

void printVertices() {int regular=0, irregular=0;
  for (int v=0; v<nv; v++) {if (Valence[v]==6) {regular++;} else {irregular++;};}; 
  println("          > "+regular+" regular & "+irregular+" irregular vertices. Total="+nv);};
  
void printVertexTypes() {int land=0, water=0, border=0;
  for (int v=0; v<nv; v++) {
    if (Tv[v]==0) {water++;};
    if (Tv[v]==1) {border++;};
    if (Tv[v]==2) {land++;};
    };
  println("          > "+water+" water, "+border+" border, "+land+" land vertices. Total="+nv);};
  
void setVertexTypes() {
  for (int c=0; c<nc; c++) {if (Land[t(c)]) {Tv[v(c)]=2;} else {Tv[v(c)]=0;};};
  for (int c=0; c<nc; c++) {if ( (Land[t(c)]&&(Tv[v(c)]==0)) ||(!Land[t(c)]&&(Tv[v(c)]==2))) {Tv[v(c)]=1;};};
  };
  
void setVertexTypesToValence() {
  for (int v=0; v<nv; v++) {if (Valence[v]==6) {Tv[v]=0;} else {Tv[v]=2;};}; 
  };