//Half-Edge Mesh
//Source: http://www.flipcode.com/archives/The_Half-Edge_Data_Structure.shtml
class HE_Mesh{
  ArrayList vertices; 
  ArrayList faces;
  ArrayList halfEdges;
  ArrayList edges;
  PVector bodyCenter;
  ArrayList faceCenters;
  ArrayList faceNormals;
  ArrayList vertexNormals;
  ArrayList edgeCenters;
  ArrayList edgeNormals;
  boolean bodyCenterUpdated;
  boolean faceCentersUpdated;
  boolean faceNormalsUpdated;
  boolean vertexNormalsUpdated;
  boolean edgeCentersUpdated;
  boolean edgeNormalsUpdated;  
  boolean indicesUpdated;
  
  boolean lastOperationFailed;


  //CONSTRUCTOR
  HE_Mesh(){
    vertices=new ArrayList(); 
    faces=new ArrayList();
    halfEdges=new ArrayList();
    edges=new ArrayList();
    faceCenters=new ArrayList();
    faceNormals=new ArrayList();
    vertexNormals=new ArrayList();
    edgeCenters=new ArrayList();
    edgeNormals=new ArrayList();
    bodyCenter=new PVector(0,0,0) ;
    bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
  }

  //DRAWING ROUTINES
  void draw(){
    if(!faceCentersUpdated) setFaceCenters();
    Iterator faceItr = faces.iterator();
    Iterator faceCtrItr = faceCenters.iterator();
    while(faceItr.hasNext()){
      HE_Face face=(HE_Face)faceItr.next();
      HE_Vertex faceCenter =(HE_Vertex)faceCtrItr.next();
      if(face.visible){
        HE_HalfEdge halfEdge=face.halfEdge;
        beginShape(TRIANGLE_STRIP);
        do{              
          HE_Vertex v=halfEdge.vert; 
          vertex(v.x,v.y,v.z);
          vertex(faceCenter.x,faceCenter.y,faceCenter.z);
          halfEdge= halfEdge.next;
        } 
        while(halfEdge!=face.halfEdge);
 HE_Vertex v=halfEdge.vert; 
vertex(v.x,v.y,v.z);
        endShape();
        //}
      }
    }
  }   


  void drawVertices(float s){
    Iterator vertexItr = vertices.iterator();
    while(vertexItr.hasNext()){
      HE_Vertex v=(HE_Vertex)vertexItr.next();
      pushMatrix();
      translate(v.x,v.y,v.z);  
      box(s);
      popMatrix();
    }
  }

  void drawFaceCenters(){
    Iterator vertexItr = faceCenters.iterator();
    while(vertexItr.hasNext()){
      HE_Vertex v=(HE_Vertex)vertexItr.next();
      pushMatrix();
      translate(v.x,v.y,v.z);  
      box(5);
      popMatrix();
    }
  }

  void drawEdgeCenters(){
    Iterator vertexItr = edgeCenters.iterator();
    while(vertexItr.hasNext()){
      HE_Vertex v=(HE_Vertex)vertexItr.next();
      pushMatrix();
      translate(v.x,v.y,v.z);  
      box(5);
      popMatrix();
    }
  }

  void drawEdges(boolean all){
    Iterator eItr = edges.iterator();
    while(eItr.hasNext()){
      HE_Edge e=(HE_Edge)eItr.next();
      if((e.visible)||(all))line(e.halfEdge.vert.x,e.halfEdge.vert.y,e.halfEdge.vert.z,e.halfEdge.next.vert.x,e.halfEdge.next.vert.y,e.halfEdge.next.vert.z);
    }
  }

  void drawEdgesSketch(int r, float d){
    for(int i=0;i<r;i++){
      Iterator eItr = edges.iterator();
      while(eItr.hasNext()){
        HE_Edge e=(HE_Edge)eItr.next();
        if(e.visible)line(e.halfEdge.vert.x+random(-d,d),e.halfEdge.vert.y+random(-d,d),e.halfEdge.vert.z+random(-d,d),e.halfEdge.next.vert.x+random(-d,d),e.halfEdge.next.vert.y+random(-d,d),e.halfEdge.next.vert.z+random(-d,d));
      }
    }
  }

  void drawHalfEdges(){
    Iterator heItr = halfEdges.iterator();
    while(heItr.hasNext()){
      HE_HalfEdge he=(HE_HalfEdge)heItr.next();
      line(he.vert.x,he.vert.y,he.vert.z,he.next.vert.x,he.next.vert.y,he.next.vert.z);
    }
  }

  void drawFaceNormals(float d){
    if(!faceCentersUpdated) setFaceCenters();
    if(!faceNormalsUpdated) setFaceNormals();
    Iterator fnItr = faceNormals.iterator();
    Iterator fcItr = faceCenters.iterator();
    while(fnItr.hasNext()){
      HE_Vertex fc=(HE_Vertex)fcItr.next();
      PVector fn=(PVector)fnItr.next();
      line(fc.x,fc.y,fc.z,(fc.x+d*fn.x),(fc.y+d*fn.y),(fc.z+d*fn.z));
    }
  }

  void drawEdgeNormals(float d){
    if(!edgeCentersUpdated) setEdgeCenters();
    if(!edgeNormalsUpdated) setEdgeNormals();
    Iterator enItr = edgeNormals.iterator();
    Iterator ecItr = edgeCenters.iterator();
    while(enItr.hasNext()){
      HE_Vertex ec=(HE_Vertex)ecItr.next();
      PVector en=(PVector)enItr.next();
      line(ec.x,ec.y,ec.z,(ec.x+d*en.x),(ec.y+d*en.y),(ec.z+d*en.z));
    }
  }

  void drawVertexNormals(float d){
    if(!vertexNormalsUpdated) setVertexNormals();
    Iterator vItr = vertices.iterator();
    Iterator vnItr = vertexNormals.iterator();
    while(vnItr.hasNext()){
      HE_Vertex v=(HE_Vertex)vItr.next();
      PVector vn=(PVector)vnItr.next();
      line(v.x,v.y,v.z,(v.x+d*vn.x),(v.y+d*vn.y),(v.z+d*vn.z));
    }
  }

  void drawCurves(){
    Iterator faceItr = faces.iterator();
    while(faceItr.hasNext()){
      HE_Face face=(HE_Face)faceItr.next();
      HE_HalfEdge halfEdge=face.halfEdge;
      if(face.visible){
        HE_Vertex u=halfEdge.vert;
        HE_Vertex v=halfEdge.next.vert;
        beginShape();
        vertex((0.5*(u.x+v.x)),(0.5*(u.y+v.y)),(0.5*(u.z+v.z)));     
        do{
          u=halfEdge.next.vert;
          v=halfEdge.next.next.vert;      
          bezierVertex(u.x,u.y,u.z,u.x,u.y,u.z,(0.5*(u.x+v.x)),(0.5*(u.y+v.y)),(0.5*(u.z+v.z)));
          halfEdge=halfEdge.next;
        }
        while(halfEdge!=face.halfEdge);
        endShape(CLOSE);
      }
    }
  }

  //MESH CREATORS

  HE_Mesh buildMesh(float[][] simpleVertices, int[][] simpleFaces){
    for(int i=0;i<simpleVertices.length;i++){  
     vertices.add(new HE_Vertex(simpleVertices[i][0],simpleVertices[i][1],simpleVertices[i][2],i)); 
    }
    //Create a loop of halfedges for each face. We assume the vertices are given in a consistent order.
    // To extend this for a non-ordered list, sort the vertices of the input vertices here.
    HE_HalfEdge he;
    for(int i=0;i<simpleFaces.length;i++){
      ArrayList faceEdges=new ArrayList();
      HE_Face hef=new HE_Face();
      hef.visible=true;
      faces.add(hef);
      //TODO: sort face vertices here
      for(int j=0;j<simpleFaces[i].length;j++){
        he=new HE_HalfEdge();
        faceEdges.add(he);
        he.face=hef;
        if (hef.halfEdge==null) hef.halfEdge=he;
        he.vert=(HE_Vertex)vertices.get(simpleFaces[i][j]);
        if(he.vert.halfEdge==null) he.vert.halfEdge=he;  
      }
      cycleHalfEdges(faceEdges,false);    
      halfEdges.addAll(faceEdges);
    }
    pairHalfEdges();
    for(int i=0;i<edges.size();i++){
      ((HE_Edge)edges.get(i)).visible=true;
    }
    return this;
  }

