// GLOBAL VARIABLES
Mesh M = new Mesh();
Mesh S = new Mesh();
int showSkeletOnly=1;
int rings=1;                               // number of rings for colorcoding
float r=10.0;                                // radius of spheres for displaying vertices

class Mesh {
// int w=6;                                   // size of initial Tmesh-grid made in setup

//  ==================================== INIT, CREATE, COPY MESH ====================================
 Mesh() {}
 void reset() {nv=0; nc=0; nt=0;}
 void declare() {
   for (int i=0; i<maxnv; i++) {G[i]=new pt(0,0,0); Nv[i]=new vec(0,0,0);};   // init vertices and normals
   for (int i=0; i<maxnt; i++) {Nt[i]=new vec(0,0,0);  };       // init triangle normals and skeleton lab els
   }
 void makeGrid (int w) { // make a 2D grid of vertices
  for (int i=0; i<w; i++) {for (int j=0; j<w; j++) { G[w*i+j].setTo(height*.8*j/(w-1)+height/10,height*.8*i/(w-1)+height/10,0);}}    
  for (int i=0; i<w-1; i++) {for (int j=0; j<w-1; j++) {                  // define the triangles for the grid
    V[(i*(w-1)+j)*6]=i*w+j;       V[(i*(w-1)+j)*6+2]=(i+1)*w+j;       V[(i*(w-1)+j)*6+1]=(i+1)*w+j+1;
    V[(i*(w-1)+j)*6+3]=i*w+j;     V[(i*(w-1)+j)*6+5]=(i+1)*w+j+1;     V[(i*(w-1)+j)*6+4]=i*w+j+1;}; };
  nv = w*w;
  nt = 2*(w-1)*(w-1); 
  nc=3*nt;
  }
 void init() { 
   for (int i=0; i<maxnt; i++) {Nt[i]=new vec(0,0,0); visible[i]=true;};  
   computeO(); 
   computeValenceAndResetNormals(); 
   computeTriNormals(); 
   computeVertexNormals(); 
   c=0; sc=0;  
  }
//  ==========================================================  VERTICES ===========================================
 int maxnv = 20000;                         //  max number of vertices
 int nv = 0;                              // current  number of vertices
 pt[] G = new pt [maxnv];                   // geometry table (vertices)
 float[] rad = new float[maxnv];                   // geometry table (vertices)
 vec[] Nv = new vec [maxnv];                 // vertex normals or laplace vectors
 int[] Mv = new int[maxnv];                  // vertex markers
 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
void computeVertexNormals() {  // computes the vertex normals as sums of the normal vectors of incident tirangles scaled by area/2
  for (int i=0; i<nv; i++) {Nv[i].setTo(0,0,0);};  // resets the valences to 0
  for (int i=0; i<3*nt; i++) {Nv[v(i)].add(Nt[t(i)]);};
  for (int i=0; i<nv; i++) {Nv[i].makeUnit();}; 
  };
int addVertex(pt P) { G[nv].setTo(P); nv++; return nv-1;};
int addVertex(pt P, float pr) { G[nv].setTo(P); rad[nv]=pr; nv++; return nv-1;};
int addVertex(float x, float y, float z) { G[nv].x=x; G[nv].y=y; G[nv].z=z; nv++; return nv-1;};
  

//  ==========================================================  TRIANGLES ===========================================
 int maxnt = maxnv*2;                       // max number of triangles
 int nt = 0;                   // current and max number of triangles
 vec[] Nt = new vec [maxnt];                // vtriangles normals
 boolean[] visible = new boolean[maxnt];    // set if triangle visible
 int[] Mt = new int[maxnt];                 // triangle markers for distance and other things   
 boolean [] VisitedT = new boolean [maxnt];  // triangle visited
pt triCenter(int i) {return(triCenterFromPts( G[V[3*i]], G[V[3*i+1]], G[V[3*i+2]] )); };  pt triCenter() {return triCenter(t());}  // computes center of triangle t(i) 
vec triNormal(int i) { return(triNormalFromPts(G[V[3*i]], G[V[3*i+1]], G[V[3*i+2]])); };  vec triNormal() {return triNormal(t());} // computes triangle t(i) normal * area / 2
void computeTriNormals() {for (int i=0; i<nt; i++) {Nt[i].setToVec(triNormal(i)); }; };             // caches normals of all tirangles
void writeTri (int i) {println("T"+i+": V = ("+V[3*i]+":"+v(o(3*i))+","+V[3*i+1]+":"+v(o(3*i+1))+","+V[3*i+2]+":"+v(o(3*i+2))+")"); };
void hitTriangle() {
   float smallestDepth=10000000;
  boolean hit=false;
  for (int t=0; t<nt; t++) {
    if (rayHitTri(eye,mark,g(3*t),g(3*t+1),g(3*t+2))) {hit=true;
      float depth = rayDistTriPlane(eye,mark,g(3*t),g(3*t+1),g(3*t+2));
      if ((depth>0)&&(depth<smallestDepth)) {smallestDepth=depth;  c=3*t;};
      }; 
    };
  if (hit) {
    pt X = eye.make(); X.addScaledVec(smallestDepth,eye.vecTo(mark));
    mark.setToPoint(X);
    float distance=X.disTo(g(c));
    int b=c;
    if (X.disTo(g(n(c)))<distance) {b=n(c); distance=X.disTo(g(b)); };
    if (X.disTo(g(p(c)))<distance) {b=p(c);};
    c=b;
    };
  }
void addTriangle(int i, int j, int k) {V[nc++]=i; V[nc++]=j; V[nc++]=k; visible[nt]=true; nt++;}
  
  
// ============================================= ISOCURVES =======================================
void makeIsosurface (float rr, Mesh S) {
  for (int cc=0; cc<nc; cc++) marked[cc]=false;
  S.reset();
  for (int cc=0; cc<nc; cc++) if (!marked[cc]&&isStick(cc,rr) ) {
    pt F = centerEdgesIsoloop(cc,rr); 
    makeFan(cc,rr,F,S);
    };
  }

void makeFan(int cc, float rr, pt F, Mesh S) {
    int f=0, l=0, n=0;
    f=S.addVertex(F);
    if (rad[v(n(cc))]<rr) cc=o(cc);
    marked[cc]=true; marked[o(cc)]=true;
    pt L = isopointOnStick(cc,rr); l=S.addVertex(L); cc=nextOnIsoloop(cc,rr); 
    int l1=l;
    int t=t(cc);
    while (t(o(cc))!=t)  {
       marked[cc]=true; marked[o(cc)]=true; 
       pt N = isopointOnStick(cc,rr); n=S.addVertex(N);  
       S.addTriangle(f,l,n);
       L.setTo(N); l=n;
       cc=nextOnIsoloop(cc,rr); } ;
       S.addTriangle(f,n,l1);
//    marked[cc]=true; marked[o(cc)]=true; 
//    pt N = isopointOnStick(cc,rr);  n=S.addVertex(N);  S.addTriangle(f,l,n);
 //   L.setTo(N); l=n;  cc=nextOnIsoloop(cc,rr);
   }
 
   
void fanIsocurves (float rr) {
  for (int cc=0; cc<nc; cc++) marked[cc]=false;
  S.reset();
  for (int cc=0; cc<nc; cc++) if (!marked[cc]&&isStick(cc,rr) ) {
    pt F = centerEdgesIsoloop(cc,rr); 
    fanIsoloop(cc,rr,F);
    };
  }

pt centerEdgesIsoloop(int cc, float rr) {
    if (rad[v(n(cc))]<rr) cc=o(cc);
    pt C=new pt(); float D=0;
    pt L = isopointOnStick(cc,rr); cc=nextOnIsoloop(cc,rr); 
    int t=t(cc);
    while (t(o(cc))!=t)  {pt N = isopointOnStick(cc,rr);  pt M = midPt(L,N); float d=L.disTo(N); M.mul(d); D+=d; C.addPt(M); L.setTo(N); cc=nextOnIsoloop(cc,rr); } ;
    pt N = isopointOnStick(cc,rr);  pt M = midPt(L,N); float d=L.disTo(N); M.mul(d); D+=d; C.addPt(M); L.setTo(N); cc=nextOnIsoloop(cc,rr);
    C.mul(1./D);
    return(C);
   }

int nextOnIsoloop(int cc, float rr) {int ncc=c=n(o(cc)); if(rad[v(o(cc))]>rr) ncc=n(ncc); return ncc;}

boolean isStick(int cc, float v) {  float a = rad[v(n(cc))];   float b = rad[v(p(cc))];  return ((b<v)&&(v<a)) ; }

void showIsopointsOnSticks (float rr) {
  for (int cc=0; cc<nc; cc++) if (isStick(cc,rr)) isopointOnStick(cc, rr).show(6);
  }

void traceIsocurves (float rr) {
  for (int cc=0; cc<nc; cc++) marked[cc]=false;
  for (int cc=0; cc<nc; cc++) if (!marked[cc]&&isStick(cc,rr) ) {
    traceIsoloop(cc,rr);
    pt F = centerEdgesIsoloop(cc,rr); F.show(5);
   };
  }
  
void fanIsoloop(int cc, float rr, pt F) {
    if (rad[v(n(cc))]<rr) cc=o(cc);
    marked[cc]=true; marked[o(cc)]=true;
    pt L = isopointOnStick(cc,rr); cc=nextOnIsoloop(cc,rr); 
    int t=t(cc);
    while (t(o(cc))!=t)  {
       marked[cc]=true; marked[o(cc)]=true; 
       pt N = isopointOnStick(cc,rr);   
       showTriangle(F,L,N);
       L.setTo(N); 
       cc=nextOnIsoloop(cc,rr); } ;
    marked[cc]=true; marked[o(cc)]=true; 
    pt N = isopointOnStick(cc,rr);  
    showTriangle(F,L,N);  
    L.setTo(N); 
    cc=nextOnIsoloop(cc,rr);
   }
   
void traceIsoloop(int cc, float rr) {
    marked[cc]=true; marked[o(cc)]=true;
    if (rad[v(n(cc))]<rr) cc=o(cc);
    beginShape();
    marked[cc]=true; marked[o(cc)]=true; isopointOnStick(cc,rr).vert(); cc=nextOnIsoloop(cc,rr); 
    int t=t(cc);
    while (t(o(cc))!=t)  { marked[cc]=true; marked[o(cc)]=true; isopointOnStick(cc,rr).vert(); cc=nextOnIsoloop(cc,rr); } ;
    endShape(CLOSE);
   }

pt centerIsoloop(int cc, float rr) {
    if (rad[v(n(cc))]<rr) cc=o(cc);
    pt C=new pt(); int n=0;
    C.addPt(isopointOnStick(cc,rr)); n++; cc=nextOnIsoloop(cc,rr); 
    int t=t(cc);
    while (t(o(cc))!=t)  {C.addPt(isopointOnStick(cc,rr));  n++; cc=nextOnIsoloop(cc,rr); } ;
    C.mul(1./n);
    return(C);
   }

pt isopointOnStick(int cc, float v) {  float a = rad[v(n(cc))];   float b = rad[v(p(cc))];  pt A = g(n(cc));   pt B = g(p(cc)); 
  pt P = A.make();  P.moveTowards((v-a)/(b-a),B);  return P;   }; 

void showIsopointOnStick(pt A, float a, float v, pt B, float b) {
  if ((min(a,b)<v)&&(v<max(a,b))) {
   pt P = A.make();  P.moveTowards((v-a)/(b-a),B);  P.show(6); 
  }; 
}
// ============================================= DISPLAY =======================================
pt Cbox = new pt(width/2,height/2,0);                   // mini-max box center
float Rbox=1000;                                        // Radius of enclosing ball
boolean showEdges=false;
boolean showNormals=false;
boolean showVertices=false;

void computeBox() {
  pt Lbox =  G[0].make();  pt Hbox =  G[0].make();
  for (int i=1; i<nv; i++) { 
    Lbox.x=min(Lbox.x,G[i].x); Lbox.y=min(Lbox.y,G[i].y); Lbox.z=min(Lbox.z,G[i].z);
    Hbox.x=max(Hbox.x,G[i].x); Hbox.y=max(Hbox.y,G[i].y); Hbox.z=max(Hbox.z,G[i].z); 
    };
  Cbox.setToPoint(midPt(Lbox,Hbox));  Rbox=Cbox.disTo(Hbox);
//  println("Box r="+Rbox); print("C="); Cbox.write(); print("L="); Lbox.write(); print("H="); Hbox.write();
  };
pt cg(int c) {pt cPt = midPt(g(c),midPt(g(c),triCenter(t(c))));  return(cPt); };   // computes point at corner
void showCorner(int c, int r) {pt cPt = midPt(g(c),midPt(g(c),triCenter(t(c))));  cPt.show(r); };   // renders corner c as small ball
void showCornerAndNormal(int c, int r) {pt cPt = midPt(g(c),midPt(g(c),triCenter(t(c))));  noStroke(); cPt.show(r); 
                     stroke(magenta); vec N = Nt[t(c)].make(); N.makeUnit();  N.mul(10*r);  N.show(cPt);};   // renders corner c as small ball
void drawEdge(int c) {line(g(p(c)).x,g(p(c)).y,g(p(c)).z,g(n(c)).x,g(n(c)).y,g(n(c)).z); };  // draws edge of t(c) opposite to corner c
void showBorder() {for (int i=0; i<nc; i++) {if (visible[t(i)]&&border(i)) {drawEdge(i);}; }; };         // draws all border edges
void shade(int i) {if(visible[i]) { beginShape(TRIANGLES); G[V[3*i]].vert(); G[V[3*i+1]].vert(); G[V[3*i+2]].vert(); endShape(); };}; // shade tris
void showTriNormals() {vec N= new vec(0,0,0); for (int i=0; i<nt; i++) { 
N.setToVec(Nt[i]); 
N.makeUnit();  
N.mul(r);  
N.show(triCenter(i)); };  };

void showVertexNormals() {vec N= new vec(0,0,0); for (int i=0; i<nv; i++) {N.setToVec(Nv[i]); N.makeUnit(); N.mul(5*r);  N.show(G[i]); };  };
void show() {
  fill(magenta); stroke(magenta); showCornerAndNormal(c,5);  
  if (showEdges) {stroke(red); for(int i=0; i<nc; i++) if(visible[t(i)]) drawEdge(i);};  
  stroke(dblue); showBorder();
  if (showVertices) {stroke(red); fill(red); for (int v=0; v<nv; v++) G[v].show(r); };
 // fill(cyan,80); 
  }

void showTriangles() {for(int t=0; t<nt; t++) shade(t); }
void showVertices(float r) {for (int v=0; v<nv; v++) G[v].show(r);  }
void showEdges() {for(int i=0; i<nc; i++) if(visible[t(i)]) drawEdge(i);}

// ============================================= CORNER OPERATORS =======================================
 int nc = nt*3;                             // current number of corners (3 per triangle)
 int c = 0;                                 // current corner shown in image and manipulated with keys: n, p, o, l, r
 int sc=0;                                  // saved value of c
 int[] V = new int [3*maxnt];               // V table (triangle/vertex indices)
 int[] O = new int [3*maxnt];               // O table (opposite corner indices)
 int[] W = new int [3*maxnt];               // mid-edge vertex indices for subdivision (associated with corner opposite to edge)
 boolean[] marked = new boolean[3*maxnt];               // corner visited

int t (int c) {int r=int(c/3); return(r);};            int t() {return t(c);}
int n (int c) {int r=3*int(c/3)+(c+1)%3; return(r);};  int n() {return n(c);}
int p (int c) {int r=3*int(c/3)+(c+2)%3; return(r);};  int p() {return p(c);}
int v (int c) {return(V[c]);};                         int v() {return v(c);}
int o (int c) {return(O[c]);};                         int o() {return o(c);}
int l (int c) {return(o(n(c)));};                      int l() {return l(c);}
int r (int c) {return(o(p(c)));};                      int r() {return r(c);}
pt g (int c) {return(G[V[c]]);}; pt g() {return g(c);}            // shortcut to get the point of the vertex v(c) of corner c
int w (int c) {return(W[c]);};               // temporary indices to mid-edge vertices associated with corners during subdivision
boolean nb(int c) {return(O[c]!=-1);};  boolean nb() {return nb(c);}     // not border
boolean border (int c) {return(O[c]==-1);};  // returns true if corner has no opposite
  
void previous() {c=p();};
void next() {c=n();};
void opposite() {if(nb()) {c=o(c);};};
void left() {next(); opposite();};
void right() {previous(); opposite();};
void back() {opposite();};
void turn() {left(); next(); };

void writeCorner (int c) {println("c="+c+", n="+n(c)+", p="+p(c)+", o="+o(c)+", v="+v(c)+", t="+t(c)+", nt="+nt+", nv="+nv ); }; 
void writeCorner () {writeCorner (c);}
void writeCorners () {for (int c=0; c<nc; c++) {println("T["+c+"]="+t(c)+", visible="+visible[t(c)]+", v="+v(c)+",  o="+o(c));};}

// ============================================= O TABLE CONSTRUCTION =========================================
void computeOnaive() {                         // sets the O table from the V table, assumes consistent orientation of triangles
  for (int i=0; i<3*nt; i++) {O[i]=-1;};  // init O table to -1: has no opposite (i.e. is a border corner)
  for (int i=0; i<3*nt; i++) {  for (int j=i+1; j<3*nt; j++) {       // for each corner i, for each other corner j
      if( (v(n(i))==v(p(j))) && (v(p(i))==v(n(j))) ) {O[i]=j; O[j]=i;};};}; // make i and j opposite if they match         
  };

void computeO() { 
  int nIC [] = new int [maxnv];                   // number of incident corners
  int maxValence=0;
  for (int c=0; c<nc; c++) {O[c]=-1;};  // init O table to -1: has no opposite (i.e. is a border corner)
  for (int v=0; v<nv; v++) {nIC[v]=0; };
  for (int c=0; c<nc; c++) {nIC[v(c)]++;}
  for (int v=0; v<nv; v++) {if(nIC[v]>maxValence) {maxValence=nIC[v]; };};
//  println(" Max valence = "+maxValence+". ");
  int IC [][] = new int [maxnv][maxValence];                   // incident corners
  for (int v=0; v<nv; v++) {nIC[v]=0; };
  for (int c=0; c<nc; c++) {IC[v(c)][nIC[v(c)]++]=c;}
  for (int c=0; c<nc; c++) {
    for (int i=0; i<nIC[v(p(c))]; i++) {
      int a = IC[v(p(c))][i];
      for (int j=0; j<nIC[v(n(c))]; j++) {
         int b = IC[v(n(c))][j];
         if ((b==n(a))&&(c!=n(b))) {O[c]=n(b); O[n(b)]=c; };
         };
      };
    };
  }


//  ==========================================================  SIMPLIFICATION ===========================================

void flip() {flip(c);}
void flip(int c) {      // fip edge opposite to corner c
    V[n(o(c))]=v(c); V[n(c)]=v(o(c));
    int co=o(c); O[co]=r(c); O[r(c)]=co; O[c]=r(co); O[r(co)]=c; O[p(c)]=p(co); O[p(co)]=p(c);
  }
  
void doFlips() {  for (int c=0; c<3*nt; c++) {
  if (nb(c)) {if (g(n(c)).disTo(g(p(c)))>g(c).disTo(g(o(c)))) {flip(c);}; }; };
  } // assumes manifold

void collapse() {collapse(c);}
void collapse(int c) {      // collapse edge opposite to corner c
   int b=p(c), oc=o(c), vnc=v(n(c));
   visible[t(c)]=false; visible[t(oc)]=false;
   if (true) return;
   for (int a=b; a!=n(oc); a=p(r(a))) {V[a]=vnc;}; V[p(c)]=vnc; V[n(oc)]=vnc; 
    O[l(c)]=r(c); O[r(c)]=l(c);     O[l(oc)]=r(oc); O[r(oc)]=l(oc); 
  }

//  ==========================================================  HOLES ===========================================
pt holeCenter = new pt (0,0,0);
vec holeNormal = new vec(0,0,1);
pt  centerOfHole() {pt C=new pt(0,0,0); int nb=0; for (int i=0; i<nc; i++) {if (visible[t(i)]&&border(i)) {nb++; C.addPt(g(p(i)));}; }; C.mul(1./nb); return C;};         // draws all border edges
vec  normalOfHole(pt C) {vec N=new vec(0,0,0); for (int i=0; i<nc; i++) {if (visible[t(i)]&&border(i)) N.add(cross(C.vecTo(g(p(i))),C.vecTo(g(n(i))))); }; N.makeUnit(); return N;};         // draws all border edges
void excludeInvisibleTriangles () {for (int b=0; b<nc; b++) {if (!visible[t(o(b))]) {O[b]=-1;};};}
void hole() {holeCenter.setTo(centerOfHole()); holeNormal.setTo(normalOfHole(holeCenter)); };
void fanHoles() {
  println("FANHOLES: nv="+nv +", nt="+nt +", nc="+nc );
  for (int t=0; t<nt; t++) {VisitedT[t]=false;};
  int lnt=nt;
  int L=0;
  for (int cc=0; cc<nc; cc++) {
   if (visible[t(cc)]&&(!VisitedT[t(cc)]) && (!nb(cc) )) {L++;
      print("<"); G[nv].setTo(0,0,0); int hl=fanHole(cc,L); G[nv].mul(1.0/float(hl)); nv++; println("> hl="+hl);
     };
   };
   for (int t=lnt; t<nt; t++) {visible[t]=true;};
   nc=3*nt;
   println("Filled "+L+" holes");
  }
  
int fanHole(int cc, int L) {
 int hl=0; int o=0;  int f=cc;
 VisitedT[t(f)]=true; hl++; o=3*nt; V[o]=nv; V[n(o)]=v(p(f)); V[p(o)]=v(n(f)); O[o]=f; O[f]=o; nt++; G[nv].addPt(g(p(f)));
 int lc=p(o); 
 f=n(f);  while(nb(f)) {f=n(o(f)); }; 
 while (f!=cc) {    
       print("."); VisitedT[t(f)]=true; hl++; o=3*nt; V[o]=nv; V[n(o)]=v(p(f)); V[p(o)]=v(n(f)); O[o]=f; O[f]=o; nt++; G[nv].addPt(g(p(f)));
      O[n(o)]=lc; O[lc]=n(o); lc=p(o);
      f=n(f);  while(nb(f)&&(f!=cc)) {f=n(o(f)); }; 
      }; 
  O[lc]=n(o(cc)); O[n(o(cc))]=lc; c=lc;
  return(hl);
  }
  
void compactVO() {  
  println("COMPACT TRIANGLES: nv="+nv +", nt="+nt +", nc="+nc );
  int[] U = new int [nc];
  int lc=-1; for (int c=0; c<nc; c++) {if (visible[t(c)]) {U[c]=++lc; }; };
  for (int c=0; c<nc; c++) {if (nb(c)) {O[c]=U[o(c)];} else {O[c]=-1;}; };
  int lt=0;
  for (int t=0; t<nt; t++) {
    if (visible[t]) {
      V[3*lt]=V[3*t]; V[3*lt+1]=V[3*t+1]; V[3*lt+2]=V[3*t+2]; 
      O[3*lt]=O[3*t]; O[3*lt+1]=O[3*t+1]; O[3*lt+2]=O[3*t+2]; 

      visible[lt]=true; 
      lt++;
      };
    };
nt=lt; nc=3*nt;    
  println("      ...  NOW: nv="+nv +", nt="+nt +", nc="+nc );
  }

void compactV() {  
  println("COMPACT VERTICES: nv="+nv +", nt="+nt +", nc="+nc );
  int[] U = new int [nv];
  boolean[] deleted = new boolean [nv];
  for (int v=0; v<nv; v++) {deleted[v]=true;};
  for (int c=0; c<nc; c++) {deleted[v(c)]=false;};
  int lv=-1; for (int v=0; v<nv; v++) {if (!deleted[v]) {U[v]=++lv; }; };
  for (int c=0; c<nc; c++) {V[c]=U[v(c)]; };
  lv=0;
  for (int v=0; v<nv; v++) {
    if (!deleted[v]) {G[lv].setToPoint(G[v]);  deleted[lv]=false; 
      lv++;
      };
    };
 nv=lv;
 println("      ...  NOW: nv="+nv +", nt="+nt +", nc="+nc );
  }

// ============================================================= SMOOTHING ============================================================
void computeValenceAndResetNormals() {      // caches valence of each vertex
  for (int i=0; i<nv; i++) {Nv[i].setTo(0,0,0); Valence[i]=0;};  // resets the valences to 0
  for (int i=0; i<nc; i++) {Valence[v(i)]++; };
  }

void computeLaplaceVectors() {  // computes the vertex normals as sums of the normal vectors of incident tirangles scaled by area/2
  computeValenceAndResetNormals();
  for (int i=0; i<3*nt; i++) {Nv[v(p(i))].add(g(p(i)).vecTo(g(n(i))));};
  for (int i=0; i<nv; i++) {Nv[i].div(Valence[i]);}; 
  };
  
void tuck(float s) {for (int i=0; i<nv; i++) {if (!Border[i]) G[i].addScaledVec(s,Nv[i]);}; };  // displaces each vertex by a fraction s of its normal

void markBorderVertices() {
  for (int i=0; i<nv; i++) {Border[i]=false;}; 
  for (int i=0; i<nc; i++) {if (border(i)) Border[v(n(i))]=true;};
 }

// ============================================================= SUBDIVISION ============================================================
void splitEdges() {            // creates a new vertex for each edge and stores its ID in the W of the corner (and of its opposite if any)
  for (int i=0; i<3*nt; i++) {  // for each corner i
    if(border(i)) {G[nv]=midPt(g(n(i)),g(p(i))); W[i]=nv++;}
    else {if(i<o(i)) {G[nv]=midPt(g(n(i)),g(p(i))); W[o(i)]=nv; W[i]=nv++; }; };  // if this corner is the first to see the edge
    };
  };
  
void bulge() {              // tweaks the new mid-edge vertices according to the Butterfly mask
  for (int i=0; i<3*nt; i++) {
    if((nb(i))&&(i<o(i))) {    // no tweak for mid-vertices of border edges
     if (nb(p(i))&&nb(n(i))&&nb(p(o(i)))&&nb(n(o(i))))
      {G[W[i]].addScaledVec(0.25,midPt(midPt(g(l(i)),g(r(i))),midPt(g(l(o(i))),g(r(o(i))))).vecTo(midPt(g(i),g(o(i))))); };
      }; 
    };
  };
  
void splitTriangles() {    // splits each tirangle into 4
  for (int i=0; i<3*nt; i=i+3) {
    V[3*nt+i]=v(i); V[n(3*nt+i)]=w(p(i)); V[p(3*nt+i)]=w(n(i));
    V[6*nt+i]=v(n(i)); V[n(6*nt+i)]=w(i); V[p(6*nt+i)]=w(p(i));
    V[9*nt+i]=v(p(i)); V[n(9*nt+i)]=w(n(i)); V[p(9*nt+i)]=w(i);
    V[i]=w(i); V[n(i)]=w(n(i)); V[p(i)]=w(p(i));
    };
  nt=4*nt; nc=3*nt;
  };

// ********************* FIT TO ISOSURFACE IN C *****************
void fit(Cube C) {
  computeVertexNormals();
  for (int i=0; i<nv; i++) {G[i].addScaledVec((C.value(G[i])-threshold)/5.0,Nv[i]);}; 
  }

  } // ==== END OF MESH CLASS