//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;


  //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;
        /*if(halfEdge.next.next.next==halfEdge){
         beginShape();
         do{              
         HE_Vertex v=halfEdge.vert; 
         vertex(v.x,v.y,v.z);
         
         halfEdge= halfEdge.next;
         } 
         while(halfEdge!=face.halfEdge);
         endShape(CLOSE);    
         }
         else{*/
        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(){
    Iterator vertexItr = vertices.iterator();
    while(vertexItr.hasNext()){
      HE_Vertex v=(HE_Vertex)vertexItr.next();
      pushMatrix();
      translate(v.x,v.y,v.z);  
      box(10);
      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(){
    Iterator eItr = edges.iterator();
    while(eItr.hasNext()){
      HE_Edge e=(HE_Edge)eItr.next();
      if(e.visible)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){
    HE_Mesh result=new HE_Mesh();
    for(int i=0;i<simpleVertices.length;i++){  
      result.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;
      result.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)result.vertices.get(simpleFaces[i][j]);
        if(he.vert.halfEdge==null) he.vert.halfEdge=he;  
      }
      result.cycleHalfEdges(faceEdges,false);    
      result.halfEdges.addAll(faceEdges);
    }
    result.pairHalfEdges(result.halfEdges, result.edges);
    for(int i=0;i<result.edges.size();i++){
      ((HE_Edge)result.edges.get(i)).visible=true;
    }
    return result;
  }

  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;
  }

  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(result.halfEdges,result.edges);
    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 getBox(PVector c, float W, float H, float D){


    float[][] vertices={
      {
        c.x+W,c.y+H,c.z+D                                              }
      ,
      {
        c.x-W,c.y+H,c.z+D                                              }
      ,
      {
        c.x-W,c.y+H,c.z-D                                              }
      ,
      {
        c.x+W,c.y+H,c.z-D                                              }
      ,
      {
        c.x+W,c.y-H,c.z+D                                              }
      ,
      {
        c.x-W,c.y-H,c.z+D                                              }
      ,
      {
        c.x-W,c.y-H,c.z-D                                              }
      ,
      {
        c.x+W,c.y-H,c.z-D                                              }
    };
    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,2,1                                          }
      ,
      {
        0,3,2                                          }
      ,
      {
        0,1,3                                          }
      ,
      {
        1,2,3                                          }

    };

    return buildMesh(vertices, faces);


  }

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

  HE_Mesh extremize(float f){
    HE_Mesh result=get();
    result.setBodyCenter();
    float dmean=0f;
    for(int i=0;i<result.vertices.size();i++){
      HE_Vertex v=(HE_Vertex)result.vertices.get(i);
      dmean+=dist(v.x,v.y,v.z,result.bodyCenter.x,result.bodyCenter.y,result.bodyCenter.z);
    }
    dmean/=result.vertices.size();
    for(int i=0;i<result.vertices.size();i++){
      HE_Vertex v=(HE_Vertex)result.vertices.get(i);
      float d=dist(v.x,v.y,v.z,result.bodyCenter.x,result.bodyCenter.y,result.bodyCenter.z);
      float ff=1+f*(d-dmean)/dmean;
      v.x=result.bodyCenter.x+ff*(v.x-result.bodyCenter.x);
      v.y=result.bodyCenter.y+ff*(v.y-result.bodyCenter.y);
      v.z=result.bodyCenter.z+ff*(v.z-result.bodyCenter.z);
    }
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    return result;
  }

  HE_Mesh spherify(float r, float f){
    HE_Mesh result=get();
    result.setBodyCenter();
    for(int i=0;i<result.vertices.size();i++){
      HE_Vertex v=(HE_Vertex)result.vertices.get(i);
      float d=dist(v.x,v.y,v.z,result.bodyCenter.x,result.bodyCenter.y,result.bodyCenter.z);
      float ff=((1-f)*d+f*r)/d;
      v.x=result.bodyCenter.x+ff*(v.x-result.bodyCenter.x);
      v.y=result.bodyCenter.y+ff*(v.y-result.bodyCenter.y);
      v.z=result.bodyCenter.z+ff*(v.z-result.bodyCenter.z);
    }
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    return result;
  }





  //GEOMETRIC OPERATIONS 
  void 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);

  }

  void triangulate(){
    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;
    }
    bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
    println(faces.size());
  }

  HE_Mesh insertPointInEdge(int id, float f){
    if(!indicesUpdated) reindex();
    HE_Mesh result=get();
    HE_Edge origE=(HE_Edge)result.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);
      result.vertices.add(v);      
      HE_HalfEdge newHE1=new HE_HalfEdge();
      result.halfEdges.add(newHE1);
      HE_HalfEdge newHE2=new HE_HalfEdge();
      result.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();
      result.edges.add(newE);
      newE.halfEdge=newHE1;
      newHE1.edge=newE;
      origHE2.edge=newE;
      newHE2.edge=origE;
      origHE1.edge=origE;
      newE.visible=origE.visible;
    }
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    return result;
  }

  void insertPointInEdgeSelf(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;
    }
    bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
  }

  HE_Mesh insertMidEdgePoints(){
    HE_Mesh result=get();
    int n=result.edges.size();
    for(int i=0;i<n;i++){
      result=result.insertPointInEdge(i,0.5f);
    }
    return result;
  }

  void insertMidEdgePointsSelf(){
    int n=edges.size();
    for(int i=0;i<n;i++){
      insertPointInEdgeSelf(i,0.5f);
    }
  }

  HE_Mesh cutSurface(Plane P){
    if(!indicesUpdated) reindex();
    HE_Mesh result=get(); 
    ArrayList verticesOnPlane=new ArrayList();
    ArrayList newVertices=new ArrayList();
    int n=result.edges.size();
    int numOfIntersections=0;
    for(int i=0;i<n;i++){
      HE_Edge e=(HE_Edge)result.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)
      {
        result.insertPointInEdgeSelf(i,u);
        newVertices.add((HE_Vertex)result.vertices.get(result.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<result.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)result.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);
        result.halfEdges.add(nhe1);
        result.halfEdges.add(nhe2);
        result.edges.add(ne);
        newFaces.add(nf);
      }
    }
    result.faces.addAll(newFaces);
    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 cutSurface(ArrayList cutPlanes){
    HE_Mesh result=get();
    for(int k=0;k<cutPlanes.size();k++){
      Plane P = (Plane)cutPlanes.get(k);          
      result=result.cutSurface(P);
    }
    return result;
  }



  HE_Mesh cutMesh(Plane P, PVector c,boolean safe){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_Mesh result=get(); 
    ArrayList verticesOnPlane=new ArrayList();
    ArrayList newVertices=new ArrayList();
    int n=result.edges.size();
    int numOfIntersections=0;
    for(int i=0;i<n;i++){
      HE_Edge e=(HE_Edge)result.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)
      {
        result.insertPointInEdgeSelf(i,u);
        newVertices.add((HE_Vertex)result.vertices.get(result.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<result.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)result.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);
        result.halfEdges.add(nhe1);
        result.halfEdges.add(nhe2);
        result.edges.add(ne);
        if(P.side(c,0.0001f)*P.side(getFaceCenter(nf),0.0001f)>=0f){
          newFaces.add(nf);
        }
        else{
          newFaces.add(f); 
        }
      }
      else{
        if(P.side(c,0.0001f)*P.side((PVector)result.faceCenters.get(i),0.0001f)>=0f)newFaces.add(f);
      }
    }
    result.faces=newFaces;        
    result.cleanUnusedElements();                    
    result.capOpenLoop();
    result.bodyCenterUpdated=false;
    result.faceCentersUpdated=false;
    result.faceNormalsUpdated=false;
    result.vertexNormalsUpdated=false;
    result.edgeCentersUpdated=false;
    result.edgeNormalsUpdated=false;
    result.indicesUpdated=false;    
    return (result.validate(false,false)||(!safe))?result:get(); 
  }


  HE_Mesh cutMesh(ArrayList cutPlanes, PVector c,boolean safe){
    HE_Mesh result=get();
    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);          
      result=result.cutMesh(P,c,safe);
    }
    return result;
  }


  HE_Mesh chamfer(float f){
    float locf=1f-f;
    if(!faceCentersUpdated)setFaceCenters();
    HE_Mesh result=new HE_Mesh();
    HE_Mesh tracker=get();
    for(int i=0;i<faces.size();i++){
      HE_Face origF=(HE_Face)faces.get(i);
      HE_Face trackF=(HE_Face)tracker.faces.get(i);
      PVector fc=(PVector)faceCenters.get(i);
      HE_HalfEdge origHe=origF.halfEdge; 
      HE_HalfEdge trackHe=trackF.halfEdge; 
      do{
        HE_Vertex vo=origHe.vert;
        HE_Vertex vn=new HE_Vertex(fc.x+locf*(vo.x-fc.x),fc.y+locf*(vo.y-fc.y),fc.z+locf*(vo.z-fc.z),-1);  
        if(!result.vertices.contains(vn))result.vertices.add(vn);
        trackHe.vert=vn;
        origHe=origHe.next;
        trackHe=trackHe.next;       
      }
      while(origHe!=origF.halfEdge); 
    } 
    for(int i=0;i<faces.size();i++){
      HE_Face trackF=(HE_Face)tracker.faces.get(i); 
      HE_HalfEdge trackHe=trackF.halfEdge; 
      HE_Face nf=new HE_Face();
      result.faces.add(nf);
      nf.visible=trackF.visible;
      nf.fixed=trackF.fixed;
      ArrayList newHalfEdges=new ArrayList();
      do{
        HE_Vertex v=trackHe.vert;
        HE_HalfEdge he=new HE_HalfEdge();
        newHalfEdges.add(he);
        he.vert=v;
        he.face=nf;
        if(nf.halfEdge==null) nf.halfEdge=he;
        trackHe=trackHe.next;       
      }
      while(trackHe!=trackF.halfEdge); 
      result.cycleHalfEdges(newHalfEdges,false);
      result.halfEdges.addAll(newHalfEdges);  
    }     
    for(int i=0;i<vertices.size();i++){
      HE_Vertex trackV=(HE_Vertex)tracker.vertices.get(i);
      HE_HalfEdge trackHe=trackV.halfEdge; 
      HE_Face nf=new HE_Face();
      result.faces.add(nf);
      nf.visible=true;
      nf.fixed=isFixed(trackV); 
      ArrayList newHalfEdges=new ArrayList();
      do{
        HE_Vertex v=trackHe.vert;
        HE_HalfEdge he=new HE_HalfEdge();
        newHalfEdges.add(he);
        he.vert=v;
        he.face=nf;      
        if(nf.halfEdge==null) nf.halfEdge=he;
        if(v.halfEdge==null) v.halfEdge=he;
        trackHe=trackHe.pair.next;       
      }
      while(trackHe!=trackV.halfEdge); 
      result.cycleHalfEdges(newHalfEdges,true);
      result.halfEdges.addAll(newHalfEdges);    
    }
    for(int i=0;i<edges.size();i++){
      HE_Edge trackE=(HE_Edge)tracker.edges.get(i);
      HE_HalfEdge trackHe1=trackE.halfEdge; 
      HE_HalfEdge trackHe2=trackE.halfEdge.pair; 
      HE_Face nf=new HE_Face();
      result.faces.add(nf);
      nf.visible=true;
      nf.fixed=isFixed(trackE); 
      HE_Vertex v4=trackHe1.vert;
      HE_Vertex v3=trackHe1.next.vert;
      HE_Vertex v2=trackHe2.vert;
      HE_Vertex v1=trackHe2.next.vert;
      HE_HalfEdge he1=new HE_HalfEdge();
      HE_HalfEdge he2=new HE_HalfEdge();
      HE_HalfEdge he3=new HE_HalfEdge();
      HE_HalfEdge he4=new HE_HalfEdge();
      nf.halfEdge=he1;
      he1.vert=v1;
      he1.face=nf;
      he1.next=he2;
      he2.vert=v2;
      he2.face=nf;
      he2.next=he3;
      he3.vert=v3;
      he3.face=nf;
      he3.next=he4;
      he4.vert=v4;
      he4.face=nf;
      he4.next=he1;
      result.halfEdges.add(he1);
      result.halfEdges.add(he2);
      result.halfEdges.add(he3);
      result.halfEdges.add(he4);
    }
    pairHalfEdges(result.halfEdges,result.edges);
    for(int i=0;i<result.edges.size();i++){
      HE_Edge e=(HE_Edge)result.edges.get(i);
      if((e.halfEdge.face.visible)||(e.halfEdge.pair.face.visible)) e.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;    
  }
  
  
  void 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(halfEdges,edges);  
    
   bodyCenterUpdated=false;
    faceCentersUpdated=false;
    faceNormalsUpdated=false;
    vertexNormalsUpdated=false;
    edgeCentersUpdated=false;
    edgeNormalsUpdated=false;
    indicesUpdated=false;
    
  }
  

   //SUBDIVISION

  HE_Mesh CatmullClark(int rep){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_Mesh result=get(); 
    for(int r=0;r<rep;r++){
      ArrayList verticesToAdjust=new ArrayList();
      verticesToAdjust.addAll(result.vertices);
      if(!result.indicesUpdated) result.reindex();
      if(!result.faceCentersUpdated) result.setFaceCenters();
      int n= result.edges.size();
      for(int i=0;i<n;i++){
        HE_Edge origE=(HE_Edge)result.edges.get(i);
        HE_HalfEdge origHE1= origE.halfEdge;
        HE_HalfEdge origHE2= origHE1.pair;
        HE_Vertex v=new HE_Vertex(0,0,0,-1);
        HE_Vertex v1=origHE1.vert;
        HE_Vertex v2=origHE2.vert;
        HE_Vertex v3=(HE_Vertex)result.faceCenters.get(origHE1.face.id);
        HE_Vertex v4=(HE_Vertex)result.faceCenters.get(origHE2.face.id);    
        if(isFixed(origE)){
          v.x=0.5f*(v1.x+v2.x);
          v.y=0.5f*(v1.y+v2.y);
          v.z=0.5f*(v1.z+v2.z);
        }
        else{
          v.x=0.25f*(v1.x+v2.x+v3.x+v4.x);
          v.y=0.25f*(v1.y+v2.y+v3.y+v4.y);
          v.z=0.25f*(v1.z+v2.z+v3.z+v4.z);
        }
        result.vertices.add(v);
        HE_Edge newE= new HE_Edge();
        newE.visible=origE.visible;
        result.edges.add(newE);
        HE_HalfEdge newHE1=new HE_HalfEdge();
        result.halfEdges.add(newHE1);
        HE_HalfEdge newHE2=new HE_HalfEdge();
        result.halfEdges.add(newHE2);
        newHE1.vert=v;    
        newHE1.next=origHE1.next;
        newHE1.face=origHE1.face;
        newHE1.pair=origHE2;
        newHE1.edge=newE;
        origHE1.next=newHE1;
        origHE1.pair=newHE2;        
        newHE2.vert=v;
        newHE2.next=origHE2.next;
        newHE2.face=origHE2.face;
        newHE2.pair=origHE1;
        newHE2.edge=origE;
        origHE2.next=newHE2;
        origHE2.pair=newHE1;
        origHE2.edge=newE;
        newE.halfEdge=newHE1;
        v.halfEdge=newHE1;
      }
      int tot=0;
      for(int i=0;i<verticesToAdjust.size();i++){
        HE_Vertex P=(HE_Vertex)verticesToAdjust.get(i);
        if(!isFixed(P)){
          int N=0;
          float Fx=0;
          float Fy=0;
          float Fz=0;
          float Rx=0;
          float Ry=0;
          float Rz=0;
          HE_HalfEdge he=P.halfEdge;
          do{
            HE_Face face=he.face;
            HE_Vertex fv=(HE_Vertex)result.faceCenters.get(he.face.id);
            Fx+=fv.x;
            Fy+=fv.y;
            Fz+=fv.z;
            N++;
            HE_Vertex ev1=he.vert;
            HE_Vertex ev2=he.pair.vert;
            Rx+=0.5f*(ev1.x+ev2.x);
            Ry+=0.5f*(ev1.y+ev2.y);
            Rz+=0.5f*(ev1.z+ev2.z);
            he=he.pair;
            he=he.next;
          }
          while(he!=P.halfEdge);
          P.x=(Fx/N+2*Rx/N+(N-3)*P.x)/N;
          P.y=(Fy/N+2*Ry/N+(N-3)*P.y)/N;
          P.z=(Fz/N+2*Rz/N+(N-3)*P.z)/N;
        }
      }
      ArrayList newFaces=new ArrayList();
      for(int i=0;i<result.faces.size();i++){
        HE_Face face=(HE_Face)result.faces.get(i);
        HE_HalfEdge startHE = face.halfEdge;
        while(startHE.vert.id!=-1) startHE=startHE.next;
        HE_HalfEdge origHE1=startHE;
        do{
          HE_Face newFace=new HE_Face();       
          newFace.visible=face.visible;
          newFaces.add(newFace);
          newFace.fixed=face.fixed;
          newFace.halfEdge=origHE1;
          HE_HalfEdge origHE2=origHE1.next;
          HE_HalfEdge origHE3=origHE2.next;
          HE_HalfEdge newHE1=new HE_HalfEdge();
          HE_HalfEdge newHE2=new HE_HalfEdge();    
          result.halfEdges.add(newHE1);
          result.halfEdges.add(newHE2);
          origHE2.next=newHE1;
          newHE1.next=newHE2;
          newHE2.next=origHE1;
          newHE1.vert=origHE3.vert;
          HE_Vertex fv=(HE_Vertex)result.faceCenters.get(origHE1.face.id);
          newHE2.vert=fv;
          if(fv.halfEdge==null) fv.halfEdge=newHE2;
          newHE1.face=newFace;
          newHE2.face=newFace;
          origHE1.face=newFace;
          origHE2.face=newFace;
          origHE1=origHE3;
        }
        while(origHE1!=startHE);
      }
      pairHalfEdges(result.halfEdges,result.edges);
      result.faces=newFaces;
      result.vertices.addAll(result.faceCenters);  
      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 CatmullClark(){
    return CatmullClark(1);  
  }

  HE_Mesh quadSubdivision(int rep){
    if(!indicesUpdated) reindex();
    if(!faceCentersUpdated) setFaceCenters();
    HE_Mesh result=get();
    for(int r=0;r<rep;r++){
      if(!result.indicesUpdated) result.reindex();
      if(!result.faceCentersUpdated) result.setFaceCenters();
      int n= result.edges.size();
      for(int i=0;i<n;i++){
        HE_Edge origE=(HE_Edge)result.edges.get(i);
        HE_HalfEdge origHE1= origE.halfEdge;
        HE_HalfEdge origHE2= origHE1.pair;
        HE_Vertex v=new HE_Vertex(0,0,0,-1);
        HE_Vertex v1=origHE1.vert;
        HE_Vertex v2=origHE2.vert;
        v.x=0.5*(v1.x+v2.x);
        v.y=0.5*(v1.y+v2.y);
        v.z=0.5*(v1.z+v2.z);
        result.vertices.add(v);
        HE_Edge newE= new HE_Edge();
        newE.visible=origE.visible;
        result.edges.add(newE);
        HE_HalfEdge newHE1=new HE_HalfEdge();
        result.halfEdges.add(newHE1);
        HE_HalfEdge newHE2=new HE_HalfEdge();
        result.halfEdges.add(newHE2);
        newHE1.vert=v;    
        newHE1.next=origHE1.next;
        newHE1.face=origHE1.face;
        newHE1.pair=origHE2;
        newHE1.edge=newE;
        origHE1.next=newHE1;
        origHE1.pair=newHE2;        
        newHE2.vert=v;
        newHE2.next=origHE2.next;
        newHE2.face=origHE2.face;
        newHE2.pair=origHE1;
        newHE2.edge=origE;
        origHE2.next=newHE2;
        origHE2.pair=newHE1;
        origHE2.edge=newE;
        newE.halfEdge=newHE1;
        v.halfEdge=newHE1;
      }
      ArrayList newFaces=new ArrayList();
      for(int i=0;i<result.faces.size();i++){
        HE_Face face=(HE_Face)result.faces.get(i);
        HE_HalfEdge startHE = face.halfEdge;
        while(startHE.vert.id!=-1) startHE=startHE.next;
        HE_HalfEdge origHE1=startHE;
        do{
          HE_Face newFace=new HE_Face();
          newFace.visible=face.visible;       
          newFaces.add(newFace);
          newFace.halfEdge=origHE1;
          HE_HalfEdge origHE2=origHE1.next;
          HE_HalfEdge origHE3=origHE2.next;
          HE_HalfEdge newHE1=new HE_HalfEdge();
          HE_HalfEdge newHE2=new HE_HalfEdge();    
          result.halfEdges.add(newHE1);
          result.halfEdges.add(newHE2);
          origHE2.next=newHE1;
          newHE1.next=newHE2;
          newHE2.next=origHE1;
          newHE1.vert=origHE3.vert;
          HE_Vertex fv=(HE_Vertex)result.faceCenters.get(origHE1.face.id);
          newHE2.vert=fv;
          if(fv.halfEdge==null) fv.halfEdge=newHE2;
          newHE1.face=newFace;
          newHE2.face=newFace;
          origHE1.face=newFace;
          origHE2.face=newFace;
          origHE1=origHE3;
        }
        while(origHE1!=startHE);
      }
      result.pairHalfEdges(result.halfEdges,result.edges);
      result.faces=newFaces;
      result.vertices.addAll(result.faceCenters);
      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 quadSubdivision(){
    return quadSubdivision(1);
  }




  HE_Mesh DooSabin(){
    return chamfer(0.5f);
  }   

  HE_Mesh DooSabin(int n){
    HE_Mesh result=get();
    for(int i=0;i<n;i++){
      result=result.chamfer(0.5f);
    }
    return result;
  }

  //HELPER FUNCTIONS

  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 srcHalfEdges, ArrayList tgtEdges){
    int n=srcHalfEdges.size();
    for(int i=0;i<n;i++){
      HE_HalfEdge he=(HE_HalfEdge)srcHalfEdges.get(i);
      if(he.pair==null){
        for(int j=0;j<n;j++){
          if(i!=j){
            HE_HalfEdge he2=(HE_HalfEdge)srcHalfEdges.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;
              tgtEdges.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(){
    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=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);
    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 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 

  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 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);

  }

}

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);
}