  HE_Mesh get(){

    if(!indicesUpdated) reindex();
    HE_Mesh result=new HE_Mesh();
    for(int i=0;i<vertices.size();i++){
      result.vertices.add(((HE_Vertex)vertices.get(i)).get());
    }
    for(int i=0;i<faces.size();i++){
      result.faces.add(new HE_Face());
    }
    for(int i=0;i<halfEdges.size();i++){
      result.halfEdges.add(new HE_HalfEdge());
    }
    for(int i=0;i<edges.size();i++){
      result.edges.add(new HE_Edge());
    }
   
    for(int i=0;i<vertices.size();i++){
      HE_Vertex sv=(HE_Vertex)vertices.get(i);
      HE_Vertex tv=(HE_Vertex)result.vertices.get(i);
      tv.halfEdge=(HE_HalfEdge)result.halfEdges.get(sv.halfEdge.id);
    }
    for(int i=0;i<faces.size();i++){
      HE_Face sf=(HE_Face)faces.get(i);
      HE_Face tf=(HE_Face)result.faces.get(i);
      tf.id=sf.id;
      tf.halfEdge=(HE_HalfEdge)result.halfEdges.get(sf.halfEdge.id);
      tf.visible=sf.visible;
      tf.fixed=sf.fixed;
    }
    for(int i=0;i<edges.size();i++){
      HE_Edge se=(HE_Edge)edges.get(i);
      HE_Edge te=(HE_Edge)result.edges.get(i);
      te.halfEdge=(HE_HalfEdge)result.halfEdges.get(se.halfEdge.id);
      te.id=se.id;
      te.visible=se.visible;
    }
    for(int i=0;i<halfEdges.size();i++){
      HE_HalfEdge she=(HE_HalfEdge)halfEdges.get(i);
      HE_HalfEdge the=(HE_HalfEdge)result.halfEdges.get(i);
      the.pair=(HE_HalfEdge)result.halfEdges.get(she.pair.id);
      the.next=(HE_HalfEdge)result.halfEdges.get(she.next.id);
      the.vert=(HE_Vertex)result.vertices.get(she.vert.id);
      the.face=(HE_Face)result.faces.get(she.face.id);
      the.edge=(HE_Edge)result.edges.get(she.edge.id);
      the.id=she.id;
    }
   result.bodyCenterUpdated=bodyCenterUpdated;
    if(bodyCenterUpdated) result.bodyCenter.set(bodyCenter);
    result.faceCentersUpdated=faceCentersUpdated;
    if(faceCentersUpdated){
      for(int i=0;i<faceCenters.size();i++){
        result.faceCenters.add(((PVector)faceCenters.get(i)).get());
      }      
    }
    result.faceNormalsUpdated=faceNormalsUpdated;
    if(faceNormalsUpdated){
      for(int i=0;i<faceNormals.size();i++){
        result.faceNormals.add(((PVector)faceNormals.get(i)).get());
      }      
    }
    result.vertexNormalsUpdated=vertexNormalsUpdated;
    if(vertexNormalsUpdated){
      for(int i=0;i<vertexNormals.size();i++){
        result.vertexNormals.add(((PVector)vertexNormals.get(i)).get());
      }      
    }
    result.edgeCentersUpdated=edgeCentersUpdated;
    if(edgeCentersUpdated){
      for(int i=0;i<edgeCenters.size();i++){
        result.edgeCenters.add(((PVector)edgeCenters.get(i)).get());
      }      
    }
    result.edgeNormalsUpdated=edgeNormalsUpdated;
    if(edgeNormalsUpdated){
      for(int i=0;i<edgeNormals.size();i++){
        result.edgeNormals.add(((PVector)edgeNormals.get(i)).get());
      }      
    }
    result.indicesUpdated=indicesUpdated;
    if(indicesUpdated) result.reindex();
    return result;
  }
  
  void set(HE_Mesh target){
    HE_Mesh result=target.get();
   
   vertices=result.vertices; 
    faces=result.faces; 
    halfEdges=result.halfEdges;
    edges=result.edges;
    faceCenters=result.faceCenters;
    faceNormals=result.faceNormals;
    vertexNormals=result.vertexNormals;
    edgeCenters=result.edgeCenters;
    edgeNormals=result.edgeNormals;
    bodyCenter=result.bodyCenter ;
    bodyCenterUpdated=result.bodyCenterUpdated;
    faceCentersUpdated=result.faceCentersUpdated;
    faceNormalsUpdated=result.faceNormalsUpdated;
    vertexNormalsUpdated=result.vertexNormalsUpdated;
    edgeCentersUpdated=result.edgeCentersUpdated;
    edgeNormalsUpdated=result.edgeNormalsUpdated;
    indicesUpdated=result.indicesUpdated;
    
  }
    

  HE_Mesh getDual(boolean rescale){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_Mesh result=new HE_Mesh();
    result.vertices.addAll(faceCenters);
    result.reindex();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      HE_HalfEdge he=v.halfEdge;
      HE_Face f= he.face;
      ArrayList faceHalfEdges=new ArrayList();
      HE_Face nf=new HE_Face();
      result.faces.add(nf);
      do{
        HE_HalfEdge hen = new HE_HalfEdge();
        faceHalfEdges.add(hen);
        hen.face=nf;
        nf.visible=true;
        hen.vert=(HE_Vertex)result.vertices.get(f.id);
        if(hen.vert.halfEdge==null) hen.vert.halfEdge=hen;
        if (nf.halfEdge==null) nf.halfEdge=hen;
        he=he.pair.next;
        f=he.face;
      }
      while(he!=v.halfEdge);
      cycleHalfEdges(faceHalfEdges,true);    
      result.halfEdges.addAll(faceHalfEdges);
    }
    result.pairHalfEdges();
    for(int i=0;i<result.edges.size();i++){
      ((HE_Edge)result.edges.get(i)).visible=true;
    }
    if(rescale){
      float dmax=0f;
      for(int i=0;i<vertices.size();i++){
        HE_Vertex v=(HE_Vertex)vertices.get(i);
        dmax=max(0,dist(v.x,v.y,v.z,0,0,0));

      }
      float rdmax=0f;
      for(int i=0;i<result.vertices.size();i++){
        HE_Vertex v=(HE_Vertex)result.vertices.get(i);
        rdmax=max(0,dist(v.x,v.y,v.z,0,0,0));
      }
      result.rescale(dmax/rdmax);
    }
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    result.indicesUpdated=false;
    return result;
  }

HE_Mesh getSurfaceOfRevolution(PVector[] curve, PVector p1, PVector p2,int n){
   int c=curve.length;
   float a=TWO_PI/n;
   HE_Mesh result=new HE_Mesh();
   for(int i=0;i<n;i++){
     for(int j=0;j<c;j++){
       HE_Vertex v=new HE_Vertex(curve[j],0);
       if(i>0) v.set(rotatePointAboutAxis(v,i*a,p1,p2));
       result.vertices.add(v);      
     }
   }
   for(int i=0;i<n;i++){
    int ip=(i+1)%n;
     for(int j=0;j<c;j++){
       int jp=(j+1)%c;
       HE_Face f=new HE_Face();
       f.visible=true;
       HE_HalfEdge he0=new HE_HalfEdge();
       HE_HalfEdge he1=new HE_HalfEdge();
       HE_HalfEdge he2=new HE_HalfEdge();
       HE_HalfEdge he3=new HE_HalfEdge();
       he0.face=f;he1.face=f;he2.face=f;he3.face=f;
       he0.next=he1;he1.next=he2;he2.next=he3;he3.next=he0;
       he0.vert=(HE_Vertex)result.vertices.get(i*c+j);
       he1.vert=(HE_Vertex)result.vertices.get(i*c+jp);
       he2.vert=(HE_Vertex)result.vertices.get(ip*c+jp);
       he3.vert=(HE_Vertex)result.vertices.get(ip*c+j);
       he0.vert.halfEdge=he0;
       he1.vert.halfEdge=he1;
       he2.vert.halfEdge=he2;
       he3.vert.halfEdge=he3;
       f.halfEdge=he0;
       result.faces.add(f);
       result.halfEdges.add(he0);
       result.halfEdges.add(he1);
       result.halfEdges.add(he2);
       result.halfEdges.add(he3);    
     }
    }

    for(int i=0;i<n;i++){
    int ip=(i+1)%n;
     for(int j=0;j<c;j++){
       int jp=(j+1)%c;
       HE_HalfEdge he1=(HE_HalfEdge)result.halfEdges.get(4*(j+i*c)+1);
       HE_HalfEdge he2=(HE_HalfEdge)result.halfEdges.get(4*(j+i*c)+2);
       
       HE_HalfEdge he0n=(HE_HalfEdge)result.halfEdges.get(4*(j+ip*c));
       HE_HalfEdge he3n=(HE_HalfEdge)result.halfEdges.get(4*(jp+i*c)+3);
       
       he2.pair=he0n;
       he0n.pair=he2;
       he1.pair=he3n;
       he3n.pair=he1;
       
       HE_Edge e0=new HE_Edge();
       HE_Edge e1=new HE_Edge();
       
       e0.halfEdge=he1;
       he1.edge=e0;
       he3n.edge=e0;
              e1.halfEdge=he2;
       he2.edge=e1;
       he0n.edge=e1;
       result.edges.add(e0);
       result.edges.add(e1);
       
     }
    }
    
      return result;
    
    
  }


  HE_Mesh getBox(PVector c, float W, float H, float D){
    float hW=0.5f*W;
float hH=0.5f*H;
float hD=0.5f*D;
    float[][] vertices={
      {
        c.x+hW,c.y+hH,c.z+hD                                                                                                                }
      ,
      {
        c.x-hW,c.y+hH,c.z+hD                                                                                                                }
      ,
      {
        c.x-hW,c.y+hH,c.z-hD                                                                                                                }
      ,
      {
        c.x+hW,c.y+hH,c.z-hD                                                                                                                }
      ,
      {
        c.x+hW,c.y-hH,c.z+hD                                                                                                                }
      ,
      {
        c.x-hW,c.y-hH,c.z+hD                                                                                                                }
      ,
      {
        c.x-hW,c.y-hH,c.z-hD                                                                                                                }
      ,
      {
        c.x+hW,c.y-hH,c.z-hD                                                                                                                }
    };
    int[][] faces={
      {
        3,2,1,0                                                                                                                }
      ,
      {
        5,4,0,1                                                                                                                }
      ,
      {
        4,7,3,0                                                                                                                }
      ,
      {
        6,5,1,2                                                                                                                }
      ,
      {
        6,7,4,5                                                                                                                }
      ,
      {
        7,6,2,3                                                                                                                }
    };


    return buildMesh(vertices, faces);


  }

  HE_Mesh getCube(PVector c, float R){

    float S=R/sqrt(3f);
    return getBox(c,S,S,S);


  }


  HE_Mesh getOctahedron(PVector c, float R){

    return getCube(c,R*sqrt(3f)).getDual(false);


  }


  HE_Mesh getIcosahedron(PVector c,float R){
    float phi=0.5f*(sqrt(5f)+1f);
    float S=R/sqrt(1+phi*phi);
    float pS=S*phi;
    float[][] vertices={
      {
        c.x+S,c.y+pS,c.z                                                                                                            }
      ,
      {
        c.x+S,c.y-pS,c.z                                                                                                            }
      ,
      {
        c.x-S,c.y-pS,c.z                                                                                                            }
      ,
      {
        c.x-S,c.y+pS,c.z                                                                                                            }
      ,
      {
        c.x+pS,c.y,c.z+S                                                                                                            }
      ,
      {
        c.x-pS,c.y,c.z+S                                                                                                            }
      ,
      {
        c.x-pS,c.y,c.z-S                                                                                                            }
      ,
      {
        c.x+pS,c.y,c.z-S                                                                                                            }
      ,
      {
        c.x,c.y+S,c.z+pS                                                                                                            }
      ,
      {
        c.x,c.y+S,c.z-pS                                                                                                            }
      ,
      {
        c.x,c.y-S,c.z-pS                                                                                                            }
      ,
      {
        c.x,c.y-S,c.z+pS                                                                                                            }

    };

    int[][] faces={
      {
        8,4,0                                                                                                            }
      ,
      {
        8,0,3                                                                                                            }
      ,
      {
        8,3,5                                                                                                            }
      ,
      {
        8,5,11                                                                                                            }
      ,
      {
        8,11,4                                                                                                            }
      ,
      {
        4,7,0                                                                                                            }
      ,
      {
        0,9,3                                                                                                            }
      ,
      {
        3,6,5                                                                                                            }
      ,
      {
        5,2,11                                                                                                            }
      ,
      {
        11,1,4                                                                                                            }
      ,
      {
        4,1,7                                                                                                            }
      ,
      {
        0,7,9                                                                                                            }
      ,
      {
        3,9,6                                                                                                            }
      ,
      {
        5,6,2                                                                                                            }
      ,
      {
        11,2,1                                                                                                            }
      ,
      {
        10,7,1                                                                                                            }
      ,
      {
        10,9,7                                                                                                            }
      ,
      {
        10,6,9                                                                                                            }
      ,
      {
        10,2,6                                                                                                            }
      ,
      {
        10,1,2                                                                                                            }



    };

    return buildMesh(vertices, faces);


  }

  HE_Mesh getDodecahedron(PVector c,float R){
    float phi=0.5f*(sqrt(5f)+1f);
    float S=sqrt(3)/(1+phi)*R*sqrt(1+phi*phi);
    return getIcosahedron(c,S).getDual(false);


  }

  HE_Mesh getTetrahedron(PVector c,float R){

    float S=R/sqrt(3f);

    float[][] vertices={  
      {
        S,S,S                                                                                                            }
      ,
      {
        S,-S,-S                                                                                                            }
      ,
      {
        -S,S,-S                                                                                                            }
      ,
      {
        -S,-S,S                                                                                                            }

    };

    int[][] faces={
      {
        0,1,2                                                                                                            }
      ,
      {
        0,2,3                                                                                                            }
      ,
      {
        0,3,1                                                                                                            }
      ,
      {
        1,3,2                                                                                                            }

    };

    return buildMesh(vertices, faces);


  }

HE_Mesh getGeodesic(PVector c,float R, int level){
  HE_Mesh result= getIcosahedron(c,R);
  PlanarMidEdgeSubdivision pmes=new PlanarMidEdgeSubdivision();
  result.subdivideSelf(pmes,level).spherify(R,1f);
  return result;


  }


  HE_Mesh getIsoSurface(float[] values, int w, int h, int d, PVector lowerCorner, PVector upperCorner, float threshold, float background){
    HE_Mesh result=new HE_Mesh();
    float[] paddedValues=new float[(w+2)*(h+2)*(d+2)];
    for(int i=0; i<(w+2)*(h+2)*(d+2);i++){
      paddedValues[i]=background;
    }
    for(int i=0; i<w;i++){
      for(int j=0; j<h;j++){
        for(int k=0; k<d;k++){
          paddedValues[index3D(i+1,j+1,k+1,w+2,h+2)]=values[index3D(i,j,k,w,h)];
        }
      }
    }

    HashMap isoVertices=new HashMap();

    // create all vertices, store in hashMap for later reference
    for(int i=1; i<w+1;i++){
      for(int j=1; j<h+1;j++){
        for(int k=1; k<d+1;k++){                    
          float v0=paddedValues[index3D(i,j,k,w+2,h+2)];
          int i0=index3D(i,j,k,w+2,h+2);
          for(int u=-1;u<2;u++){
            for(int v=-1;v<2;v++){
              for(int x=-1;x<2;x++){
                if(!((u==0)&&(v==0)&&(w==0))){
                  float v1=paddedValues[index3D(i+u,j+v,k+x,w+2,h+2)];
                  if(between(threshold,v0,v1)){
                    float f=(threshold-v0)/(v1-v0); 
                    int i1=index3D(i+u,j+v,k+x,w+2,h+2);
                    String id=str(min(i0,i1))+"_"+str(max(i0,i1));
                    if(isoVertices.get(id)==null){
                    HE_Vertex vn=new HE_Vertex(i+u*f,j+v*f,k+x*f,0);
                                                
                    isoVertices.put(str(min(i0,i1))+"_"+str(max(i0,i1)),vn);
                    }
                  }
                }
              }
            }
          }
        }
           
      } 
    }
    for(int i=0; i<w+1;i++){
      for(int j=0; j<h+1;j++){
        for(int k=0; k<d+1;k++){    
          //indices of the 8 cornerpoints of cube
          int[] indices={
            index3D(i,j,k,w+2,h+2),index3D(i+1,j,k,w+2,h+2),index3D(i+1,j,k+1,w+2,h+2),index3D(i,j,k+1,w+2,h+2),index3D(i,j+1,k,w+2,h+2),
            index3D(i+1,j+1,k,w+2,h+2),
            index3D(i+1,j+1,k+1,w+2,h+2),
            index3D(i,j+1,k+1,w+2,h+2)                                        };
          //6 tetrahedra per cube, 4 corners per tetrahedron             
          int[] c0={
            0,0,0,0,0,5                                                            };
          int[] c1={
            2,2,4,6,6,6                                                            };
          int[] c2={
            3,7,6,1,4,1                                                            };
          int[] c3={
            7,6,7,2,1,4                                                            };

          for (int t=0;t<6;t++){

            int type=0;
            int i0=indices[c0[t]];
            int i1=indices[c1[t]];
            int i2=indices[c2[t]];
            int i3=indices[c3[t]];

            float v0=paddedValues[i0];
            float v1=paddedValues[i1];
            float v2=paddedValues[i2];
            float v3=paddedValues[i3];
            if(v0>threshold) type+=1;
            if(v1>threshold) type+=2;
            if(v2>threshold) type+=4;
            if(v3>threshold) type+=8;
            HE_Face f=new HE_Face();
            
            HE_HalfEdge he0=new HE_HalfEdge();
            HE_HalfEdge he1=new HE_HalfEdge();
            HE_HalfEdge he2=new HE_HalfEdge();
            switch(type) {
            case 0:
            case 15:
              break;
            case 1:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 14:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 2:
            he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 13:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 3:
            he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);             

              break;
            case 12:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);             

              break;
            case 4:
            he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 11:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 5:
             he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);             
              break;
             case 10:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i2))+"_"+str(max(i1,i2)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);             
              break;
            case 6:
he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);  
              break;
            case 9:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i1))+"_"+str(max(i0,i1)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);            
              f=new HE_Face();
              he0=new HE_HalfEdge();
              he1=new HE_HalfEdge();
              he2=new HE_HalfEdge();
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i0,i2))+"_"+str(max(i0,i2)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);             

              break;
            case 7:
he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he0.next=he2;
              he1.next=he0;
              he2.next=he1;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);
              break;
            case 8:
              he0.vert=(HE_Vertex)isoVertices.get(str(min(i0,i3))+"_"+str(max(i0,i3)));
              he1.vert=(HE_Vertex)isoVertices.get(str(min(i1,i3))+"_"+str(max(i1,i3)));
              he2.vert=(HE_Vertex)isoVertices.get(str(min(i2,i3))+"_"+str(max(i2,i3)));
              he0.next=he1;
              he1.next=he2;
              he2.next=he0;
              he0.face=f;
              he1.face=f;
              he2.face=f;
              f.halfEdge=he0;
              he0.vert.halfEdge=he0;
              he1.vert.halfEdge=he1;
              he2.vert.halfEdge=he2;
              result.faces.add(f);
              result.halfEdges.add(he0);
              result.halfEdges.add(he1);
              result.halfEdges.add(he2);

              break;
            default:
              break;
            }
          }
        }
      }
    }
    for(int i=0;i<result.halfEdges.size();i++){
     HE_HalfEdge he=(HE_HalfEdge)result.halfEdges.get(i);
      if(!result.vertices.contains(he.vert)) result.vertices.add(he.vert);            
    }  

    for(int i=0;i<result.vertices.size();i++){
      HE_Vertex v=(HE_Vertex)result.vertices.get(i);
      v.x=(v.x-1)/(w-1)*(upperCorner.x-lowerCorner.x)+lowerCorner.x;
      v.y=(v.y-1)/(h-1)*(upperCorner.y-lowerCorner.y)+lowerCorner.y;
      v.z=(v.z-1)/(d-1)*(upperCorner.z-lowerCorner.z)+lowerCorner.z;
    }
    
    result.pairHalfEdges();
    
    for(int i=0;i<result.faces.size();i++){
     ((HE_Face)result.faces.get(i)).visible=true;
    }
        for(int i=0;i<result.edges.size();i++){
     ((HE_Edge)result.edges.get(i)).visible=true;
    }
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    result.indicesUpdated=false;
    return result;
  }


HE_Mesh getVoronoiCell(PVector[] points, int i, int n,HE_Mesh container){
  HE_Mesh result=container.get();
   ArrayList cutPlanes=new ArrayList();
   for(int j=0;j<n;j++){
    if(i!=j){
   
        PVector N=PVector.sub(points[j],points[i]); 
        N.normalize();
        PVector O=PVector.add(points[j],points[i]); // plane origin=point halfway between point i and point j
        O.mult(0.5f);
        Plane P=new Plane(O,N);
        cutPlanes.add(P);   
    }        
        
     }             
  return result.cutMesh(cutPlanes,points[i],true,true,true);  
  
}

HE_Mesh getVoronoiCell(PVector[] points, int i, int n){
 return  getVoronoiCell(points, i,n,this);
  
}


HE_Mesh[] getAllVoronoiCells(PVector[] points, int n){
  HE_Mesh[] result=new HE_Mesh[n];
  for(int i=0;i<n;i++){
   result[i]=getVoronoiCell(points,i,n,this) ;
    
  }
  return result;
}

HE_Mesh[] getAllVoronoiCells(PVector[] points, int n, HE_Mesh container){
  HE_Mesh result[]=new HE_Mesh[n];
  for(int i=0;i<n;i++){
   result[i]=getVoronoiCell(points,i,n,container) ;
    
  }
  return result;
}

  //MODIFIERS
  HE_Mesh rescale(float f){
    if(!bodyCenterUpdated) setBodyCenter();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      v.x=bodyCenter.x+f*(v.x-bodyCenter.x);
      v.y=bodyCenter.y+f*(v.y-bodyCenter.y);
      v.z=bodyCenter.z+f*(v.z-bodyCenter.z);
    }
    faceCentersUpdated=false;
    edgeCentersUpdated=false;
    return this;
  }

  HE_Mesh rescale(float fx,float fy, float fz){
    if(!bodyCenterUpdated) setBodyCenter();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      v.x=bodyCenter.x+fx*(v.x-bodyCenter.x);
      v.y=bodyCenter.y+fy*(v.y-bodyCenter.y);
      v.z=bodyCenter.z+fz*(v.z-bodyCenter.z);
    }
    faceCentersUpdated=false;
    edgeCentersUpdated=false;
    return this;
  }

  HE_Mesh extremize(float f){
    if(!bodyCenterUpdated) setBodyCenter();
    float dmean=0f;
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      dmean+=dist(v.x,v.y,v.z,bodyCenter.x,bodyCenter.y,bodyCenter.z);
    }
    dmean/=vertices.size();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      float d=dist(v.x,v.y,v.z,bodyCenter.x,bodyCenter.y,bodyCenter.z);
      float ff=1+f*(d-dmean)/dmean;
      v.x=bodyCenter.x+ff*(v.x-bodyCenter.x);
      v.y=bodyCenter.y+ff*(v.y-bodyCenter.y);
      v.z=bodyCenter.z+ff*(v.z-bodyCenter.z);
    }
    clearSecondaries();
    return this;
  }

  HE_Mesh spherify(float r, float f){
    if(!bodyCenterUpdated) setBodyCenter();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      float d=dist(v.x,v.y,v.z,bodyCenter.x,bodyCenter.y,bodyCenter.z);
      float ff=((1-f)*d+f*r)/d;
      v.x=bodyCenter.x+ff*(v.x-bodyCenter.x);
      v.y=bodyCenter.y+ff*(v.y-bodyCenter.y);
      v.z=bodyCenter.z+ff*(v.z-bodyCenter.z);
    }
    clearSecondaries();
    return this;
  }


  //GEOMETRIC OPERATIONS 
  HE_Mesh collapseVertex(int id){
    HE_Vertex v=(HE_Vertex)vertices.get(id);
    HE_HalfEdge he=v.halfEdge;
    float x=0f;
    float y=0f;
    float z=0f;
    int c=0;
    do{
      c++;
      x+=he.next.vert.x;
      y+=he.next.vert.y;
      z+=he.next.vert.z;
      he=he.pair.next;
    }      
    while(he!=v.halfEdge);
    v.set(x/c,y/c,z/c);
    return this;

  }

  HE_Mesh triangulate(){
    if(!indicesUpdated) reindex();
    ArrayList quadFaces=new ArrayList();
    ArrayList polyFaces=new ArrayList();
    for(int i=0;i<faces.size();i++){
      HE_Face f=(HE_Face)faces.get(i);
      HE_HalfEdge he=f.halfEdge;
      int cv=0;
      do{
        cv++;
        he=he.next;
      }
      while(he!=f.halfEdge);
      if(cv==4) quadFaces.add(f);
      if(cv>4) polyFaces.add(f);
    }
    for(int i=0;i<quadFaces.size();i++){
      HE_Face f=(HE_Face)quadFaces.get(i);
      HE_HalfEdge he1=f.halfEdge;
      HE_HalfEdge he2=he1.next;
      HE_HalfEdge he3=he2.next;
      HE_HalfEdge he4=he3.next;
      HE_HalfEdge hen1=new HE_HalfEdge(); 
      HE_HalfEdge hen2=new HE_HalfEdge();      
      HE_Edge en=new HE_Edge();      
      HE_Face fn=new HE_Face();
      halfEdges.add(hen1);
      halfEdges.add(hen2); 
      edges.add(en); 
      faces.add(fn);
      en.halfEdge=hen1;
      en.visible=true;
      fn.halfEdge=hen2;
      fn.visible=f.visible;
      fn.fixed=f.fixed;
      hen1.vert=he3.vert;
      hen2.vert=he1.vert;
      hen1.pair=hen2;
      hen2.pair=hen1;
      hen1.edge=en;
      hen2.edge=en;
      hen1.face=f;
      hen2.face=fn;
      he3.face=fn;
      he4.face=fn;
      hen1.next=he1;
      hen2.next=he3;
      he2.next=hen1;
      he4.next=hen2;
    }
    clearSecondaries();
    return this;
  }

  HE_Mesh splitEdge(int id, float f){
    
    HE_Edge origE=(HE_Edge)edges.get(id);
    if((f>0f)&&(f<1f)){
      HE_HalfEdge origHE1= origE.halfEdge;
      HE_HalfEdge origHE2= origHE1.pair;
      HE_Vertex v1=origHE1.vert;
      HE_Vertex v2=origHE2.vert;
      HE_Vertex v;
      v=new HE_Vertex(0,0,0,0); 
      v.x=v1.x+f*(v2.x-v1.x);
      v.y=v1.y+f*(v2.y-v1.y);
      v.z=v1.z+f*(v2.z-v1.z);
      vertices.add(v);      
      HE_HalfEdge newHE1=new HE_HalfEdge();
      halfEdges.add(newHE1);
      HE_HalfEdge newHE2=new HE_HalfEdge();
      halfEdges.add(newHE2);
      newHE1.vert=v;    
      newHE1.next=origHE1.next;
      newHE1.face=origHE1.face;
      newHE1.pair=origHE2;
      origHE1.next=newHE1;
      origHE1.pair=newHE2;        
      newHE2.vert=v;
      newHE2.next=origHE2.next;
      newHE2.face=origHE2.face;
      newHE2.pair=origHE1;   
      origHE2.next=newHE2;
      origHE2.pair=newHE1;
      v.halfEdge=newHE1;
      HE_Edge newE= new HE_Edge();
      edges.add(newE);
      newE.halfEdge=newHE1;
      newHE1.edge=newE;
      origHE2.edge=newE;
      newHE2.edge=origE;
      origHE1.edge=origE;
      newE.visible=origE.visible;
    }
clearSecondaries();
    return this;
  }

  HE_Mesh divideEdge(int id, int n){
    if(n>1){
    
    HE_Edge origE=(HE_Edge)edges.get(id);
    float f=1f/n;
   
      HE_HalfEdge origHE1= origE.halfEdge;
      HE_HalfEdge origHE2= origHE1.pair;
      HE_Vertex v1=origHE1.vert;
      HE_Vertex v2=origHE2.vert;
      HE_Vertex v;
      v=new HE_Vertex(0,0,0,0); 
      v.x=v1.x+f*(v2.x-v1.x);
      v.y=v1.y+f*(v2.y-v1.y);
      v.z=v1.z+f*(v2.z-v1.z);
      vertices.add(v);      
      HE_HalfEdge newHE1=new HE_HalfEdge();
      halfEdges.add(newHE1);
      HE_HalfEdge newHE2=new HE_HalfEdge();
      halfEdges.add(newHE2);
      newHE1.vert=v;    
      newHE1.next=origHE1.next;
      newHE1.face=origHE1.face;
      newHE1.pair=origHE2;
      origHE1.next=newHE1;
      origHE1.pair=newHE2;        
      newHE2.vert=v;
      newHE2.next=origHE2.next;
      newHE2.face=origHE2.face;
      newHE2.pair=origHE1;   
      origHE2.next=newHE2;
      origHE2.pair=newHE1;
      v.halfEdge=newHE1;
      HE_Edge newE= new HE_Edge();
      edges.add(newE);
      newE.halfEdge=newHE1;
      newHE1.edge=newE;
      origHE2.edge=newE;
      newHE2.edge=origE;
      origHE1.edge=origE;
      newE.visible=origE.visible;
    divideEdge(edges.size()-1, n-1);
    clearSecondaries();
    }
    return this;
  }




  HE_Mesh cutSurface(Plane P){
    if(!indicesUpdated) reindex();
    ArrayList verticesOnPlane=new ArrayList();
    ArrayList newVertices=new ArrayList();
    int n=edges.size();
    int numOfIntersections=0;
    for(int i=0;i<n;i++){
      HE_Edge e=(HE_Edge)edges.get(i);
      float u=e.intersectionWithPlane(P,.0001f);
      if(u==0f){
        if(!verticesOnPlane.contains(e.halfEdge.vert))verticesOnPlane.add(e.halfEdge.vert); 
        numOfIntersections++;
      }
      else if(u==1f){
        if(!verticesOnPlane.contains(e.halfEdge.pair.vert))verticesOnPlane.add(e.halfEdge.pair.vert); 
        numOfIntersections++;
      }
      else if(u==-2f){
        if(!verticesOnPlane.contains(e.halfEdge.vert))verticesOnPlane.add(e.halfEdge.vert); 
        if(!verticesOnPlane.contains(e.halfEdge.pair.vert))verticesOnPlane.add(e.halfEdge.pair.vert); 
        numOfIntersections++;
      }
      else if(u>0f)
      {
        splitEdge(i,u);
        newVertices.add((HE_Vertex)vertices.get(vertices.size()-1));
        numOfIntersections++;
      }
    }

    for(int i=0;i<verticesOnPlane.size();i++){
      ((HE_Vertex)verticesOnPlane.get(i)).id=-2;
    }
    for(int i=0;i<newVertices.size();i++){
      ((HE_Vertex)newVertices.get(i)).id=-1;
    }
    ArrayList newFaces=new ArrayList();
    for(int i=0;i<faces.size();i++){
      HE_Vertex cv1=null;
      HE_Vertex cv2=null;
      HE_HalfEdge she1=null;
      HE_HalfEdge she2=null;
      HE_HalfEdge ehe1=null;
      HE_HalfEdge ehe2=null;
      HE_Face f=(HE_Face)faces.get(i);
      HE_HalfEdge he=f.halfEdge;
      do{
        if(he.vert.id<0){
          if(cv1==null){
            cv1=he.vert;
            she1=he;
          }
          else{
            cv2=he.vert;
            she2=he;
          } 
        }
        he=he.next;
      }
      while(he!= f.halfEdge);
      if((cv1!=null)&&(cv2!=null)&&!((cv1.id==-2)&&(cv2.id==-2))){
        he=f.halfEdge;
        do{
          if(he.next==she2) ehe1=he;
          if(he.next==she1) ehe2=he;
          he=he.next;
        }
        while(he!= f.halfEdge);                                                        
        f.halfEdge=she1;
        HE_HalfEdge nhe1=new HE_HalfEdge();
        HE_HalfEdge nhe2=new HE_HalfEdge();
        nhe1.vert=cv2;
        nhe2.vert=cv1;
        nhe1.pair=nhe2;
        nhe2.pair=nhe1;
        nhe1.next=she1;
        ehe1.next=nhe1;
        nhe2.next=she2;
        ehe2.next=nhe2;
        HE_Edge ne= new HE_Edge();
        ne.visible=true;
        ne.halfEdge=nhe1;
        nhe1.edge=ne;
        nhe2.edge=ne;
        HE_Face nf=new HE_Face();
        nf.halfEdge=she2;
        nf.visible=f.visible;
        nf.fixed=f.fixed;
        he=she1;
        do{       
          he.face=f;
          he=he.next;
        }
        while(he!=she1);
        he=she2;
        do{       
          he.face=nf;
          he=he.next;
        }
        while(he!=she2);
       halfEdges.add(nhe1);
        halfEdges.add(nhe2);
        edges.add(ne);
        newFaces.add(nf);
      }
    }
    faces.addAll(newFaces);
    clearSecondaries();
    return this; 
  }

  HE_Mesh cutSurface(ArrayList cutPlanes){
    
    for(int k=0;k<cutPlanes.size();k++){
      Plane P = (Plane)cutPlanes.get(k);          
      cutSurface(P);
    }
    return this;
  }

 HE_Mesh cutMesh(Plane P, boolean safe, boolean closed){ 
 return cutMesh(P, safe, closed, false);
 }

  HE_Mesh cutMesh(Plane P, boolean safe, boolean closed, boolean reverse){    
    float pol=(reverse)?-1:1;
    lastOperationFailed=false;
    if(vertices.size()==0) return this;
    
  
    float checkref=P.side((HE_Vertex)vertices.get(0),0.0001f);
    float check=0f;
    for(int i=1;i<vertices.size();i++){
      check=checkref*P.side((HE_Vertex)vertices.get(i),0.0001f);
      if (check<=0f) break;      
    }
    if (check>0f){
      if(!bodyCenterUpdated)setBodyCenter();
       if(pol*P.side(bodyCenter,0.0001f)>=0f){
         return this;
       }else{
         HE_Mesh nothing=new HE_Mesh();
         set(nothing);
         return this;
       }
    }
       
      HE_Mesh backup=new HE_Mesh();
    if(safe) backup=get();
    backup.lastOperationFailed=true;
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    ArrayList verticesOnPlane=new ArrayList();
    ArrayList newVertices=new ArrayList();
    int n=edges.size();
    int numOfIntersections=0;
    for(int i=0;i<n;i++){
      HE_Edge e=(HE_Edge)edges.get(i);
      float u=e.intersectionWithPlane(P,.0001f);
      if(u==0f){
        if(!verticesOnPlane.contains(e.halfEdge.vert))verticesOnPlane.add(e.halfEdge.vert); 
        numOfIntersections++;
      }
      else if(u==1f){
        if(!verticesOnPlane.contains(e.halfEdge.pair.vert))verticesOnPlane.add(e.halfEdge.pair.vert); 
        numOfIntersections++;
      }
      else if(u==-2f){
        if(!verticesOnPlane.contains(e.halfEdge.vert))verticesOnPlane.add(e.halfEdge.vert); 
        if(!verticesOnPlane.contains(e.halfEdge.pair.vert))verticesOnPlane.add(e.halfEdge.pair.vert); 
        numOfIntersections++;
      }
      else if(u>0f)
      {
        splitEdge(i,u);
        newVertices.add((HE_Vertex)vertices.get(vertices.size()-1));
        numOfIntersections++;
      }
    }


    for(int i=0;i<verticesOnPlane.size();i++){
      ((HE_Vertex)verticesOnPlane.get(i)).id=-2;
    }
    for(int i=0;i<newVertices.size();i++){
      ((HE_Vertex)newVertices.get(i)).id=-1;
    }
    ArrayList newFaces=new ArrayList();
    for(int i=0;i<faces.size();i++){
      HE_Vertex cv1=null;
      HE_Vertex cv2=null;
      HE_HalfEdge she1=null;
      HE_HalfEdge she2=null;
      HE_HalfEdge ehe1=null;
      HE_HalfEdge ehe2=null;
      HE_Face f=(HE_Face)faces.get(i);
      HE_HalfEdge he=f.halfEdge;
      do{
        if(he.vert.id<0){
          if(cv1==null){
            cv1=he.vert;
            she1=he;
          }
          else{
            cv2=he.vert;
            she2=he;
          } 
        }
        he=he.next;
      }
      while(he!= f.halfEdge);
      if((cv1!=null)&&(cv2!=null)&&!((cv1.id==-2)&&(cv2.id==-2))){
        he=f.halfEdge;
        do{
          if(he.next==she2) ehe1=he;
          if(he.next==she1) ehe2=he;
          he=he.next;
        }
        while(he!= f.halfEdge);                                                        
        f.halfEdge=she1;
        HE_HalfEdge nhe1=new HE_HalfEdge();
        HE_HalfEdge nhe2=new HE_HalfEdge();
        nhe1.vert=cv2;
        cv2.halfEdge=nhe1;
        nhe2.vert=cv1;
        cv1.halfEdge=nhe2;
        nhe1.pair=nhe2;
        nhe2.pair=nhe1;
        nhe1.next=she1;
        ehe1.next=nhe1;
        nhe2.next=she2;
        ehe2.next=nhe2;
        HE_Edge ne= new HE_Edge();
        ne.visible=true;
        ne.halfEdge=nhe1;
        nhe1.edge=ne;
        nhe2.edge=ne;
        HE_Face nf=new HE_Face();
        nf.visible=f.visible;
        nf.fixed=f.fixed;
        nf.halfEdge=she2;
        he=she1;
        do{       
          he.face=f;
          he=he.next;
        }
        while(he!=she1);
        he=she2;
        do{       
          he.face=nf;
          he=he.next;
        }
        while(he!=she2);
       halfEdges.add(nhe1);
        halfEdges.add(nhe2);
       edges.add(ne);
        if(pol*P.side(getFaceCenter(nf),0.0001f)>=0f){
          newFaces.add(nf);
        }
        else{
          newFaces.add(f); 
        }
      }
      else{
        if(pol*P.side((PVector)faceCenters.get(i),0.0001f)>=0f)newFaces.add(f);
      }
    }
   faces=newFaces;        
    cleanUnusedElements();                    
    capOpenLoop(closed);
    clearSecondaries();
    if(safe){
    return (validate(false, false))?this:backup;  
    }
    else
    {    
    return this; 
    }
  }



HE_Mesh[] splitMesh(final Plane P, boolean safe, boolean closed){    

 
 HE_Mesh[] result=new HE_Mesh[2];
 result[0]=get();
 
 result[1]=get();
 result[0].cutMesh(P,safe,closed);
 result[1].cutMesh(P,safe,closed, true);
 
  return result;
}


  HE_Mesh cutMesh(ArrayList cutPlanes, PVector c, boolean safe, boolean closed, boolean retry){
   ArrayList cutPlanesToRetry=new ArrayList();
    float[] r=new float[cutPlanes.size()];
    for(int k=0;k<cutPlanes.size();k++){
      Plane P = (Plane)cutPlanes.get(k);
      r[k]=dist(P.o.x,P.o.y,P.o.z,c.x,c.y,c.z);
    }
    for (int k = cutPlanes.size(); --k >= 0; ) {
      for (int m = 0; m < k; m++) {
        Plane pl = (Plane)cutPlanes.get(m);       
        Plane pll = (Plane)cutPlanes.get(m+1);
        if (r[m] > r[m+1]) {
          Collections.swap(cutPlanes,m,m+1);
          float tmp=r[m];
          r[m]=r[m+1];
          r[m+1]=tmp;
        }
      }
    }
    for(int k=0;k<cutPlanes.size();k++){
      Plane P = (Plane)cutPlanes.get(k);          
      cutMesh(P,safe,closed);
      if(lastOperationFailed) {
        cutPlanesToRetry.add(P);
         P.o.add(random(-.1,.1), random(-.1,.1),random(-.1,.1));
   P.update();      
      }
    }
     Collections.shuffle(cutPlanesToRetry);
    if((retry)&&(cutPlanesToRetry.size()>0)){
      
      cutMesh(cutPlanesToRetry,c,safe, closed, false);
      
    }
    
    return this;
  }
 

  HE_Mesh extrudeFace(int id, float d, float t){
    if(!faceNormalsUpdated) setFaceNormals();
    if(!faceCentersUpdated) setFaceCenters();
    HE_Face f=(HE_Face)faces.get(id); 
    PVector n=(PVector)faceNormals.get(id);
    PVector fc=(PVector)faceCenters.get(id);
    n.mult(d);   
    HE_HalfEdge he=f.halfEdge;
    ArrayList faceVertices=new ArrayList();
    ArrayList extFaceVertices=new ArrayList();
    do{
      faceVertices.add(he.vert);
      extFaceVertices.add(he.vert.get());       
      he=he.next;                     
    }
    while(he!=f.halfEdge);  

    for (int i=0;i<extFaceVertices.size();i++){
      HE_Vertex v=(HE_Vertex)extFaceVertices.get(i);      
      PVector diff=PVector.sub(fc,v);
      diff.mult(t);
      v.add(diff);
      v.add(n);
      vertices.add(v);
    }
    he=f.halfEdge;
    int c=0;
    do{ 
      HE_Face fNew=new HE_Face();
      fNew.visible=true;
      HE_HalfEdge heOrig1=he;
      HE_HalfEdge heOrig2=he.pair;
      HE_HalfEdge heNew1=new HE_HalfEdge();
      HE_HalfEdge heNew2=new HE_HalfEdge();
      HE_HalfEdge heNew3=new HE_HalfEdge();
      HE_HalfEdge heNew4=new HE_HalfEdge();
      HE_Edge eNew=new HE_Edge();
      eNew.visible=true;
      HE_Vertex v1=(HE_Vertex)faceVertices.get(c);
      HE_Vertex v2=(HE_Vertex)faceVertices.get((c+1)%faceVertices.size());
      HE_Vertex v4=(HE_Vertex)extFaceVertices.get(c);
      HE_Vertex v3=(HE_Vertex)extFaceVertices.get((c+1)%extFaceVertices.size());

      heNew1.vert=v1;
      v1.halfEdge=heNew1;
      heNew1.face=fNew;

      fNew.halfEdge=heNew1;
      heNew1.edge=heOrig2.edge;
      heNew1.edge.halfEdge=heNew1;
      heNew1.pair=heOrig2;
      heOrig2.pair=heNew1;
      heNew1.next=heNew2;

      heNew2.vert=v2;
      v2.halfEdge=heNew2;
      heNew2.face=fNew;
      heNew2.next=heNew3; 

      heNew3.vert=v3;
      v3.halfEdge=heNew3;
      heNew3.face=fNew;
      heNew3.edge=eNew;
      heOrig1.edge=eNew;
      eNew.halfEdge=heNew3;
      heNew3.pair=heOrig1;
      heOrig1.pair=heNew3;
      heNew3.next=heNew4;

      heNew4.vert=v4;
      v4.halfEdge=heNew4;
      heNew4.face=fNew;
      heNew4.next=heNew1;

      heOrig1.vert=v4;
      faces.add(fNew);
      edges.add(eNew);
      halfEdges.add(heNew1);
      halfEdges.add(heNew2);      
      halfEdges.add(heNew3);
      halfEdges.add(heNew4);

      he=he.next;
      c++;
    }
    while(he!=f.halfEdge);
    pairHalfEdges();  

    clearSecondaries();
    return this;

  }


  //SUBDIVISION



/*static*/ HE_Mesh subdivide(final HE_Mesh mesh, final HE_Subdivision subDivision, int rep){
 return  subDivision.divide(mesh,rep);  
}

/*static*/ HE_Mesh subdivide(final HE_Mesh mesh,final HE_Subdivision subDivision){  
 return  subDivision.divide(mesh,1);  
}

HE_Mesh subdivide(final HE_Subdivision subDivision, int rep){
 return  subDivision.divide(this,rep);  
}

HE_Mesh subdivide(final HE_Subdivision subDivision){  
 return  subDivision.divide(this,1);  
}

HE_Mesh subdivideSelf(final HE_Subdivision subDivision, int rep){
 set(subDivision.divide(this,rep));  
 return this;
}

HE_Mesh subdivideSelf(final HE_Subdivision subDivision){  
 set(subDivision.divide(this,1));  
 return this;
}





  
  //HELPER FUNCTIONS


  int index3D(int i, int j, int k, int w, int h){
    return(i+w*j+w*h*k); 

  }

  void cycleHalfEdges(ArrayList srcHalfEdges, boolean reverse){   
    HE_HalfEdge he;
    int n=srcHalfEdges.size();
    if(n>0){
      if(!reverse){
        for(int j=0;j<n-1;j++){
          he=(HE_HalfEdge)srcHalfEdges.get(j);
          he.next=(HE_HalfEdge)srcHalfEdges.get(j+1);
        }
        he=(HE_HalfEdge)srcHalfEdges.get(n-1);
        he.next=(HE_HalfEdge)srcHalfEdges.get(0);
      }
      else{
        he=(HE_HalfEdge)srcHalfEdges.get(0);
        he.next=(HE_HalfEdge)srcHalfEdges.get(n-1);
        for(int j=1;j<n;j++){
          he=(HE_HalfEdge)srcHalfEdges.get(j);
          he.next=(HE_HalfEdge)srcHalfEdges.get(j-1);
        }
      }
    }
  }

  //go through all the halfedges and find matching pairs
  void pairHalfEdges(){
    ArrayList unpairedHalfEdges=new ArrayList();
    for(int i=0;i<halfEdges.size();i++){
    HE_HalfEdge he=(HE_HalfEdge)halfEdges.get(i);
    if (he.pair==null) unpairedHalfEdges.add(he);
      
    }
    
    for(int i=0;i<unpairedHalfEdges.size();i++){
      HE_HalfEdge he=(HE_HalfEdge)unpairedHalfEdges.get(i);
      if(he.pair==null){
        for(int j=i+1;j<unpairedHalfEdges.size();j++){
            HE_HalfEdge he2=(HE_HalfEdge)unpairedHalfEdges.get(j);
            if((he2.pair==null)&&(he.vert==he2.next.vert)&&(he2.vert==he.next.vert)){
              he.pair=he2;
              he2.pair=he;
              HE_Edge e=new HE_Edge();
              e.halfEdge=he;
              he.edge=e;
              he2.edge=e;
              edges.add(e);
              break;              
            
          }
        }
      }
    }
  }

  void reindex(){
    for(int i=0;i<vertices.size();i++){
      ((HE_Vertex)vertices.get(i)).id=i;
    }
    for(int i=0;i<faces.size();i++){
      ((HE_Face)faces.get(i)).id=i;
    }
    for(int i=0;i<halfEdges.size();i++){
      ((HE_HalfEdge)halfEdges.get(i)).id=i;
    }
    for(int i=0;i<edges.size();i++){
      ((HE_Edge)edges.get(i)).id=i;
    }
    indicesUpdated=true;
  }


  void cleanUnusedElements(){
    ArrayList cleanedVertices=new ArrayList();
    ArrayList cleanedHalfEdges=new ArrayList(); 

    ArrayList cleanedEdges=new ArrayList();    

    for(int i=0;i<faces.size();i++){
      HE_Face f=(HE_Face)faces.get(i);
      HE_HalfEdge he=f.halfEdge;
      do{
        if(!cleanedVertices.contains(he.vert)) cleanedVertices.add(he.vert);
        if(!cleanedHalfEdges.contains(he)) cleanedHalfEdges.add(he);
        if(he.edge.visible) he.id=-3;
        he.edge=null;
        he=he.next;

      }
      while(he!=f.halfEdge);
    }
    for(int i=0;i<cleanedHalfEdges.size();i++){
      HE_HalfEdge he=(HE_HalfEdge)cleanedHalfEdges.get(i);
      if(!cleanedHalfEdges.contains(he.pair)){
        he.pair=null;      

        he.vert.halfEdge=he;
      }
      else{
        if(he.edge==null){
          HE_Edge e=new HE_Edge();
          e.halfEdge=he;
          he.edge=e;
          he.pair.edge=e;
          if(he.id==-3) e.visible=true;
          cleanedEdges.add(e); 
        }
      }
    }
    vertices=cleanedVertices;
    halfEdges=cleanedHalfEdges;
    edges=cleanedEdges;
    bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
  }

  void capOpenLoop(boolean vis){
    ArrayList unpairedEdges=new ArrayList();
    for(int i=0;i<halfEdges.size();i++){
      HE_HalfEdge he=(HE_HalfEdge)halfEdges.get(i);
      if(he.pair==null) unpairedEdges.add(he);
    } 
    while(unpairedEdges.size()>0){
      ArrayList loopedHalfEdges=new ArrayList();
      HE_HalfEdge start = (HE_HalfEdge)unpairedEdges.get(0); 
      loopedHalfEdges.add(start);
      HE_HalfEdge he=start;
      HE_HalfEdge hen=start;
      boolean stuck=false;    
      do{
        for(int i=0;i<unpairedEdges.size();i++){
          hen=(HE_HalfEdge)unpairedEdges.get(i);
          if(hen.vert==he.next.vert){            
            if(!loopedHalfEdges.contains(hen)){
              loopedHalfEdges.add(hen);
            }
            else{
              stuck=true;
            }
            break;
          }
        }
        if(hen.vert!=he.next.vert) stuck=true;
        he=hen;
      }
      while((hen.next.vert!=start.vert)&&(!stuck));
      unpairedEdges.removeAll(loopedHalfEdges);
      HE_Face nf=new HE_Face();
      nf.visible=vis;
     // nf.fixed=true;
      faces.add(nf);
      ArrayList newHalfEdges=new ArrayList();
      for(int i=0;i<loopedHalfEdges.size();i++){
        HE_HalfEdge phe=(HE_HalfEdge)loopedHalfEdges.get(i);
        HE_HalfEdge nhe=new HE_HalfEdge();
        halfEdges.add(nhe);
        newHalfEdges.add(nhe);
        nhe.vert=phe.next.vert;
        nhe.pair=phe;
        phe.pair=nhe;
        nhe.face=nf;
        if(nf.halfEdge==null)nf.halfEdge=nhe;
        HE_Edge ne=new HE_Edge();
        edges.add(ne);
        ne.visible=true;
        ne.halfEdge=nhe;
        nhe.edge=ne;
        phe.edge=ne;

      }
      cycleHalfEdges(newHalfEdges,true);
    }
    bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
  }



  boolean isFixed(HE_Vertex v){
    HE_HalfEdge he=v.halfEdge;
    boolean result=false;
    do{
      result=(result)||(he.face.fixed);
      he=he.pair.next;
    }
    while((he!=v.halfEdge)&&(!result));
    return result;
  }

  boolean isFixed(HE_Edge e){ 
    return ((e.halfEdge.face.fixed)||(e.halfEdge.pair.face.fixed));
  }

  boolean isFixed(HE_Face f){
    HE_HalfEdge he=f.halfEdge;
    boolean result=true;
    do{
      result=result&&(he.face.fixed||he.pair.face.fixed);
      he=he.next;
    }
    while(he!=f.halfEdge);
    return result;
  }




  //SECONDARIES

  void clearSecondaries(){
    
   bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
  }

  void setBodyCenter(){
    bodyCenter=new PVector(0,0,0);
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);   
      bodyCenter.add(v);
    } 
    bodyCenter.mult(1f/vertices.size());
    bodyCenterUpdated=true;
  }

  void setFaceCenters(){
    faceCenters.clear();
    for(int i=0;i<faces.size();i++){
      HE_Face face=(HE_Face)faces.get(i);
      faceCenters.add(getFaceCenter(face));     
    }
    faceCentersUpdated=true;
  }

  PVector getFaceCenter(HE_Face f){
    HE_HalfEdge he= f.halfEdge;      
    float x=0;
    float y=0;
    float z=0;
    int c=0;
    do{         
      x+=he.vert.x; 
      y+=he.vert.y; 
      z+=he.vert.z; 
      c++;
      he = he.next;
    } 
    while(he!=f.halfEdge);
    return(new HE_Vertex(x/c,y/c,z/c,-1));
  }

  void setFaceNormals(){
    faceNormals.clear();
    for(int i=0;i<faces.size();i++){
      HE_Face face=(HE_Face)faces.get(i);
      HE_HalfEdge halfEdge= face.halfEdge;      
      float x=0;
      float y=0;
      float z=0;
      int c=0;
      do{
        Plane P=new Plane(halfEdge.vert,halfEdge.next.vert,halfEdge.next.next.vert);         
        x+=P.A; 
        y+=P.B; 
        z+=P.C; 
        c++;
        halfEdge = halfEdge.next;
      } 
      while(halfEdge!=face.halfEdge);
      PVector n=new PVector(x/c,y/c,z/c);
      n.normalize();
      faceNormals.add(n);     
    }
    faceNormalsUpdated=true;
  }

  PVector getAverageAdjacentFaceCenter(HE_Vertex v){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_HalfEdge halfEdge= v.halfEdge;      
    float x=0;
    float y=0;
    float z=0;
    int c=0;
    do{
      PVector fn=(PVector)faceCenters.get(halfEdge.face.id);         
      x+=fn.x; 
      y+=fn.y; 
      z+=fn.z; 
      c++;
      halfEdge = halfEdge.pair.next;
    } 
    while(halfEdge!=v.halfEdge);
    PVector n=new PVector(x/c,y/c,z/c);
    return(n);     
  }

  PVector getAverageAdjacentFaceCenter(HE_Edge e){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_HalfEdge halfEdge= e.halfEdge;      
    PVector f1=(PVector)faceCenters.get(halfEdge.face.id); 
    PVector f2=(PVector)faceCenters.get(halfEdge.pair.face.id); 
    PVector n=new PVector(0.5f*(f1.x+f2.x),0.5f*(f1.y+f2.y),0.5f*(f1.z+f2.z));
    return(n);     
  }

  void setVertexNormals(){
    if(!faceNormalsUpdated) setFaceNormals();
    if(!indicesUpdated) reindex();
    vertexNormals.clear();
    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);
      HE_HalfEdge halfEdge= v.halfEdge;      
      float x=0;
      float y=0;
      float z=0;
      int c=0;
      do{
        PVector fn=(PVector)faceNormals.get(halfEdge.face.id);         
        x+=fn.x; 
        y+=fn.y; 
        z+=fn.z; 
        c++;
        halfEdge = halfEdge.pair.next;
      } 
      while(halfEdge!=v.halfEdge);
      PVector n=new PVector(x/c,y/c,z/c);
      n.normalize();
      vertexNormals.add(n);     
    }
    vertexNormalsUpdated=true;
  }

  void setEdgeCenters(){
    edgeCenters.clear();
    for(int i=0;i<edges.size();i++){
      HE_Edge e=(HE_Edge)edges.get(i);
      HE_HalfEdge he1= e.halfEdge;      
      HE_HalfEdge he2= e.halfEdge.pair; 
      float x=0.5*(he1.vert.x+he2.vert.x);
      float y=0.5*(he1.vert.y+he2.vert.y);
      float z=0.5*(he1.vert.z+he2.vert.z);
      edgeCenters.add(new HE_Vertex(x,y,z,-1));     
    }
    edgeCentersUpdated=true;
  }

  void setEdgeNormals(){
    if(!faceNormalsUpdated) setFaceNormals();
    if(!indicesUpdated) reindex();
    edgeNormals.clear();
    for(int i=0;i<edges.size();i++){
      HE_Edge e=(HE_Edge)edges.get(i);
      HE_HalfEdge he= e.halfEdge;      
      HE_Face f1=he.face;
      HE_Face f2=he.pair.face;
      PVector n1=(PVector)faceNormals.get(f1.id);
      PVector n2=(PVector)faceNormals.get(f2.id);
      float x=0.5*(n1.x+n2.x);
      float y=0.5*(n1.y+n2.y);
      float z=0.5*(n1.z+n2.z);
      PVector n=new PVector(x,y,z);
      n.normalize();
      edgeNormals.add(n);     
    }
    edgeNormalsUpdated=true;
  }

  void updateMesh(){   
    reindex();
    setBodyCenter();
    setFaceCenters();
    setFaceNormals();
    setVertexNormals();
    setEdgeCenters();
    setEdgeNormals();
  }

  //VALIDATION 

void analyse(){
 int minVertices=0;
 int maxVertices=0;

for(int i=0;i<faces.size();i++) {
 HE_Face f=(HE_Face)faces.get(i);
HE_HalfEdge he=f.halfEdge;
int c=0;
do{
  c++;
  he=he.next;
  
}while(he!=f.halfEdge);
  minVertices=(i==0)?c:min(c,minVertices);
  maxVertices=(i==0)?c:max(c,minVertices);
}

println(minVertices+" "+maxVertices);
  
}



  boolean validate(boolean verbose, boolean force){
    boolean result =true;
    if(verbose) println("Checking face ("+ faces.size()+") properties");
    for(int i=0;i<faces.size();i++){
      HE_Face face=(HE_Face)faces.get(i);
      if(face.halfEdge==null) {
        if(verbose) println("Null reference in face "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }

      }
      else{
        if(!halfEdges.contains(face.halfEdge)) {
          if(verbose) println("External reference in face "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }        
        }
        else
        {
          if(face.halfEdge.face!=null){
            if(face.halfEdge.face!=face) {
              if(verbose) println("Wrong reference in face "+i+".");
              if(force){
                result=false;
              }
              else{
                return false;
              }
            }
          }
        }
      }
    }

    if(verbose) println("Checking vertex ("+ vertices.size()+") properties");

    for(int i=0;i<vertices.size();i++){
      HE_Vertex v=(HE_Vertex)vertices.get(i);   
      if(v.halfEdge==null) {
        if(verbose) println("Null reference in vertex  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else
      {
        if(!halfEdges.contains(v.halfEdge)) {
          if(verbose) println("External reference in vertex  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
        if(v.halfEdge.vert!=null){
          if(v.halfEdge.vert!=v){ 
            if(verbose) println("Wrong reference in vertex  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }
          }
        }
      }
    }

    if(verbose) println("Checking edge ("+ edges.size()+") properties");

    for(int i=0;i<edges.size();i++){
      HE_Edge e=(HE_Edge)edges.get(i);   
      if(e.halfEdge==null) {
        if(verbose) println("Null reference in edge  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else{
        if(!halfEdges.contains(e.halfEdge)) {
          if(verbose) println("External reference in edge  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
        if(e.halfEdge.edge!=null){
          if(e.halfEdge.edge!=e) {
            if(verbose) println("Wrong reference in edge  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }
          }
        }
      }
    }

    if(verbose) println("Checking half edge ("+ halfEdges.size()+") properties");

    for(int i=0;i<halfEdges.size();i++){
      HE_HalfEdge he=(HE_HalfEdge)halfEdges.get(i);
      if(he.next==null) {
        if(verbose) println("External reference (next) in half edge  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else{
        if(!halfEdges.contains(he.next)) {
          if(verbose) println("External reference (next) in half edge  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
        if((he.face!=null)&&(he.next.face!=null)){
          if(he.face!=he.next.face) {
            if(verbose) println("Incosistent reference (face) in half edge  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }
          }
        }
      }
      if(he.pair==null) {
        if(verbose) println("Null reference (pair) in half edge  "+i+".");  
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }   
      else{
        if(!halfEdges.contains(he.pair)) {
          if(verbose) println("External reference (pair) in half edge  "+i+"."); 
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }   
        if(he.pair.pair==null){
          if(verbose) println("No pair reference back to half edge  "+i+"."); 
        }
        else
        {
          if(he.pair.pair!=he) {
            if(verbose) println("Wrong pair reference back to half edge  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }
          }
        }
        if((he.edge!=null)&&(he.pair.edge!=null)){
          if(he.edge!=he.pair.edge){
            if(verbose) println("Inconsistent reference (edge) in half edge  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }

          }
        }
      }

      if((he.next!=null)&&(he.pair!=null)){
        if((he.next.vert!=null)&&(he.pair.vert!=null)){
          if(he.next.vert!=he.pair.vert) {
            if(verbose) println("Inconsistent reference (pair)/(next) in half edge  "+i+".");
            if(force){
              result=false;
            }
            else{
              return false;
            }
          }
        }
      }
      if(he.face==null) {
        if(verbose) println("Null reference (face) in half edge  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else{
        if(!faces.contains(he.face)) {
          if(verbose) println("External reference (face) in half edge  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
      }

      if(he.vert==null) {
        if(verbose) println("Null reference (vert) in half edge  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else{
        if(!vertices.contains(he.vert)) {
          if(verbose) println("External reference (vert) in half edge  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
      }
      if(he.edge==null) {
        if(verbose) println("Null reference (edge) in half edge  "+i+".");
        if(force){
          result=false;
        }
        else{
          return false;
        }
      }
      else{
        if(!edges.contains(he.edge)) {
          if(verbose) println("External reference (edge) in half edge  "+i+".");
          if(force){
            result=false;
          }
          else{
            return false;
          }
        }
      }
    }
    if(verbose) println("Validation complete!");
    return result;
  }






}

class HE_HalfEdge{
  int id;
  HE_Vertex vert;
  HE_HalfEdge pair;
  HE_Face face;
  HE_HalfEdge next;
  HE_Edge edge;

  HE_HalfEdge (){
  }
}


class HE_Edge{
  int id;
  boolean visible;
  HE_HalfEdge halfEdge;

  HE_Edge (){
  }

  float intersectionWithPlane(final Plane P, float cutoff){
    PVector p0=halfEdge.vert;
    PVector p1=halfEdge.pair.vert;
    float denom=P.A*(p0.x-p1.x)+P.B*(p0.y-p1.y)+P.C*(p0.z-p1.z);
    if (abs(denom)<cutoff){// edge parallel to plane
      if((P.distanceToPlane(p0)<cutoff)&&(P.distanceToPlane(p1)<cutoff)){// edge in plane
        return -2f;//special case: both endpoints can be considered as intersections
      }
      else{
        return -1f;//no intersection
      }
    }

    float u=(P.A*p0.x+P.B*p0.y+P.C*p0.z+P.D)/denom;
    if ((u<-cutoff)||(u>1f+cutoff)) return -1f;// intersection beyond endpoints
    if(abs(u)<=cutoff) return 0f;// intersection at startpoint
    if((u>=1f-cutoff)&&(u<=1f+cutoff))return 1f;//intersection at endpoint

    return u;//intersection on edge
  }
}




class HE_Vertex extends PVector{
  int id;
  HE_HalfEdge halfEdge;

  HE_Vertex(float x, float y, float z, int id){
    super(x,y,z); 
    this.id=id;
  }
  
  HE_Vertex(PVector v, int id){
    super(v.x,v.y,v.z); 
    this.id=id;
  }


  HE_Vertex get(){
    return new HE_Vertex(x,y,z,id);   

  }

  void projectToPlane(Plane P){
    float d=P.A*x+P.B*y+P.C*z+P.D;
    float x0=x-P.A*d;
    float y0=y-P.B*d;
    float z0=z-P.C*d;
    set(x0,y0,z0);

  }

  void interpolateSelf(HE_Vertex v1, HE_Vertex v2, float f){
    x=v1.x+f*(v2.x-v1.x);
    y=v1.y+f*(v2.y-v1.y);
    z=v1.z+f*(v2.z-v1.z);   

  }

}

class HE_Face{
  int id;
  boolean visible;
  boolean fixed;

  HE_HalfEdge halfEdge; 

  HE_Face(){

  }
}


float sign(float v){
  return (v<0)?-1:((v>0)?1:0);
}



boolean between(float vt, float v1, float v2){
  return(vt<max(v1,v2))&&(vt>min(v1,v2));

}