WebSVN – gelsvn – Diff – /branches/GEL_github/src/HMesh/Manifold.cpp

 /* ----------------------------------------------------------------------- *
  * This file is part of GEL, http://www.imm.dtu.dk/GEL
  * Copyright (C) the authors and DTU Informatics
  * For license and list of authors, see ../../doc/intro.pdf
  * ----------------------------------------------------------------------- */
 #include <iterator>
 #include "Manifold.h"
 #include <iostream>
 #include <vector>
-#include <map>
+#include <unordered_map>
 #include <iterator>
 #include "../Geometry/TriMesh.h"
 namespace HMesh
 {
     using namespace std;
     using namespace Geometry;
     using namespace CGLA;
     namespace
     {
         /************************************************************
 		 * Edgekeys and Edges for halfedge identification during build
 		 ************************************************************/
         class EdgeKey
         {
         public:
             EdgeKey(VertexID va, VertexID vb)
             {
                 if(va < vb){
                     v0 = va;
                     v1 = vb;
                 }
                 else{
                     v0 = vb;
                     v1 = va;
                 }
             }
             bool operator<(const EdgeKey& k2) const
             {
                 if(v0 < k2.v0){
                     return true;
                 }
                 else if( k2.v0 < v0){
                     return false;
                 }
                 else{
                     return v1 < k2.v1;
                 }
             }
-        private:
-            VertexID v0;
+            bool operator==(const EdgeKey& k2) const {return k2.v0 == v0 && k2.v1==v1;}
-            VertexID v1;
-        };
+            size_t hash() const {return v0.get_index()*3125*49+v1.get_index()*3125+7;}
+        private:
-        struct Edge
+            VertexID v0;
-        {
+            VertexID v1;
-            HalfEdgeID h0;
+        };
-            HalfEdgeID h1;
-            int count;
+        struct Edge
-            Edge() : count(0){}
+        {
-        };
+            HalfEdgeID h0;
-    }
+            HalfEdgeID h1;
+            int count;
-    /*********************************************
+            Edge() : count(0){}
-	 * Public functions
+        };
-	 *********************************************/
+    }
-    void Manifold::build(const TriMesh& mesh)
-    {
+    /*********************************************
-        // A vector of 3's - used to tell build how many indices each face consists of
+	 * Public functions
-        vector<int> faces(mesh.geometry.no_faces(), 3);
+	 *********************************************/
+    void Manifold::build(const TriMesh& mesh)
-        build_template( static_cast<size_t>(mesh.geometry.no_vertices()),
+    {
-					   reinterpret_cast<const float*>(&mesh.geometry.vertex(0)),
+        // A vector of 3's - used to tell build how many indices each face consists of
-					   static_cast<size_t>(faces.size()),
+        vector<int> faces(mesh.geometry.no_faces(), 3);
-					   static_cast<const int*>(&faces[0]),
-					   reinterpret_cast<const int*>(&mesh.geometry.face(0)));
+        build_template( static_cast<size_t>(mesh.geometry.no_vertices()),
-    }
+					   reinterpret_cast<const float*>(&mesh.geometry.vertex(0)),
+					   static_cast<size_t>(faces.size()),
-    void Manifold::build(   size_t no_vertices,
+					   static_cast<const int*>(&faces[0]),
-						 const float* vertvec,
+					   reinterpret_cast<const int*>(&mesh.geometry.face(0)));
-						 size_t no_faces,
+    }
-						 const int* facevec,
-						 const int* indices)
+    void Manifold::build(   size_t no_vertices,
-    {
+						 const float* vertvec,
-        build_template(no_vertices, vertvec, no_faces, facevec, indices);
+						 size_t no_faces,
-    }
+						 const int* facevec,
+						 const int* indices)
-    void Manifold::build(   size_t no_vertices,
+    {
-						 const double* vertvec,
+        build_template(no_vertices, vertvec, no_faces, facevec, indices);
-						 size_t no_faces,
+    }
-						 const int* facevec,
-						 const int* indices)
+    void Manifold::build(   size_t no_vertices,
-    {
+						 const double* vertvec,
-        build_template(no_vertices, vertvec, no_faces, facevec, indices);
+						 size_t no_faces,
-    }
+						 const int* facevec,
+						 const int* indices)
-    FaceID Manifold::add_face(std::vector<Manifold::Vec> points)
+    {
-    {
+        build_template(no_vertices, vertvec, no_faces, facevec, indices);
-        int F = points.size();
+    }
-        vector<int> indices;
-        for(size_t i=0;i<points.size(); ++i)
+    FaceID Manifold::add_face(std::vector<Manifold::Vec> points)
-            indices.push_back(i);
+    {
-        FaceID fid = *faces_end();
+        int F = points.size();
-        build(points.size(), reinterpret_cast<double*>(&points[0]), 1, &F, &indices[0]);
+        vector<int> indices;
-        return fid;
+        for(size_t i=0;i<points.size(); ++i)
-    }
+            indices.push_back(i);
+        FaceID fid = *faces_end();
-    bool Manifold::remove_face(FaceID fid)
+        build(points.size(), reinterpret_cast<double*>(&points[0]), 1, &F, &indices[0]);
-    {
+        return fid;
-        if(!in_use(fid))
+    }
-            return false;
+    bool Manifold::remove_face(FaceID fid)
-        HalfEdgeID he = kernel.last(fid);
+    {
-        HalfEdgeID last = he;
+        if(!in_use(fid))
+            return false;
-        vector<HalfEdgeID> halfedges;
-        do
+        HalfEdgeID he = kernel.last(fid);
-        {
+        HalfEdgeID last = he;
-            halfedges.push_back(he);
-            kernel.set_face(he, InvalidFaceID);
+        vector<HalfEdgeID> halfedges;
-            he = kernel.next(he);
+        do
-        }
+        {
-        while(he != last);
+            halfedges.push_back(he);
+            kernel.set_face(he, InvalidFaceID);
-        vector<HalfEdgeID> halfedges_garbage;
+            he = kernel.next(he);
-        vector<VertexID> vertices_garbage;
+        }
-        for(size_t i=0;i<halfedges.size(); ++i)
+        while(he != last);
-        {
-            HalfEdgeID h = halfedges[i];
+        vector<HalfEdgeID> halfedges_garbage;
-            HalfEdgeID hopp = kernel.opp(h);
+        vector<VertexID> vertices_garbage;
-            if(kernel.face(hopp) == InvalidFaceID)
+        for(size_t i=0;i<halfedges.size(); ++i)
-            {
+        {
-                halfedges_garbage.push_back(h);
+            HalfEdgeID h = halfedges[i];
-                VertexID v0 = kernel.vert(hopp);
+            HalfEdgeID hopp = kernel.opp(h);
-                if(valency(*this, v0) <= 2 )
+            if(kernel.face(hopp) == InvalidFaceID)
-                    vertices_garbage.push_back(v0);
+            {
-                else {
+                halfedges_garbage.push_back(h);
-                    link(kernel.prev(h), kernel.next(hopp));
+                VertexID v0 = kernel.vert(hopp);
-                    kernel.set_out(v0, kernel.opp(kernel.prev(h)));
+                if(valency(*this, v0) <= 2 )
-                }
+                    vertices_garbage.push_back(v0);
-                VertexID v1 = kernel.vert(h);
+                else {
-                if(valency(*this, v1)>2)
+                    link(kernel.prev(h), kernel.next(hopp));
-                {
+                    kernel.set_out(v0, kernel.opp(kernel.prev(h)));
-                    link(kernel.prev(hopp),kernel.next(h));
+                }
-                    kernel.set_out(v1, kernel.opp(kernel.prev(hopp)));
+                VertexID v1 = kernel.vert(h);
-                }
+                if(valency(*this, v1)>2)
+                {
-            }
+                    link(kernel.prev(hopp),kernel.next(h));
-        }
+                    kernel.set_out(v1, kernel.opp(kernel.prev(hopp)));
-        for(size_t i=0;i<halfedges_garbage.size();++i){
+                }
-            kernel.remove_halfedge(kernel.opp(halfedges_garbage[i]));
-            kernel.remove_halfedge(halfedges_garbage[i]);
+            }
         }
-        for(size_t i=0;i<vertices_garbage.size(); ++i)
+        for(size_t i=0;i<halfedges_garbage.size();++i){
-            kernel.remove_vertex(vertices_garbage[i]);
+            kernel.remove_halfedge(kernel.opp(halfedges_garbage[i]));
+            kernel.remove_halfedge(halfedges_garbage[i]);
-        kernel.remove_face(fid);
+        }
-        return true;
+        for(size_t i=0;i<vertices_garbage.size(); ++i)
-    }
+            kernel.remove_vertex(vertices_garbage[i]);
+        kernel.remove_face(fid);
-    bool Manifold::remove_edge(HalfEdgeID hid)
+        return true;
-    {
+    }
-        if(!in_use(hid))
-            return false;
+    bool Manifold::remove_edge(HalfEdgeID hid)
-        FaceID f0 = kernel.face(hid);
+    {
-        FaceID f1 = kernel.face(kernel.opp(hid));
+        if(!in_use(hid))
+            return false;
-        remove_face(f0);
-        remove_face(f1);
+        FaceID f0 = kernel.face(hid);
+        FaceID f1 = kernel.face(kernel.opp(hid));
-        return true;
-    }
+        remove_face(f0);
+        remove_face(f1);
+        return true;
-    bool Manifold::remove_vertex(VertexID vid)
+    }
-    {
-        if(!in_use(vid))
-            return false;
+    bool Manifold::remove_vertex(VertexID vid)
-        vector<FaceID> faces;
+    {
-        int N = circulate_vertex_ccw(*this, vid, (std::function<void(FaceID)>)[&](FaceID f) {
+        if(!in_use(vid))
-            faces.push_back(f);
+            return false;
-        });
-        for(size_t i=0;i<N;++i)
+        vector<FaceID> faces;
-            remove_face(faces[i]);
+        int N = circulate_vertex_ccw(*this, vid, (std::function<void(FaceID)>)[&](FaceID f) {
+            faces.push_back(f);
-        return true;
+        });
-    }
+        for(size_t i=0;i<N;++i)
+            remove_face(faces[i]);
+        return true;
-    void Manifold::collapse_edge(HalfEdgeID h, bool avg_vertices)
+    }
-    {
-        HalfEdgeID ho = kernel.opp(h);
-        VertexID hv = kernel.vert(h);
-        VertexID hov = kernel.vert(ho);
+    void Manifold::collapse_edge(HalfEdgeID h, bool avg_vertices)
-        HalfEdgeID hn = kernel.next(h);
+    {
-        HalfEdgeID hp = kernel.prev(h);
+        HalfEdgeID ho = kernel.opp(h);
-        HalfEdgeID hon = kernel.next(ho);
+        VertexID hv = kernel.vert(h);
-        HalfEdgeID hop = kernel.prev(ho);
+        VertexID hov = kernel.vert(ho);
-		FaceID f = kernel.face(h);
+        HalfEdgeID hn = kernel.next(h);
-		FaceID fo = kernel.face(ho);
+        HalfEdgeID hp = kernel.prev(h);
+        HalfEdgeID hon = kernel.next(ho);
-        // average the vertex positions
+        HalfEdgeID hop = kernel.prev(ho);
-        pos(hv) = avg_vertices ? (0.5f * (pos(hov) + pos(hv))) : pos(hv);
+		FaceID f = kernel.face(h);
+		FaceID fo = kernel.face(ho);
-        // update all halfedges pointing to hov to point to hv, effectively removing hov from all loops
-        HalfEdgeID he = kernel.out(hov);
+        // average the vertex positions
-        HalfEdgeID last = he;
+        pos(hv) = avg_vertices ? (0.5f * (pos(hov) + pos(hv))) : pos(hv);
-		do {
-			assert(kernel.vert(kernel.opp(he)) == hov);
+        // update all halfedges pointing to hov to point to hv, effectively removing hov from all loops
-            kernel.set_vert(kernel.opp(he), hv);
+        HalfEdgeID he = kernel.out(hov);
-            he = kernel.next(kernel.opp(he));
+        HalfEdgeID last = he;
-        }
+		do {
-		while(he != last);
+			assert(kernel.vert(kernel.opp(he)) == hov);
-        kernel.set_out(hv, hn);
+            kernel.set_vert(kernel.opp(he), hv);
+            he = kernel.next(kernel.opp(he));
-        // link hp and hn, effectively removing h from opposite loop
+        }
-        link(hp, hn);
+		while(he != last);
-        // make face owning h own hn instead
+        kernel.set_out(hv, hn);
-        if(kernel.face(h) != InvalidFaceID)
-            kernel.set_last(f, hn);
+        // link hp and hn, effectively removing h from opposite loop
+        link(hp, hn);
-        // link hop and hon, effectively removing h from opposite face loop
+        // make face owning h own hn instead
-        link(hop, hon);
+        if(kernel.face(h) != InvalidFaceID)
-        // make opposite face owning h own hon instead
+            kernel.set_last(f, hn);
-        if(kernel.face(ho) != InvalidFaceID)
-            kernel.set_last(fo, hon);
+        // link hop and hon, effectively removing h from opposite face loop
+        link(hop, hon);
-        // remove the obsolete entities
+        // make opposite face owning h own hon instead
-        kernel.remove_vertex(hov);
+        if(kernel.face(ho) != InvalidFaceID)
-        kernel.remove_halfedge(h);
+            kernel.set_last(fo, hon);
-        kernel.remove_halfedge(ho);
+        // remove the obsolete entities
-        // verify that remaining faces haven't become degenerate because of collapse
+        kernel.remove_vertex(hov);
-        remove_face_if_degenerate(hn);
+        kernel.remove_halfedge(h);
-        remove_face_if_degenerate(hon);
+        kernel.remove_halfedge(ho);
-        // verify that v is sane after collapse
+        // verify that remaining faces haven't become degenerate because of collapse
-        ensure_boundary_consistency(hv);
+        remove_face_if_degenerate(hn);
-    }
+        remove_face_if_degenerate(hon);
-    FaceID Manifold::split_face_by_edge(FaceID f, VertexID v0, VertexID v1)
+        // verify that v is sane after collapse
-    {
+        ensure_boundary_consistency(hv);
-		assert(!connected(*this, v0, v1));
+    }
-        HalfEdgeID he = kernel.last(f);
-        HalfEdgeID last = he;
+    FaceID Manifold::split_face_by_edge(FaceID f, VertexID v0, VertexID v1)
-        int steps = 0;
+    {
-        while(he != last || steps == 0){
+		if(connected(*this, v0, v1))
-            ++steps;
+            return InvalidFaceID;
-            he = kernel.next(he);
-        }
+        HalfEdgeID he = kernel.last(f);
+        HalfEdgeID last = he;
-        // make sure this is not a triangle
+        int steps = 0;
-        assert(steps > 3);
+        while(he != last || steps == 0){
-        // make sure we are not trying to connect a vertex to itself
+            ++steps;
-        assert(v0 != v1);
+            he = kernel.next(he);
+        }
-		HalfEdgeID h0 = kernel.out(v0);
-		for(Walker w = walker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){
+        // make sure this is not a triangle
-			if(w.face() == f){
+        assert(steps > 3);
-				h0 = w.halfedge();
+        // make sure we are not trying to connect a vertex to itself
-				break;
+        assert(v0 != v1);
-			}
-		}
+		HalfEdgeID h0 = kernel.out(v0);
-		assert(kernel.face(h0) != InvalidFaceID);
+		for(Walker w = walker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){
-        assert(kernel.face(h0) == f);
+			if(w.face() == f){
+				h0 = w.halfedge();
-        // the halfedge belonging to f, going out from v0, is denoted h. Move along h until we hit v1.
+				break;
-        // now we have the halfedge which belongs to f and points to v1.
+			}
-        HalfEdgeID h = h0;
+		}
-        while(kernel.vert(h) != v1){
+		assert(kernel.face(h0) != InvalidFaceID);
-            h = kernel.next(h);
+        assert(kernel.face(h0) == f);
-        }
-        assert(h != h0);
+        // the halfedge belonging to f, going out from v0, is denoted h. Move along h until we hit v1.
+        // now we have the halfedge which belongs to f and points to v1.
-        // create a new halfedge ha which connects v1 and v0 closing the first loop.
+        HalfEdgeID h = h0;
-        HalfEdgeID h1 = kernel.next(h);
+        while(kernel.vert(h) != v1){
-        HalfEdgeID ha = kernel.add_halfedge();
+            h = kernel.next(h);
-        link(h, ha);
+        }
-        link(ha, h0);
+        assert(h != h0);
-        kernel.set_face(ha, f);
-        kernel.set_vert(ha, v0);
+        // create a new halfedge ha which connects v1 and v0 closing the first loop.
-        kernel.set_last(f, ha);
+        HalfEdgeID h1 = kernel.next(h);
+        HalfEdgeID ha = kernel.add_halfedge();
-        // create a new face, f2, and set all halfedges in the remaining part of the polygon to point to this face.
+        link(h, ha);
-        h = h1;
+        link(ha, h0);
-        FaceID f2 = kernel.add_face();
+        kernel.set_face(ha, f);
-        while(kernel.vert(h) != v0){
+        kernel.set_vert(ha, v0);
-            kernel.set_face(h, f2);
+        kernel.set_last(f, ha);
-            h = kernel.next(h);
-        }
+        // create a new face, f2, and set all halfedges in the remaining part of the polygon to point to this face.
-        kernel.set_face(h, f2);
+        h = h1;
-        assert(h != h1);
+        FaceID f2 = kernel.add_face();
+        while(kernel.vert(h) != v0){
-        // create a new halfedge hb to connect v0 and v1.
+            kernel.set_face(h, f2);
-        HalfEdgeID hb = kernel.add_halfedge();
+            h = kernel.next(h);
-        link(h, hb);
+        }
-        link(hb, h1);
+        kernel.set_face(h, f2);
-        kernel.set_face(hb, f2);
+        assert(h != h1);
-        kernel.set_vert(hb, v1);
-        kernel.set_last(f2, hb);
+        // create a new halfedge hb to connect v0 and v1.
+        HalfEdgeID hb = kernel.add_halfedge();
-        // complete the operation by gluing the two new halfedges
+        link(h, hb);
-        glue(ha, hb);
+        link(hb, h1);
+        kernel.set_face(hb, f2);
-        // assert sanity of operation
+        kernel.set_vert(hb, v1);
-        assert(kernel.next(kernel.opp(kernel.prev(h1))) == h0);
+        kernel.set_last(f2, hb);
-        assert(kernel.next(kernel.opp(kernel.prev(h0))) == h1);
-        assert(kernel.face(kernel.next(hb)) == f2);
+        // complete the operation by gluing the two new halfedges
-        assert(kernel.face(kernel.next(kernel.next(hb))) == f2);
+        glue(ha, hb);
-        assert(kernel.face(hb) == f2);
+        // assert sanity of operation
-        // return handle to newly created face
+        assert(kernel.next(kernel.opp(kernel.prev(h1))) == h0);
-        return f2;
+        assert(kernel.next(kernel.opp(kernel.prev(h0))) == h1);
-    }
+        assert(kernel.face(kernel.next(hb)) == f2);
+        assert(kernel.face(kernel.next(kernel.next(hb))) == f2);
-    VertexID Manifold::split_face_by_vertex(FaceID f)
+        assert(kernel.face(hb) == f2);
-    {
-        //create the new vertex, with the barycenter of the face as position
+        // return handle to newly created face
-        Manifold::Vec p(0.0f);
+        return f2;
-        HalfEdgeID last_he = kernel.last(f);
+    }
-        HalfEdgeID he = last_he;
-        int steps = 0;
+    VertexID Manifold::split_face_by_vertex(FaceID f)
+    {
-        while(he != last_he || steps == 0){
+        //create the new vertex, with the barycenter of the face as position
-            p += positions[kernel.vert(he)];
+        Manifold::Vec p(0.0f);
-            ++steps;
+        HalfEdgeID last_he = kernel.last(f);
-            he = kernel.next(he);
+        HalfEdgeID he = last_he;
-        }
+        int steps = 0;
-        p /= steps;
+        while(he != last_he || steps == 0){
+            p += positions[kernel.vert(he)];
-        VertexID v = kernel.add_vertex();
+            ++steps;
-        positions[v]  = p;
+            he = kernel.next(he);
+        }
-        //circulate the face, create halfedges and connect to vertex
-        vector<HalfEdgeID> eout;
+        p /= steps;
-        vector<HalfEdgeID> ein;
-        last_he = kernel.last(f);
+        VertexID v = kernel.add_vertex();
-        he = last_he;
+        positions[v]  = p;
-        do{
+        //circulate the face, create halfedges and connect to vertex
-            HalfEdgeID hn = kernel.next(he);
+        vector<HalfEdgeID> eout;
+        vector<HalfEdgeID> ein;
-            HalfEdgeID ho = kernel.add_halfedge();
+        last_he = kernel.last(f);
-            HalfEdgeID hi = kernel.add_halfedge();
+        he = last_he;
-            FaceID fn = kernel.add_face();
+        do{
-            kernel.set_face(hi, fn);
+            HalfEdgeID hn = kernel.next(he);
-            kernel.set_vert(hi, v);
-            kernel.set_face(ho, fn);
+            HalfEdgeID ho = kernel.add_halfedge();
-            kernel.set_vert(ho, kernel.vert(kernel.opp(he)));
+            HalfEdgeID hi = kernel.add_halfedge();
-            kernel.set_face(he, fn);
-            kernel.set_last(fn, ho);
+            FaceID fn = kernel.add_face();
+            kernel.set_face(hi, fn);
-            link(hi, ho);
+            kernel.set_vert(hi, v);
-            link(ho, he);
+            kernel.set_face(ho, fn);
-            link(he, hi);
+            kernel.set_vert(ho, kernel.vert(kernel.opp(he)));
+            kernel.set_face(he, fn);
-            eout.push_back(ho);
+            kernel.set_last(fn, ho);
-            ein.push_back(hi);
+            link(hi, ho);
-            he = hn;
+            link(ho, he);
-        }
+            link(he, hi);
-        while(he != last_he);
+            eout.push_back(ho);
-        // glue new halfedges together
+            ein.push_back(hi);
-        size_t N = eout.size();
-        for(size_t i = 0; i < N; ++i){
+            he = hn;
-            glue(ein[i], eout[(i+1)%N]);
+        }
-        }
+        while(he != last_he);
-        kernel.set_out(v, eout[0]);
+        // glue new halfedges together
-        //remove the now replaced face
+        size_t N = eout.size();
-        kernel.remove_face(f);
+        for(size_t i = 0; i < N; ++i){
+            glue(ein[i], eout[(i+1)%N]);
-        return v;
+        }
-    }
+        kernel.set_out(v, eout[0]);
-    VertexID Manifold::split_edge(HalfEdgeID h)
-    {
+        //remove the now replaced face
-        HalfEdgeID ho = kernel.opp(h);
+        kernel.remove_face(f);
-        VertexID v = kernel.vert(h);
-        VertexID vo = kernel.vert(ho);
+        return v;
+    }
-        //create the new vertex with middle of edge as position and update connectivity
+    VertexID Manifold::split_edge(HalfEdgeID h)
-        VertexID vn = kernel.add_vertex();
+    {
-        positions[vn] = .5f * (positions[v] + positions[vo]);
+        HalfEdgeID ho = kernel.opp(h);
-        kernel.set_out(vn, h);
+        VertexID v = kernel.vert(h);
+        VertexID vo = kernel.vert(ho);
-        //create two necessary halfedges, and update connectivity
-        HalfEdgeID hn = kernel.add_halfedge();
+        //create the new vertex with middle of edge as position and update connectivity
-        HalfEdgeID hno = kernel.add_halfedge();
+        VertexID vn = kernel.add_vertex();
+        positions[vn] = .5f * (positions[v] + positions[vo]);
-        kernel.set_out(vo, hn);
+        kernel.set_out(vn, h);
-        kernel.set_out(v, hno);
+        //create two necessary halfedges, and update connectivity
-        glue(h, hno);
+        HalfEdgeID hn = kernel.add_halfedge();
-        glue(hn, ho);
+        HalfEdgeID hno = kernel.add_halfedge();
-        link(kernel.prev(h), hn);
+        kernel.set_out(vo, hn);
-        link(hn, h);
+        kernel.set_out(v, hno);
-        kernel.set_vert(hn, vn);
+        glue(h, hno);
-        link(kernel.prev(ho), hno);
+        glue(hn, ho);
-        link(hno, ho);
-        kernel.set_vert(hno, vn);
+        link(kernel.prev(h), hn);
+        link(hn, h);
-        // update faces in case of non boundary edge
+        kernel.set_vert(hn, vn);
-        if(kernel.face(h) != InvalidFaceID)
-            kernel.set_last(kernel.face(h), hn);
+        link(kernel.prev(ho), hno);
+        link(hno, ho);
-        if(kernel.face(ho) != InvalidFaceID)
+        kernel.set_vert(hno, vn);
-            kernel.set_last(kernel.face(ho), ho);
+        // update faces in case of non boundary edge
-        // update new edge with faces from existing edge
+        if(kernel.face(h) != InvalidFaceID)
-        kernel.set_face(hn, kernel.face(h));
+            kernel.set_last(kernel.face(h), hn);
-        kernel.set_face(hno, kernel.face(ho));
+        if(kernel.face(ho) != InvalidFaceID)
-        //if split occurs on a boundary, consistency must be ensured.
+            kernel.set_last(kernel.face(ho), ho);
-        ensure_boundary_consistency(vn);
-        ensure_boundary_consistency(v);
+        // update new edge with faces from existing edge
-        ensure_boundary_consistency(vo);
+        kernel.set_face(hn, kernel.face(h));
+        kernel.set_face(hno, kernel.face(ho));
-        return vn;
-    }
+        //if split occurs on a boundary, consistency must be ensured.
+        ensure_boundary_consistency(vn);
-    size_t link_intersection(const Manifold& m, VertexID v0, VertexID v1, vector<VertexID>& lisect)
+        ensure_boundary_consistency(v);
-    {
+        ensure_boundary_consistency(vo);
-        // get the one-ring of v0
-        vector<VertexID> link0;
+        return vn;
-        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID vn) {
+    }
-            link0.push_back(vn);
-        });
+    size_t link_intersection(const Manifold& m, VertexID v0, VertexID v1, vector<VertexID>& lisect)
+    {
-        // get the one-ring of v1
+        // get the one-ring of v0
-        vector<VertexID> link1;
+        vector<VertexID> link0;
-        circulate_vertex_ccw(m, v1, (std::function<void(VertexID)>)[&](VertexID vn) {
+        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID vn) {
-            link1.push_back(vn);
+            link0.push_back(vn);
         });
-        // sort the vertices of the two rings
+        // get the one-ring of v1
-        sort(link0.begin(), link0.end());
+        vector<VertexID> link1;
-        sort(link1.begin(), link1.end());
+        circulate_vertex_ccw(m, v1, (std::function<void(VertexID)>)[&](VertexID vn) {
+            link1.push_back(vn);
-        // get the intersection of the shared vertices from both rings
+        });
-        std::back_insert_iterator<vector<VertexID> > lii(lisect);
-        set_intersection(link0.begin(), link0.end(),
+        // sort the vertices of the two rings
-						 link1.begin(), link1.end(),
+        sort(link0.begin(), link0.end());
-						 lii);
+        sort(link1.begin(), link1.end());
-        return lisect.size();
+        // get the intersection of the shared vertices from both rings
-    }
+        std::back_insert_iterator<vector<VertexID> > lii(lisect);
+        set_intersection(link0.begin(), link0.end(),
-    bool Manifold::stitch_boundary_edges(HalfEdgeID h0, HalfEdgeID h1)
+						 link1.begin(), link1.end(),
-    {
+						 lii);
-        // Cannot stitch an edge with itself
-        if(h0 == h1)
+        return lisect.size();
-            return false;
+    }
-        // Only stitch a pair of boundary edges.
+    bool Manifold::stitch_boundary_edges(HalfEdgeID h0, HalfEdgeID h1)
-        if(kernel.face(h0) == InvalidFaceID && kernel.face(h1) == InvalidFaceID)
+    {
-        {
+        // Cannot stitch an edge with itself
-            HalfEdgeID h0o = kernel.opp(h0);
+        if(h0 == h1)
-            HalfEdgeID h1o = kernel.opp(h1);
+            return false;
-            VertexID v0a = kernel.vert(h0);
-            VertexID v0b = kernel.vert(kernel.opp(h1));
+        // Only stitch a pair of boundary edges.
-            VertexID v1b = kernel.vert(h1);
+        if(kernel.face(h0) == InvalidFaceID && kernel.face(h1) == InvalidFaceID)
-            VertexID v1a = kernel.vert(kernel.opp(h0));
+        {
+            HalfEdgeID h0o = kernel.opp(h0);
-            // Case below implies that h0 and h1 are the same edge with different ID
+            HalfEdgeID h1o = kernel.opp(h1);
-            // That should not happen.
+            VertexID v0a = kernel.vert(h0);
-            assert(!(v0a == v0b && v1a == v1b));
+            VertexID v0b = kernel.vert(kernel.opp(h1));
+            VertexID v1b = kernel.vert(h1);
-            if(connected(*this, v0a, v0b))
+            VertexID v1a = kernel.vert(kernel.opp(h0));
-                return false;
-            if(connected(*this, v1a, v1b))
+            // Case below implies that h0 and h1 are the same edge with different ID
-                return false;
+            // That should not happen.
+            assert(!(v0a == v0b && v1a == v1b));
-            if(v0b != v0a)
+            if(connected(*this, v0a, v0b))
-            {
+                return false;
+            if(connected(*this, v1a, v1b))
-                // Check the link intersection v0a and v0b are welded together
+                return false;
-                // if they share a neighbouring vertex, it will appear twice in the combined
-                // one ring unless it v1a and v1a==v1b
-                vector<VertexID> lisect;
+            if(v0b != v0a)
-                if(link_intersection(*this, v0a, v0b, lisect))
+            {
-                {
-                    vector<VertexID>::iterator iter;
+                // Check the link intersection v0a and v0b are welded together
+                // if they share a neighbouring vertex, it will appear twice in the combined
-                    if(v1a == v1b)
+                // one ring unless it v1a and v1a==v1b
-                    {
+                vector<VertexID> lisect;
-                        iter = find(lisect.begin(), lisect.end(), v1a);
+                if(link_intersection(*this, v0a, v0b, lisect))
-                        if(iter != lisect.end())
+                {
-                            lisect.erase(iter);
+                    vector<VertexID>::iterator iter;
-                    }
-                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h0)));
+                    if(v1a == v1b)
-                    if(iter != lisect.end())
+                    {
-                        lisect.erase(iter);
+                        iter = find(lisect.begin(), lisect.end(), v1a);
-                    if(lisect.size() > 0)
+                        if(iter != lisect.end())
-                        return false;
+                            lisect.erase(iter);
-                }
+                    }
-            }
+                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h0)));
+                    if(iter != lisect.end())
-            if(v1a != v1b)
+                        lisect.erase(iter);
-            {
+                    if(lisect.size() > 0)
-                // Check the same for the other endpoints.
+                        return false;
-                vector<VertexID> lisect;
+                }
-                if(link_intersection(*this, v1a, v1b, lisect))
+            }
-                {
-                    vector<VertexID>::iterator iter;
+            if(v1a != v1b)
+            {
-                    if(v0a == v0b)
+                // Check the same for the other endpoints.
-                    {
+                vector<VertexID> lisect;
-                        iter = find(lisect.begin(), lisect.end(), v1a);
+                if(link_intersection(*this, v1a, v1b, lisect))
-                        if(iter != lisect.end())
+                {
-                            lisect.erase(iter);
+                    vector<VertexID>::iterator iter;
-                    }
-                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h1)));
+                    if(v0a == v0b)
-                    if(iter != lisect.end())
+                    {
-                        lisect.erase(iter);
+                        iter = find(lisect.begin(), lisect.end(), v0a);
-                    if(lisect.size() > 0)
+                        if(iter != lisect.end())
-                        return false;
+                            lisect.erase(iter);
-                }
+                    }
-            }
+                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h1)));
+                    if(iter != lisect.end())
+                        lisect.erase(iter);
-            if(v0b != v0a)
+                    if(lisect.size() > 0)
-                circulate_vertex_ccw(*this, v0b, (std::function<void(Walker&)>)[&](Walker hew) {
+                        return false;
-                    kernel.set_vert(hew.opp().halfedge(), v0a);
+                }
-                });
+            }
-            if(v1b != v1a)
-                circulate_vertex_ccw(*this, v1b, (std::function<void(Walker&)>)[&](Walker hew) {
+            if(v0b != v0a)
-                    kernel.set_vert(hew.opp().halfedge(), v1a);
+                circulate_vertex_ccw(*this, v0b, (std::function<void(Walker&)>)[&](Walker hew) {
-                });
+                    kernel.set_vert(hew.opp().halfedge(), v0a);
+                });
-            if(v0a != v0b)
-            {
+            if(v1b != v1a)
-                HalfEdgeID h1p = kernel.prev(h1);
+                circulate_vertex_ccw(*this, v1b, (std::function<void(Walker&)>)[&](Walker hew) {
-                HalfEdgeID h0n = kernel.next(h0);
+                    kernel.set_vert(hew.opp().halfedge(), v1a);
+                });
-                if(kernel.next(h0n) == h1p)
-                {
+            if(v0a != v0b)
-                    glue(kernel.opp(h0n), kernel.opp(h1p));
+            {
-                    kernel.set_out(kernel.vert(h0n),kernel.opp(h0n));
+                HalfEdgeID h1p = kernel.prev(h1);
-                    kernel.remove_halfedge(h0n);
+                HalfEdgeID h0n = kernel.next(h0);
-                    kernel.remove_halfedge(h1p);
-                }
+                if(kernel.next(h0n) == h1p)
-                else
+                {
-                    link(h1p, h0n);
+                    glue(kernel.opp(h0n), kernel.opp(h1p));
-                kernel.remove_vertex(v0b);
+                    kernel.set_out(kernel.vert(h0n),kernel.opp(h0n));
-            }
+                    kernel.remove_halfedge(h0n);
+                    kernel.remove_halfedge(h1p);
-            if(v1a != v1b)
+                }
-            {
+                else
-                HalfEdgeID h0p = kernel.prev(h0);
+                    link(h1p, h0n);
-                HalfEdgeID h1n = kernel.next(h1);
+                kernel.remove_vertex(v0b);
-                if(kernel.next(h1n) == h0p)
+            }
-                {
-                    glue(kernel.opp(h0p), kernel.opp(h1n));
+            if(v1a != v1b)
-                    kernel.set_out(kernel.vert(h1n), kernel.opp(h1n));
+            {
-                    kernel.remove_halfedge(h0p);
+                HalfEdgeID h0p = kernel.prev(h0);
-                    kernel.remove_halfedge(h1n);
+                HalfEdgeID h1n = kernel.next(h1);
-                }
+                if(kernel.next(h1n) == h0p)
-                else
+                {
-                    link(h0p, h1n);
+                    glue(kernel.opp(h0p), kernel.opp(h1n));
-                kernel.remove_vertex(v1b);
+                    kernel.set_out(kernel.vert(h1n), kernel.opp(h1n));
-            }
+                    kernel.remove_halfedge(h0p);
-            glue(h0o, h1o);
+                    kernel.remove_halfedge(h1n);
+                }
-            kernel.remove_halfedge(h0);
+                else
-            kernel.remove_halfedge(h1);
+                    link(h0p, h1n);
+                kernel.remove_vertex(v1b);
-            kernel.set_out(v1a, h1o);
+            }
-            kernel.set_out(v0a, h0o);
+            glue(h0o, h1o);
-            ensure_boundary_consistency(v1a);
+            kernel.remove_halfedge(h0);
-            ensure_boundary_consistency(v0a);
+            kernel.remove_halfedge(h1);
-            return true;
+            kernel.set_out(v1a, h1o);
-        }
+            kernel.set_out(v0a, h0o);
-        return false;
-    }
+            ensure_boundary_consistency(v1a);
+            ensure_boundary_consistency(v0a);
+            return true;
-    FaceID Manifold::merge_one_ring(VertexID v, float max_loop_length)
+        }
-    {
+        return false;
-        // If the vertex is either not in use or has just
+    }
-        // one incident edge (or less), bail out.
-        int val = valency(*this,v);
-        if(!in_use(v) || val<2)
-            return InvalidFaceID;
+    FaceID Manifold::merge_one_ring(VertexID v, float max_loop_length)
+    {
-        // If the vertex is  a boundary vertex, we preparte the walker so that the
+        // If the vertex is either not in use or has just
-        // first face visited is not the invalid face outside the boundary. If the boundary
+        // one incident edge (or less), bail out.
-        // vertex is adjacent to only one vertex, there is little to do and we bail out.
+        int val = valency(*this,v);
-        bool vertex_is_boundary = false;
+        if(!in_use(v) || val<2)
-        Walker hew = walker(v);
+            return InvalidFaceID;
-        if(boundary(*this, v))
-        {
+        // If the vertex is  a boundary vertex, we preparte the walker so that the
-            if(val==2) return InvalidFaceID;
+        // first face visited is not the invalid face outside the boundary. If the boundary
-            vertex_is_boundary = true;
+        // vertex is adjacent to only one vertex, there is little to do and we bail out.
-            hew = hew.circulate_vertex_ccw();
+        bool vertex_is_boundary = false;
-        }
+        Walker hew = walker(v);
+        if(boundary(*this, v))
-        // Prepare some vectors for taking out the trash: We remove all old faces and all orphaned edges
+        {
-        // and vertices
+            if(val==2) return InvalidFaceID;
-        vector<HalfEdgeID> garbage_halfedges;
+            vertex_is_boundary = true;
-        vector<FaceID> garbage_faces;
+            hew = hew.circulate_vertex_ccw();
-        vector<VertexID> garbage_vertices;
+        }
-        vector<HalfEdgeID> loop; // The halfedges which form the outer loop of the merged one ring.
+        // Prepare some vectors for taking out the trash: We remove all old faces and all orphaned edges
+        // and vertices
-        // Below we loop over all faces adjacent to the vertex and add their halfedges to a big loop
+        vector<HalfEdgeID> garbage_halfedges;
-        // which will form the loop of the new merged face. Below we remove from the loop edges
+        vector<FaceID> garbage_faces;
-        // that appear twice (as each other's opposite).
+        vector<VertexID> garbage_vertices;
-        for(;!hew.full_circle(); hew = hew.circulate_vertex_ccw())
-        {
+        vector<HalfEdgeID> loop; // The halfedges which form the outer loop of the merged one ring.
-            garbage_faces.push_back(hew.face());
-            for(Walker hew2 = walker(hew.halfedge());
+        // Below we loop over all faces adjacent to the vertex and add their halfedges to a big loop
-                !hew2.full_circle(); hew2 = hew2.circulate_face_ccw())
+        // which will form the loop of the new merged face. Below we remove from the loop edges
-                loop.push_back(hew2.halfedge());
+        // that appear twice (as each other's opposite).
-        }
+        for(;!hew.full_circle(); hew = hew.circulate_vertex_ccw())
+        {
+            garbage_faces.push_back(hew.face());
-        // Now merge the loops. We iteratively remove pairs of adjacent halfedges from the loop
+            for(Walker hew2 = walker(hew.halfedge());
-        // if we find that the second is the opposite of the first since this is a degenerate
+                !hew2.full_circle(); hew2 = hew2.circulate_face_ccw())
-        // situation. However, we stop at four remaining halfedges since otherwise the loop degenerates
+                loop.push_back(hew2.halfedge());
-        // to zero after the next step if these four are also pairwise each other's opposites.
+        }
-        int did_work;
-        do
-        {
+        // Now merge the loops. We iteratively remove pairs of adjacent halfedges from the loop
-            did_work = 0;
+        // if we find that the second is the opposite of the first since this is a degenerate
-            vector<HalfEdgeID> loop_tmp(0);
+        // situation. However, we stop at four remaining halfedges since otherwise the loop degenerates
-            for(size_t i=0;i<loop.size();++i)
+        // to zero after the next step if these four are also pairwise each other's opposites.
-                if(walker(loop[i]).opp().halfedge() == loop[(i+1)%loop.size()])
+        int did_work;
-                {
+        do
-                    VertexID vid = walker(loop[i]).vertex();
+        {
-                    if(vid != v)
+            did_work = 0;
-                        garbage_vertices.push_back(walker(loop[i]).vertex());
+            vector<HalfEdgeID> loop_tmp(0);
-                    garbage_halfedges.push_back(loop[i]);
+            for(size_t i=0;i<loop.size();++i)
-                    garbage_halfedges.push_back(loop[(i+1)%loop.size()]);
+                if(walker(loop[i]).opp().halfedge() == loop[(i+1)%loop.size()])
-                    loop[i] = InvalidHalfEdgeID;
+                {
-                    loop[(i+1)%loop.size()] = InvalidHalfEdgeID;
+                    VertexID vid = walker(loop[i]).vertex();
-                    ++did_work;
+                    if(vid != v)
-                    ++i;
+                        garbage_vertices.push_back(walker(loop[i]).vertex());
-                }
+                    garbage_halfedges.push_back(loop[i]);
-            for(size_t i=0;i<loop.size();++i)
+                    garbage_halfedges.push_back(loop[(i+1)%loop.size()]);
-                if(loop[i] != InvalidHalfEdgeID)
+                    loop[i] = InvalidHalfEdgeID;
-                    loop_tmp.push_back(loop[i]);
+                    loop[(i+1)%loop.size()] = InvalidHalfEdgeID;
-            loop = loop_tmp;
+                    ++did_work;
-        } while(did_work > 0 && loop.size() > 4);
+                    ++i;
+                }
-        // Check whether the loop is too long
+            for(size_t i=0;i<loop.size();++i)
-        float loop_len=0.0;
+                if(loop[i] != InvalidHalfEdgeID)
-        for(size_t i=0;i<loop.size();++i)
+                    loop_tmp.push_back(loop[i]);
-            loop_len += length(*this, loop[i]);
+            loop = loop_tmp;
-        if(loop_len > max_loop_length)
+        } while(did_work > 0 && loop.size() > 4);
-            return InvalidFaceID;
+        // Check whether the loop is too long
-        // The following block checks wheteher the same halfedge appears twice. In this
+        float loop_len=0.0;
-        // case we fail since it means that the one ring contains the same face twice.
+        for(size_t i=0;i<loop.size();++i)
-        vector<HalfEdgeID> loop_tmp = loop;
+            loop_len += length(*this, loop[i]);
-        sort(loop_tmp.begin(), loop_tmp.end());
+        if(loop_len > max_loop_length)
-        vector<HalfEdgeID>::iterator end_iter = unique(loop_tmp.begin(), loop_tmp.end());
+            return InvalidFaceID;
-        if(end_iter != loop_tmp.end())
-            return InvalidFaceID;
+        // The following block checks wheteher the same halfedge appears twice. In this
+        // case we fail since it means that the one ring contains the same face twice.
-        // Remove all faces and connected halfedges and the original vertex v.
+        vector<HalfEdgeID> loop_tmp = loop;
-        for(size_t i=0;i<garbage_vertices.size(); ++i)
+        sort(loop_tmp.begin(), loop_tmp.end());
-            kernel.remove_vertex(garbage_vertices[i]);
+        vector<HalfEdgeID>::iterator end_iter = unique(loop_tmp.begin(), loop_tmp.end());
-        for(size_t i=0;i<garbage_faces.size(); ++i)
+        if(end_iter != loop_tmp.end())
-            kernel.remove_face(garbage_faces[i]);
+            return InvalidFaceID;
-        for(size_t i=0;i<garbage_halfedges.size(); ++i)
-            kernel.remove_halfedge(garbage_halfedges[i]);
+        // Remove all faces and connected halfedges and the original vertex v.
-        if(!vertex_is_boundary)
+        for(size_t i=0;i<garbage_vertices.size(); ++i)
-            kernel.remove_vertex(v);
+            kernel.remove_vertex(garbage_vertices[i]);
+        for(size_t i=0;i<garbage_faces.size(); ++i)
-        // Create a new face for the merged one ring and link up all the halfedges
+            kernel.remove_face(garbage_faces[i]);
-        // in the loop
+        for(size_t i=0;i<garbage_halfedges.size(); ++i)
-        FaceID f = kernel.add_face();
+            kernel.remove_halfedge(garbage_halfedges[i]);
-        kernel.set_last(f,loop[0]);
+        if(!vertex_is_boundary)
-        for(size_t i=0;i<loop.size(); ++i)
+            kernel.remove_vertex(v);
-        {
-            kernel.set_face(loop[i], f);
+        // Create a new face for the merged one ring and link up all the halfedges
-            Walker hw = walker(loop[i]);
+        // in the loop
-            kernel.set_out(hw.vertex(),hw.opp().halfedge());
+        FaceID f = kernel.add_face();
-            link(loop[i],loop[(i+1)%loop.size()]);
+        kernel.set_last(f,loop[0]);
-            assert(hw.vertex() == walker(loop[(i+1)%loop.size()]).opp().vertex());
+        for(size_t i=0;i<loop.size(); ++i)
-        }
+        {
+            kernel.set_face(loop[i], f);
-        // Finally, we ensure boundary consitency for all vertices in the loop.
+            Walker hw = walker(loop[i]);
-        for(size_t i=0;i<loop.size(); ++i)
+            kernel.set_out(hw.vertex(),hw.opp().halfedge());
-        {
+            link(loop[i],loop[(i+1)%loop.size()]);
-            Walker hw = walker(loop[i]);
+            assert(hw.vertex() == walker(loop[(i+1)%loop.size()]).opp().vertex());
-            ensure_boundary_consistency(hw.vertex());
+        }
-        }
+        // Finally, we ensure boundary consitency for all vertices in the loop.
-        // Return the brand new merged face.
+        for(size_t i=0;i<loop.size(); ++i)
-        return f;
+        {
-    }
+            Walker hw = walker(loop[i]);
+            ensure_boundary_consistency(hw.vertex());
+        }
-    bool Manifold::merge_faces(FaceID f, HalfEdgeID h)
+        // Return the brand new merged face.
-    {
+        return f;
-        //assert that we're merging a valid face with the corresponding halfedge
+    }
-        assert(kernel.face(h) == f);
-        HalfEdgeID ho = kernel.opp(h);
-        FaceID fo = kernel.face(ho);
+    bool Manifold::merge_faces(FaceID f, HalfEdgeID h)
-        HalfEdgeID hn = kernel.next(h);
+    {
-        HalfEdgeID hon = kernel.next(ho);
+        //assert that we're merging a valid face with the corresponding halfedge
+        assert(kernel.face(h) == f);
-        if(fo == f)
-            return false;
+        HalfEdgeID ho = kernel.opp(h);
+        FaceID fo = kernel.face(ho);
-        //boundary case
+        HalfEdgeID hn = kernel.next(h);
-        if(kernel.vert(hn) == kernel.vert(hon))
+        HalfEdgeID hon = kernel.next(ho);
-            return false;
+        if(fo == f)
-        HalfEdgeID hp = kernel.prev(h);
+            return false;
-        HalfEdgeID hop = kernel.prev(ho);
-        VertexID v = kernel.vert(h);
+        //boundary case
-        VertexID vo = kernel.vert(ho);
+        if(kernel.vert(hn) == kernel.vert(hon))
+            return false;
-        if(valency(*this, v) < 3 || valency(*this, vo) < 3)
-            return false;
+        HalfEdgeID hp = kernel.prev(h);
+        HalfEdgeID hop = kernel.prev(ho);
-        link(hop, hn);
+        VertexID v = kernel.vert(h);
-        link(hp, hon);
+        VertexID vo = kernel.vert(ho);
-        kernel.set_out(vo, hon);
-        kernel.set_out(v, hn);
+        if(valency(*this, v) < 3 || valency(*this, vo) < 3)
-        kernel.set_last(f, hn);
+            return false;
-        HalfEdgeID hx = hon;
+        link(hop, hn);
+        link(hp, hon);
-        assert(kernel.face(hx) == fo);
+        kernel.set_out(vo, hon);
-        while(kernel.face(hx) != f){
+        kernel.set_out(v, hn);
-            kernel.set_face(hx, f);
+        kernel.set_last(f, hn);
-            hx = kernel.next(hx);
-        }
+        HalfEdgeID hx = hon;
-        ensure_boundary_consistency(v);
+        assert(kernel.face(hx) == fo);
-        ensure_boundary_consistency(vo);
+        while(kernel.face(hx) != f){
+            kernel.set_face(hx, f);
-        kernel.remove_halfedge(h);
+            hx = kernel.next(hx);
-        kernel.remove_halfedge(ho);
+        }
-        kernel.remove_face(fo);
+        ensure_boundary_consistency(v);
-        return true;
+        ensure_boundary_consistency(vo);
-    }
+        kernel.remove_halfedge(h);
-    FaceID Manifold::close_hole(HalfEdgeID h)
+        kernel.remove_halfedge(ho);
-    {
+        kernel.remove_face(fo);
-        // invalid face is a hole
-        if(kernel.face(h) == InvalidFaceID){
+        return true;
-            FaceID f = kernel.add_face();
+    }
-            kernel.set_last(f, h);
-            do{
+    FaceID Manifold::close_hole(HalfEdgeID h)
-                kernel.set_face(h, f);
+    {
-                h = kernel.next(h);
+        // invalid face is a hole
-            }
+        if(kernel.face(h) == InvalidFaceID){
-            while(kernel.face(h) != f);
+            FaceID f = kernel.add_face();
-            return f;
+            kernel.set_last(f, h);
-        }
+            do{
-        return kernel.face(h);
+                kernel.set_face(h, f);
-    }
+                h = kernel.next(h);
+            }
-    VertexID Manifold::slit_vertex(VertexID v, HalfEdgeID h_in, HalfEdgeID h_out)
+            while(kernel.face(h) != f);
-    {
+            return f;
-        assert(kernel.face(h_in) != InvalidFaceID);
+        }
-        assert(kernel.face(h_out) != InvalidFaceID);
+        return kernel.face(h);
-        assert(kernel.opp(h_out) != h_in);
+    }
-        VertexID v_new = kernel.add_vertex();
+    VertexID Manifold::slit_vertex(VertexID v, HalfEdgeID h_in, HalfEdgeID h_out)
-        pos(v_new) = pos(v);
+    {
-        HalfEdgeID h = kernel.prev(h_out);
+        assert(kernel.face(h_in) != InvalidFaceID);
-        kernel.set_vert(h, v_new);
+        assert(kernel.face(h_out) != InvalidFaceID);
-        while ( h != h_in) {
+        assert(kernel.opp(h_out) != h_in);
-            h = kernel.prev(kernel.opp(h));
-            kernel.set_vert(h, v_new);
+        VertexID v_new = kernel.add_vertex();
-        }
+        pos(v_new) = pos(v);
+        HalfEdgeID h = kernel.prev(h_out);
-        HalfEdgeID h_in_opp = kernel.opp(h_in);
+        kernel.set_vert(h, v_new);
-        HalfEdgeID hn_in, hn_in_opp;
+        while ( h != h_in) {
-        if(kernel.face(h_in_opp) != InvalidFaceID)
+            h = kernel.prev(kernel.opp(h));
-        {
+            kernel.set_vert(h, v_new);
-            hn_in = kernel.add_halfedge();
+        }
-            kernel.set_face(hn_in, InvalidFaceID);
-            glue(h_in_opp, hn_in);
+        HalfEdgeID h_in_opp = kernel.opp(h_in);
+        HalfEdgeID hn_in, hn_in_opp;
-            hn_in_opp = kernel.add_halfedge();
+        if(kernel.face(h_in_opp) != InvalidFaceID)
-            kernel.set_face(hn_in_opp, InvalidFaceID);
+        {
-            glue(h_in, hn_in_opp);
+            hn_in = kernel.add_halfedge();
+            kernel.set_face(hn_in, InvalidFaceID);
-            link(hn_in_opp, hn_in);
+            glue(h_in_opp, hn_in);
-            VertexID v_i = kernel.vert(h_in_opp);
+            hn_in_opp = kernel.add_halfedge();
-            kernel.set_vert(hn_in_opp, v_i);
+            kernel.set_face(hn_in_opp, InvalidFaceID);
-            kernel.set_out(v_i, hn_in);
+            glue(h_in, hn_in_opp);
-        }
-        else
+            link(hn_in_opp, hn_in);
-        {
-            hn_in_opp = h_in_opp;
+            VertexID v_i = kernel.vert(h_in_opp);
-            hn_in = kernel.prev(hn_in_opp);
+            kernel.set_vert(hn_in_opp, v_i);
-            h_in_opp = kernel.opp(hn_in);
+            kernel.set_out(v_i, hn_in);
         }
+        else
-        HalfEdgeID h_out_opp = kernel.opp(h_out);
+        {
-        HalfEdgeID hn_out,hn_out_opp;
+            hn_in_opp = h_in_opp;
-        if(kernel.face(h_out_opp) != InvalidFaceID)
+            hn_in = kernel.prev(hn_in_opp);
-        {
+            h_in_opp = kernel.opp(hn_in);
-            hn_out = kernel.add_halfedge();
+        }
-            kernel.set_face(hn_out, InvalidFaceID);
-            glue(h_out_opp, hn_out);
+        HalfEdgeID h_out_opp = kernel.opp(h_out);
+        HalfEdgeID hn_out,hn_out_opp;
-            hn_out_opp = kernel.add_halfedge();
+        if(kernel.face(h_out_opp) != InvalidFaceID)
-            kernel.set_face(hn_out_opp, InvalidFaceID);
+        {
-            glue(h_out, hn_out_opp);
+            hn_out = kernel.add_halfedge();
+            kernel.set_face(hn_out, InvalidFaceID);
-            link(hn_out, hn_out_opp);
+            glue(h_out_opp, hn_out);
-            VertexID v_o = kernel.vert(h_out);
+            hn_out_opp = kernel.add_halfedge();
-            kernel.set_vert(hn_out, v_o);
+            kernel.set_face(hn_out_opp, InvalidFaceID);
-            kernel.set_out(v_o, hn_out_opp);
+            glue(h_out, hn_out_opp);
-        }
-        else
+            link(hn_out, hn_out_opp);
-        {
-            hn_out_opp = h_out_opp;
+            VertexID v_o = kernel.vert(h_out);
-            hn_out = kernel.next(hn_out_opp);
+            kernel.set_vert(hn_out, v_o);
-            h_out_opp = kernel.opp(hn_out);
+            kernel.set_out(v_o, hn_out_opp);
         }
+        else
-        link(hn_out_opp, hn_in_opp);
+        {
-        link(hn_in, hn_out);
+            hn_out_opp = h_out_opp;
+            hn_out = kernel.next(hn_out_opp);
-        kernel.set_vert(hn_in, v);
+            h_out_opp = kernel.opp(hn_out);
-        kernel.set_vert(hn_out_opp, v_new);
+        }
-        kernel.set_out(v, hn_out);
+        link(hn_out_opp, hn_in_opp);
-        kernel.set_out(v_new, hn_in_opp);
+        link(hn_in, hn_out);
-        return v_new;
+        kernel.set_vert(hn_in, v);
-    }
+        kernel.set_vert(hn_out_opp, v_new);
+        kernel.set_out(v, hn_out);
-    HalfEdgeID Manifold::slit_edges(VertexAttributeVector<int>& insel)
+        kernel.set_out(v_new, hn_in_opp);
-    {
-        HalfEdgeID h;
+        return v_new;
-        for(auto vid : vertices())
+    }
-        {
-            if(insel[vid])
-            {
+    HalfEdgeID Manifold::slit_edges(VertexAttributeVector<int>& insel)
-                HalfEdgeID h_in = InvalidHalfEdgeID, h_out = InvalidHalfEdgeID;
+    {
-                Walker w = walker(vid);
+        HalfEdgeID h;
-                while(!w.full_circle())
+        for(auto vid : vertices())
-                {
+        {
-                    if(insel[w.vertex()]) {
+            if(insel[vid])
-                        if(h_in == InvalidHalfEdgeID) {
+            {
-                            if(w.opp().face() == InvalidFaceID)
+                HalfEdgeID h_in = InvalidHalfEdgeID, h_out = InvalidHalfEdgeID;
-                                h_in = w.opp().next().opp().halfedge();
+                Walker w = walker(vid);
-                            else
+                while(!w.full_circle())
-                                h_in = w.opp().halfedge();
+                {
-                        }
+                    if(insel[w.vertex()]) {
-                        else {
+                        if(h_in == InvalidHalfEdgeID) {
-                            if(w.face() == InvalidFaceID)
+                            if(w.opp().face() == InvalidFaceID)
-                                h_out = w.prev().opp().halfedge();
+                                h_in = w.opp().next().opp().halfedge();
                             else
-                                h_out = w.halfedge();
+                                h_in = w.opp().halfedge();
-                            break;
+                        }
-                        }
+                        else {
-                    }
+                            if(w.face() == InvalidFaceID)
-                    w = w.circulate_vertex_ccw();
+                                h_out = w.prev().opp().halfedge();
-                }
+                            else
-                if(h_in != InvalidHalfEdgeID &&
+                                h_out = w.halfedge();
-                   h_out != InvalidHalfEdgeID) {
+                            break;
-                    VertexID v_new = slit_vertex(vid, h_in, h_out);
+                        }
-                    h = walker(v_new).halfedge();
+                    }
-                }
+                    w = w.circulate_vertex_ccw();
-            }
+                }
-        }
+                if(h_in != InvalidHalfEdgeID &&
-        return h;
+                   h_out != InvalidHalfEdgeID) {
-    }
+                    VertexID v_new = slit_vertex(vid, h_in, h_out);
+                    h = walker(v_new).halfedge();
+                }
-    void Manifold::flip_edge(HalfEdgeID h)
+            }
-    {
+        }
-        HalfEdgeID hn = kernel.next(h);
+        return h;
-        HalfEdgeID hp = kernel.prev(h);
+    }
-        HalfEdgeID ho = kernel.opp(h);
-        HalfEdgeID hon = kernel.next(ho);
-        HalfEdgeID hop = kernel.prev(ho);
+    void Manifold::flip_edge(HalfEdgeID h)
+    {
-        FaceID hf = kernel.face(h);
+        HalfEdgeID hn = kernel.next(h);
-        FaceID hof = kernel.face(ho);
+        HalfEdgeID hp = kernel.prev(h);
+        HalfEdgeID ho = kernel.opp(h);
-        VertexID hv = kernel.vert(h);
+        HalfEdgeID hon = kernel.next(ho);
-        VertexID hnv = kernel.vert(hn);
+        HalfEdgeID hop = kernel.prev(ho);
-        VertexID hov = kernel.vert(ho);
-        VertexID honv = kernel.vert(hon);
+        FaceID hf = kernel.face(h);
+        FaceID hof = kernel.face(ho);
-        // update face connectivity of halfedges changing face
-        kernel.set_face(hop, hf);
+        VertexID hv = kernel.vert(h);
-        kernel.set_face(hp, hof);
+        VertexID hnv = kernel.vert(hn);
+        VertexID hov = kernel.vert(ho);
-        // connectivity of faces with halfedges of flipped edge
+        VertexID honv = kernel.vert(hon);
-        kernel.set_last(hf, h);
-        kernel.set_last(hof, ho);
+        // update face connectivity of halfedges changing face
+        kernel.set_face(hop, hf);
-        // connectivity of halfedges of first face after flip
+        kernel.set_face(hp, hof);
-        link(hn, h);
-        link(h, hop);
+        // connectivity of faces with halfedges of flipped edge
-        link(hop, hn);
+        kernel.set_last(hf, h);
+        kernel.set_last(hof, ho);
-        // connectivity of halfedges of second face after flip
-        link(hon, ho);
+        // connectivity of halfedges of first face after flip
-        link(ho, hp);
+        link(hn, h);
-        link(hp, hon);
+        link(h, hop);
+        link(hop, hn);
-        // connectivity of flipped edge and corresponding vertices
-        kernel.set_vert(h, honv);
+        // connectivity of halfedges of second face after flip
-        kernel.set_vert(ho, hnv);
+        link(hon, ho);
+        link(ho, hp);
-        if(kernel.out(hv) == ho)
+        link(hp, hon);
-            kernel.set_out(hv, hn);
-        if(kernel.out(hov) == h)
+        // connectivity of flipped edge and corresponding vertices
-            kernel.set_out(hov, hon);
+        kernel.set_vert(h, honv);
+        kernel.set_vert(ho, hnv);
-        //
-        //        // if the flip occurs next to a boundary, ensure the boundary consistency
+        if(kernel.out(hv) == ho)
-        //        ensure_boundary_consistency(hv);
+            kernel.set_out(hv, hn);
-        //        ensure_boundary_consistency(hov);
+        if(kernel.out(hov) == h)
-    }
+            kernel.set_out(hov, hon);
+        //
-    /**********************************************
+        //        // if the flip occurs next to a boundary, ensure the boundary consistency
-     * Private functions
+        //        ensure_boundary_consistency(hv);
-     **********************************************/
+        //        ensure_boundary_consistency(hov);
+    }
-    template<typename size_type, typename float_type, typename int_type>
-    void Manifold::build_template(  size_type no_vertices,
-                                  const float_type* vertvec,
+    /**********************************************
-                                  size_type no_faces,
+     * Private functions
-                                  const int_type* facevec,
+     **********************************************/
-                                  const int_type* indices)
-    {
+    template<typename size_type, typename float_type, typename int_type>
-        vector<VertexID> vids(no_vertices);
+    void Manifold::build_template(  size_type no_vertices,
+                                  const float_type* vertvec,
-        // create vertices corresponding to positions stored in vertvec
+                                  size_type no_faces,
-        for(size_t i=0;i<no_vertices;++i)
+                                  const int_type* facevec,
-        {
+                                  const int_type* indices)
-            const float_type* v0 = &vertvec[i*3];
+    {
-            pos(vids[i] = kernel.add_vertex()) = Manifold::Vec(v0[0], v0[1], v0[2]);
+        vector<VertexID> vids(no_vertices);
-        }
+        // create vertices corresponding to positions stored in vertvec
-        //map over the created edges - needed to preserve edge uniqueness
+        for(size_t i=0;i<no_vertices;++i)
-        typedef map<EdgeKey, Edge> EdgeMap;
+        {
-        EdgeMap edge_map;
+            const float_type* v0 = &vertvec[i*3];
+            pos(vids[i] = kernel.add_vertex()) = Manifold::Vec(v0[0], v0[1], v0[2]);
-        // counter that jumps between faces in indices vector
+        }
-        int_type n  = 0;
+        auto hash_fun = [](const EdgeKey& k) {return k.hash();};
-        // create faces correspponding to faces stored in facevec
+        //map over the created edges - needed to preserve edge uniqueness
-        for(size_type i = 0; i < no_faces; ++i){
+        typedef unordered_map<EdgeKey, Edge, decltype(hash_fun)> EdgeMap;
-            //amount of vertices in current face
+        EdgeMap edge_map(no_vertices+no_faces,hash_fun);
-            size_type N = facevec[i];
-            vector<HalfEdgeID> fh;
+        // counter that jumps between faces in indices vector
+        int_type n  = 0;
-            //each face indice results corresponds to 1 edge, 2 halfedges
-            for(size_type j = 0; j < N; ++j){
+        // create faces correspponding to faces stored in facevec
-                // ensure indice integrity
+        for(size_type i = 0; i < no_faces; ++i){
+            //amount of vertices in current face
-                assert(static_cast<size_type>(indices[j + n]) < no_vertices);
+            size_type N = facevec[i];
-                assert(static_cast<size_type>(indices[(j + 1) % N + n]) < no_vertices);
+            vector<HalfEdgeID> fh;
+            //each face indice results corresponds to 1 edge, 2 halfedges
-                // each iteration uses two indices from the face
+            for(size_type j = 0; j < N; ++j){
-                VertexID v0 = vids[static_cast<size_type>(indices[j + n])];
+                // ensure indice integrity
-                VertexID v1 = vids[static_cast<size_type>(indices[(j + 1) % N + n])];
+                assert(static_cast<size_type>(indices[j + n]) < no_vertices);
-                // create key and search map for edge
+                assert(static_cast<size_type>(indices[(j + 1) % N + n]) < no_vertices);
-                EdgeKey ek(v0, v1);
-                EdgeMap::iterator em_iter = edge_map.find(ek);
+                // each iteration uses two indices from the face
-                // current edge has not been created
+                VertexID v0 = vids[static_cast<size_type>(indices[j + n])];
-                if(em_iter == edge_map.end()){
+                VertexID v1 = vids[static_cast<size_type>(indices[(j + 1) % N + n])];
-                    // create edge for map
-                    Edge e;
+                // create key and search map for edge
-                    e.h0 = kernel.add_halfedge();
+                EdgeKey ek(v0, v1);
-                    e.h1 = kernel.add_halfedge();
+                typename EdgeMap::iterator em_iter = edge_map.find(ek);
-                    e.count = 1;
+                // current edge has not been created
-                    // glue operation: 1 edge = 2 glued halfedges
+                if(em_iter == edge_map.end()){
-                    glue(e.h0, e.h1);
+                    // create edge for map
+                    Edge e;
-                    // update vertices with their outgoing halfedges
+                    e.h0 = kernel.add_halfedge();
-                    kernel.set_out(v0, e.h0);
+                    e.h1 = kernel.add_halfedge();
-                    kernel.set_out(v1, e.h1);
+                    e.count = 1;
-                    // update halfedges with the vertices they point to
+                    // glue operation: 1 edge = 2 glued halfedges
-                    kernel.set_vert(e.h0, v1);
+                    glue(e.h0, e.h1);
-                    kernel.set_vert(e.h1, v0);
+                    // update vertices with their outgoing halfedges
-                    // update map
+                    kernel.set_out(v0, e.h0);
-                    edge_map[ek] = e;
+                    kernel.set_out(v1, e.h1);
-                    // update vector of halfedges belonging to current face
+                    // update halfedges with the vertices they point to
-                    fh.push_back(e.h0);
+                    kernel.set_vert(e.h0, v1);
-                }
+                    kernel.set_vert(e.h1, v0);
-                else{
-                    // update current face with halfedge from edge
+                    // update map
-                    fh.push_back(em_iter->second.h1);
+                    edge_map[ek] = e;
-                    // asserting that a halfedge is visited exactly twice;
+                    // update vector of halfedges belonging to current face
-                    // once for each face on either side of the edge.
+                    fh.push_back(e.h0);
-                    em_iter->second.count++;
+                }
-                    assert(em_iter->second.count == 2);
+                else{
-                }
+                    // update current face with halfedge from edge
-            }
+                    fh.push_back(em_iter->second.h1);
-            FaceID fid = kernel.add_face();
+                    // asserting that a halfedge is visited exactly twice;
-            for(size_type j = 0; j < N; ++j){
+                    // once for each face on either side of the edge.
-                // update halfedge with face
+                    em_iter->second.count++;
-                kernel.set_face(fh[j], fid);
+                    assert(em_iter->second.count == 2);
+                }
-                // link operation: link two halfedges in the same face
+            }
-                link(fh[j], fh[(j + 1) % N]);
-            }
+            FaceID fid = kernel.add_face();
-            //update face with the first halfedge created
+            for(size_type j = 0; j < N; ++j){
-            kernel.set_last(fid, fh[0]);
+                // update halfedge with face
+                kernel.set_face(fh[j], fid);
-            // step to first index of next face
-            n += N;
+                // link operation: link two halfedges in the same face
-        }
+                link(fh[j], fh[(j + 1) % N]);
+            }
-        // test for unused vertices
+            //update face with the first halfedge created
-        for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
+            kernel.set_last(fid, fh[0]);
-            assert( (*v) != InvalidVertexID);
-            if(kernel.out(*v) == InvalidHalfEdgeID)
+            // step to first index of next face
-                kernel.remove_vertex(*v);
+            n += N;
         }
-        // boundary check while avoiding unused vertices
+        // test for unused vertices
         for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
-            if((*v) != InvalidVertexID)
+            assert( (*v) != InvalidVertexID);
-                ensure_boundary_consistency(*v);
+            if(kernel.out(*v) == InvalidHalfEdgeID)
-        }
+                kernel.remove_vertex(*v);
-    }
+        }
-    void Manifold::link(HalfEdgeID h0, HalfEdgeID h1)
-    {
+        // boundary check while avoiding unused vertices
-        kernel.set_next(h0, h1);
+        for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
-        kernel.set_prev(h1, h0);
+            if((*v) != InvalidVertexID && in_use(*v))
-    }
+                ensure_boundary_consistency(*v);
-    void Manifold::glue(HalfEdgeID h0, HalfEdgeID h1)
+        }
-    {
+    }
-        kernel.set_opp(h0, h1);
+    void Manifold::link(HalfEdgeID h0, HalfEdgeID h1)
-        kernel.set_opp(h1, h0);
+    {
-    }
+        kernel.set_next(h0, h1);
-    void Manifold::ensure_boundary_consistency(VertexID v)
+        kernel.set_prev(h1, h0);
-    {
+    }
-        // boundary consistency check by performing two vertex circulations
+    void Manifold::glue(HalfEdgeID h0, HalfEdgeID h1)
-        HalfEdgeID h = kernel.out(v);
+    {
-        HalfEdgeID last = h;
+        kernel.set_opp(h0, h1);
+        kernel.set_opp(h1, h0);
-        // step 1: circle through edges pointing away from vertex until reaching a null face
+    }
-        while(kernel.face(h) != InvalidFaceID){
+    void Manifold::ensure_boundary_consistency(VertexID v)
-            h = kernel.opp(kernel.prev(h));
+    {
-            if(h == last) // We came full circle - vertex not boundary - return.
+        // boundary consistency check by performing two vertex circulations
-                return;
+        HalfEdgeID h = kernel.out(v);
-        }
+        HalfEdgeID last = h;
-        // null face encountered, we update our vertex with half edge index and prepare for step 2
-        kernel.set_out(v, h);
+        int c = 0;
-        HalfEdgeID ho = kernel.opp(h);
+        // step 1: circle through edges pointing away from vertex until reaching a null face
+        while(kernel.face(h) != InvalidFaceID){
-        // step 2: circle through edges pointing towards vertex until reaching a null face
+            h = kernel.opp(kernel.prev(h));
-        while(kernel.face(ho) != InvalidFaceID){
+            if(h == last || ++c == 1e6) // We came full circle - vertex not boundary - return.
-            ho = kernel.opp(kernel.next(ho));
+                return;
         }
-        // null face encountered again, we update our edge with vertex index
+        // null face encountered, we update our vertex with half edge index and prepare for step 2
-        kernel.set_vert(ho, v);
+        kernel.set_out(v, h);
+        HalfEdgeID ho = kernel.opp(h);
-        // remaining step is to make the in and out going edges link to each other
-        link(ho, h);
+        // step 2: circle through edges pointing towards vertex until reaching a null face
-    }
+        while(kernel.face(ho) != InvalidFaceID){
-    void Manifold::remove_face_if_degenerate(HalfEdgeID h)
+            ho = kernel.opp(kernel.next(ho));
-    {
+        }
-        // face is degenerate if there is only two halfedges in face loop
+        // null face encountered again, we update our edge with vertex index
-        if(kernel.next(kernel.next(h)) == h)
+        kernel.set_vert(ho, v);
-        {
-            HalfEdgeID hn = kernel.next(h);
+        // remaining step is to make the in and out going edges link to each other
-            HalfEdgeID ho = kernel.opp(h);
+        link(ho, h);
-            HalfEdgeID hno = kernel.opp(hn);
+    }
-            VertexID hv = kernel.vert(h);
+    void Manifold::remove_face_if_degenerate(HalfEdgeID h)
-            VertexID hnv = kernel.vert(hn);
+    {
-            FaceID f = kernel.face(h);
+        // face is degenerate if there is only two halfedges in face loop
+        if(kernel.next(kernel.next(h)) == h)
-            assert(ho != hn);
+        {
-            assert(h != hno);
+            HalfEdgeID hn = kernel.next(h);
-            assert(hv != hnv);
+            HalfEdgeID ho = kernel.opp(h);
+            HalfEdgeID hno = kernel.opp(hn);
-            // glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
+            VertexID hv = kernel.vert(h);
-            glue(ho, hno);
+            VertexID hnv = kernel.vert(hn);
+            FaceID f = kernel.face(h);
-            // make v and vn point to opposite and next opposite halfedge, obsoleting h and hn from loop
-            kernel.set_out(hnv, hno);
+            assert(ho != hn);
-            kernel.set_out(hv, ho);
+            assert(h != hno);
+            assert(hv != hnv);
-            // if face owning h is valid, remove face
-            if(f != InvalidFaceID)
+            // glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
-                kernel.remove_face(f);
+            glue(ho, hno);
-            // remove the two invalid halfedges and the invalid face loop
-            kernel.remove_halfedge(h);
+            // make v and vn point to opposite and next opposite halfedge, obsoleting h and hn from loop
-            kernel.remove_halfedge(hn);
+            kernel.set_out(hnv, hno);
+            kernel.set_out(hv, ho);
-            // verify that v and vn is sane after removing the degenerate face
-            ensure_boundary_consistency(hv);
+            // if face owning h is valid, remove face
-            ensure_boundary_consistency(hnv);
+            if(f != InvalidFaceID)
-        }
+                kernel.remove_face(f);
-    }
+            // remove the two invalid halfedges and the invalid face loop
+            kernel.remove_halfedge(h);
-    vector<HalfEdgeID> Manifold::bridge_faces(FaceID f0, FaceID f1, const vector<pair<VertexID, VertexID> >& pairs)
+            kernel.remove_halfedge(hn);
-    {
-        // Let N be the number of vertex pairs to connect by edges
+            // verify that v and vn is sane after removing the degenerate face
-        size_t N = pairs.size();
+            ensure_boundary_consistency(hv);
+            ensure_boundary_consistency(hnv);
-        // We now create N pairs of edges, glue each pair and let them point to the appropriate
+        }
-        // vertices.
+    }
-        vector<HalfEdgeID> new_halfedges(N);
-        vector<HalfEdgeID> new_halfedges_opp(N);
+    vector<HalfEdgeID> Manifold::bridge_faces(FaceID f0, FaceID f1, const vector<pair<VertexID, VertexID> >& pairs)
-        for(size_t i=0;i<N; ++i)
+    {
-        {
+        // Let N be the number of vertex pairs to connect by edges
-            new_halfedges[i] = kernel.add_halfedge();
+        size_t N = pairs.size();
-            new_halfedges_opp[i] = kernel.add_halfedge();
-            glue(new_halfedges[i], new_halfedges_opp[i]);
+        // We now create N pairs of edges, glue each pair and let them point to the appropriate
-            kernel.set_vert(new_halfedges[i], pairs[i].second);
+        // vertices.
-            kernel.set_vert(new_halfedges_opp[i], pairs[i].first);
+        vector<HalfEdgeID> new_halfedges(N);
-        }
+        vector<HalfEdgeID> new_halfedges_opp(N);
+        for(size_t i=0;i<N; ++i)
-        // We now maintain some halfedge indices to right before
+        {
-        // and right after the point we are trying to connect on
+            new_halfedges[i] = kernel.add_halfedge();
-        // each of the two loops.
+            new_halfedges_opp[i] = kernel.add_halfedge();
-        HalfEdgeID h0p = kernel.last(f0);
+            glue(new_halfedges[i], new_halfedges_opp[i]);
-        HalfEdgeID h1n = kernel.last(f1);
+            kernel.set_vert(new_halfedges[i], pairs[i].second);
+            kernel.set_vert(new_halfedges_opp[i], pairs[i].first);
-        // loop over all connections and link the new halfedges with the old
+        }
-        // ones.
-        for(size_t i=0;i<N; ++i)
+        // We now maintain some halfedge indices to right before
-        {
+        // and right after the point we are trying to connect on
-            while(kernel.vert(h0p) != pairs[i].first)
+        // each of the two loops.
-                h0p = kernel.next(h0p);
+        HalfEdgeID h0p = kernel.last(f0);
-            while(kernel.vert(kernel.prev(h1n)) != pairs[i].second)
+        HalfEdgeID h1n = kernel.last(f1);
-                h1n = kernel.prev(h1n);
+        // loop over all connections and link the new halfedges with the old
-            HalfEdgeID h0n = kernel.next(h0p);
+        // ones.
-            HalfEdgeID h1p = kernel.prev(h1n);
+        for(size_t i=0;i<N; ++i)
+        {
-            link(h0p, new_halfedges[i]);
+            while(kernel.vert(h0p) != pairs[i].first)
-            link(new_halfedges[i],h1n);
+                h0p = kernel.next(h0p);
-            link(h1p, new_halfedges_opp[i]);
+            while(kernel.vert(kernel.prev(h1n)) != pairs[i].second)
-            link(new_halfedges_opp[i],h0n);
+                h1n = kernel.prev(h1n);
-            h0p = new_halfedges_opp[i];
+            HalfEdgeID h0n = kernel.next(h0p);
-            h1n = new_halfedges_opp[i];
+            HalfEdgeID h1p = kernel.prev(h1n);
-        }
+            link(h0p, new_halfedges[i]);
-        // Create the faces and their connections
+            link(new_halfedges[i],h1n);
-        for(size_t i=0;i<N; ++i)
+            link(h1p, new_halfedges_opp[i]);
-        {
+            link(new_halfedges_opp[i],h0n);
-            HalfEdgeID last = new_halfedges[i];
-            FaceID f = kernel.add_face();
+            h0p = new_halfedges_opp[i];
-            kernel.set_last(f, last);
+            h1n = new_halfedges_opp[i];
-            HalfEdgeID h = last;
+        }
-            do
-            {
+        // Create the faces and their connections
-                kernel.set_face(h, f);
+        for(size_t i=0;i<N; ++i)
-                h = kernel.next(h);
+        {
-            } while(h != last);
+            HalfEdgeID last = new_halfedges[i];
-        }
+            FaceID f = kernel.add_face();
+            kernel.set_last(f, last);
-        // Delete the old faces that were bridged.
+            HalfEdgeID h = last;
-        kernel.remove_face(f0);
+            do
-        kernel.remove_face(f1);
+            {
+                kernel.set_face(h, f);
-        new_halfedges.insert(new_halfedges.end(), new_halfedges_opp.begin(), new_halfedges_opp.end());
+                h = kernel.next(h);
-        return new_halfedges;
+            } while(h != last);
-    }
+        }
-    /***************************************************
+        // Delete the old faces that were bridged.
-     * Namespace functions
+        kernel.remove_face(f0);
-     ***************************************************/
+        kernel.remove_face(f1);
-    bool valid(const Manifold& m)
-    {
+        new_halfedges.insert(new_halfedges.end(), new_halfedges_opp.begin(), new_halfedges_opp.end());
-        // Verify components of halfedges
+        return new_halfedges;
-        for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
+    }
-            Walker j = m.walker(*h);
+    /***************************************************
-            if(j.vertex() == InvalidVertexID){
+     * Namespace functions
-                cout << "Halfedge lacks vert" << endl;
+     ***************************************************/
-                return false;
+    bool valid(const Manifold& m)
-            }
+    {
-            if(j.next().halfedge() == InvalidHalfEdgeID){
+        // Verify components of halfedges
-                cout << "Halfedge lacks next" << endl;
+        for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
-                return false;
+            Walker j = m.walker(*h);
-            }
-            if(j.prev().halfedge() == InvalidHalfEdgeID){
+            if(j.vertex() == InvalidVertexID){
-                cout << "Halfedge lacks prev" << endl;
+                cout << "Halfedge lacks vert" << endl;
                 return false;
             }
-            if(j.opp().halfedge() == InvalidHalfEdgeID){
+            if(j.next().halfedge() == InvalidHalfEdgeID){
-                cout << "Halfedge lacks opp" << endl;
+                cout << "Halfedge lacks next" << endl;
                 return false;
             }
+            if(j.prev().halfedge() == InvalidHalfEdgeID){
-        }
+                cout << "Halfedge lacks prev" << endl;
-        // Verify components of vertices
+                return false;
-        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
+            }
-            vector<VertexID> link;
+            if(j.opp().halfedge() == InvalidHalfEdgeID){
+                cout << "Halfedge lacks opp" << endl;
-            // circulate the halfedges of vertex
+                return false;
-            for(Walker j = m.walker(*v); !j.full_circle(); j = j.circulate_vertex_cw()){
+            }
-                // test halfedges around v
-                if(j.halfedge() == InvalidHalfEdgeID){
+        }
-                    cout << "Vertex circulation produced invalid halfedge" << endl;
+        // Verify components of vertices
-                    return false;
+        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
-                }
+            vector<VertexID> link;
-                VertexID ring_v = j.vertex();
+            // circulate the halfedges of vertex
-                // test one-ring for multiple occurences of vertex
+            for(Walker j = m.walker(*v); !j.full_circle(); j = j.circulate_vertex_cw()){
-                if(find(link.begin(), link.end(), ring_v) != link.end()){
+                // test halfedges around v
-                    cout << "Vertex appears two times in one-ring of vertex" << endl;
+                if(j.halfedge() == InvalidHalfEdgeID){
-                    return false;
+                    cout << "Vertex circulation produced invalid halfedge" << endl;
-                }
+                    return false;
-                link.push_back(ring_v);
+                }
+                VertexID ring_v = j.vertex();
-                // test for infinite loop around vertex
-                if(static_cast<size_t>(j.no_steps()) > m.no_vertices()){
+                // test one-ring for multiple occurences of vertex
-                    cout << "Vertex loop contains more vertices than manifold" << endl;
+                if(find(link.begin(), link.end(), ring_v) != link.end()){
-                    return false;
+                    cout << "Vertex appears two times in one-ring of vertex" << endl;
-                }
+                    return false;
-            }
+                }
+                link.push_back(ring_v);
-            // test one_ring size for boundary consistency
-            if(link.size() <= 2){
+                // test for infinite loop around vertex
-                Walker j = m.walker(*v);
+                if(static_cast<size_t>(j.no_steps()) > m.no_vertices()){
+                    cout << "Vertex loop CW contains more vertices than manifold" << endl;
-                if(j.face() != InvalidFaceID && j.opp().face() != InvalidFaceID)
+                    return false;
-                {
+                }
-                    if(link.size()==1)
+            }
-                        cout << "Vertex contains only a single incident edge" << endl;
-                    //
+            for(Walker j = m.walker(*v); !j.full_circle(); j = j.circulate_vertex_ccw()) {
-                    //                    cout << "Vertex contains only " << link.size() <<" incident edges" << endl;
+                if(static_cast<size_t>(j.no_steps()) > m.no_vertices()){
-                }
+                    cout << "Vertex loop CCW contains more vertices than manifold" << endl;
-                else
+                    return false;
-                    cout << "Boundary vertex contains only " << link.size() <<" incident edges" << endl;
+                }
             }
-        }
-        // verify components of faces
+            // test one_ring size for boundary consistency
-        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
+            if(link.size() <= 2){
-            // count edges on face
+                Walker j = m.walker(*v);
-            Walker j = m.walker(*f);
+                if(j.face() != InvalidFaceID && j.opp().face() != InvalidFaceID)
-            for(; !j.full_circle(); j = j.circulate_face_cw()){
+                {
-                // test that all halfedges in faces bind properly to their face
+                    if(link.size()==1)
-                if(j.face() != *f){
+                        cout << "Vertex contains only a single incident edge" << endl;
-                    cout << "Face is inconsistent, halfedge is not bound to face" << endl;
+                    //
-                    return false;
+                    //                    cout << "Vertex contains only " << link.size() <<" incident edges" << endl;
                 }
-            }
+                else
-            // test faces for valid geometrical properties
+                    cout << "Boundary vertex contains only " << link.size() <<" incident edges" << endl;
-            if(j.no_steps() < 3){
+            }
-                cout << "Face contains less than 3 edges" << endl;
+        }
-                return false;
+        // verify components of faces
-            }
+        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
-            // test for infinite loop around face
+            // count edges on face
-            if(j.no_steps() > m.no_halfedges() * 0.5f){
+            Walker j = m.walker(*f);
-                cout << "Face loop contains more halfedges than manifold" << endl;
-                return false;
+            for(; !j.full_circle(); j = j.circulate_face_cw()){
-            }
+                // test that all halfedges in faces bind properly to their face
-        }
+                if(j.face() != *f){
-        return true;
+                    cout << "Face is inconsistent, halfedge is not bound to face" << endl;
-    }
+                    return false;
+                }
-    void bbox(const Manifold& m, Manifold::Vec& pmin, Manifold::Vec& pmax)
+            }
-    {
+            // test faces for valid geometrical properties
-        if(m.no_vertices()==0)
+            if(j.no_steps() < 3){
-            return;
+                cout << "Face contains less than 3 edges" << endl;
-        VertexIDIterator v = m.vertices_begin();
+                return false;
-        pmin = pmax = m.pos(*v);
+            }
-        ++v;
+            // test for infinite loop around face
-        for(; v != m.vertices_end(); ++v){
+            if(j.no_steps() > m.no_halfedges() * 0.5f){
-            pmin = v_min(m.pos(*v), pmin);
+                cout << "Face loop contains more halfedges than manifold" << endl;
-            pmax = v_max(m.pos(*v), pmax);
+                return false;
-        }
+            }
-    }
+        }
+        return true;
-    void bsphere(const Manifold& m, Manifold::Vec& c, float& r)
+    }
-    {
-        Manifold::Vec p0, p7;
+    void bbox(const Manifold& m, Manifold::Vec& pmin, Manifold::Vec& pmax)
-        bbox(m, p0, p7);
+    {
-        Manifold::Vec rad = (p7 - p0) * 0.5f;
+        if(m.no_vertices()==0)
-        c = p0 + rad;
+            return;
-        r = rad.length();
+        VertexIDIterator v = m.vertices_begin();
-    }
+        pmin = pmax = m.pos(*v);
+        ++v;
+        for(; v != m.vertices_end(); ++v){
+            pmin = v_min(m.pos(*v), pmin);
-    bool precond_collapse_edge(const Manifold& m, HalfEdgeID h)
+            pmax = v_max(m.pos(*v), pmax);
-    {
+        }
-        Walker hew = m.walker(h);
+    }
-        HalfEdgeID ho = hew.opp().halfedge();
-        VertexID v0 = hew.opp().vertex();
+    void bsphere(const Manifold& m, Manifold::Vec& c, float& r)
-        VertexID v1 = hew.vertex();
+    {
+        Manifold::Vec p0, p7;
-        // get the one-ring of v0
+        bbox(m, p0, p7);
-        vector<VertexID> link0;
+        Manifold::Vec rad = (p7 - p0) * 0.5f;
-        vector<FaceID> faces0;
+        c = p0 + rad;
-        for(Walker vj = m.walker(h);
+        r = rad.length();
-            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
+    }
-            link0.push_back(vj.vertex());
-            if(vj.halfedge() != h)
-                faces0.push_back(vj.face());
-        }
+    bool precond_collapse_edge(const Manifold& m, HalfEdgeID h)
-        assert(link0.size() > 1);
+    {
+        Walker hew = m.walker(h);
-        // get the one-ring of v1
+        HalfEdgeID ho = hew.opp().halfedge();
-        vector<VertexID> link1;
+        VertexID v0 = hew.opp().vertex();
-        vector<FaceID> faces1;
+        VertexID v1 = hew.vertex();
-        for(Walker vj = m.walker(ho);
-            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
+        // get the one-ring of v0
-            link1.push_back(vj.vertex());
+        vector<VertexID> link0;
-            if(vj.halfedge() != ho)
+        vector<FaceID> faces0;
-                faces1.push_back(vj.face());
+        int k = 0;
-        }
+        for(Walker vj = m.walker(h);
-        assert(link1.size() > 1);
+            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
+            link0.push_back(vj.vertex());
-        // sort the vertices of the two rings
+            if(vj.halfedge() != h)
-        sort(link0.begin(), link0.end());
+                faces0.push_back(vj.face());
-        sort(link1.begin(), link1.end());
+            if(++k>m.no_vertices())
+            {
-        // get the intersection of the shared vertices from both rings
+                cout << "mesh is corrupted" << endl;
-        vector<VertexID> lisect;
+                return false;
-        std::back_insert_iterator<vector<VertexID> > lii(lisect);
+            }
-        set_intersection(link0.begin(), link0.end(),
+        }
-                         link1.begin(), link1.end(),
+        assert(link0.size() > 1);
-                         lii);
+        // get the one-ring of v1
-        // sort the faces of the two rings
+        vector<VertexID> link1;
-        sort(faces0.begin(), faces0.end());
+        vector<FaceID> faces1;
-        sort(faces1.begin(), faces1.end());
+        k=0;
+        for(Walker vj = m.walker(ho);
-        // get the intersection of the shared faces from both rings
+            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
-        vector<FaceID> fisect;
+            link1.push_back(vj.vertex());
-        std::back_insert_iterator<vector<FaceID> > fii(fisect);
+            if(vj.halfedge() != ho)
-        set_intersection(faces0.begin(), faces0.end(),
+                faces1.push_back(vj.face());
-                         faces1.begin(), faces1.end(),
+            if(++k>m.no_vertices())
-                         fii);
+            {
-        if(fisect.size() > 0)
+                cout << "mesh is corrupted" << endl;
-            return false;
+                return false;
+            }
-        int k = 0;
+        }
-        // if the adjacent face is a triangle (see 2)
+        assert(link1.size() > 1);
-        if(hew.next().next().next().halfedge() == h){
-            VertexID v = hew.next().vertex();
+        // sort the vertices of the two rings
+        sort(link0.begin(), link0.end());
-            // valency test (see 5)
+        sort(link1.begin(), link1.end());
-            if(valency(m, v) < 3)
-                return false;
+        // get the intersection of the shared vertices from both rings
+        vector<VertexID> lisect;
-            // remove the vertex shared by the two rings from the intersection
+        std::back_insert_iterator<vector<VertexID> > lii(lisect);
-            vector<VertexID>::iterator iter;
+        set_intersection(link0.begin(), link0.end(),
-            iter = find(lisect.begin(), lisect.end(), v);
+                         link1.begin(), link1.end(),
-            assert(iter != lisect.end());
+                         lii);
-            lisect.erase(iter);
-            ++k;
+        // sort the faces of the two rings
-        }
+        sort(faces0.begin(), faces0.end());
-        // if the adjacent face is a triangle (see 2)
+        sort(faces1.begin(), faces1.end());
-        if(hew.opp().next().next().next().halfedge() == hew.opp().halfedge()){
-            VertexID v = hew.opp().next().vertex();
+        // get the intersection of the shared faces from both rings
+        vector<FaceID> fisect;
-            // valency test (see 5)
+        std::back_insert_iterator<vector<FaceID> > fii(fisect);
-            if(valency(m, v) < 3)
+        set_intersection(faces0.begin(), faces0.end(),
-                return false;
+                         faces1.begin(), faces1.end(),
+                         fii);
-            // remove the vertex shared by the two rings from the intersection
+        if(fisect.size() > 0)
-            vector<VertexID>::iterator iter;
+            return false;
-            iter = find(lisect.begin(), lisect.end(), v);
-            assert(iter != lisect.end());
+         k = 0;
-            lisect.erase(iter);
+        // if the adjacent face is a triangle (see 2)
-            ++k;
+        if(hew.next().next().next().halfedge() == h){
-        }
+            VertexID v = hew.next().vertex();
-        // double edge test (see 3)
-        if(lisect.size() != 0)
+            // valency test (see 5)
-            return false;
+            if(valency(m, v) < 3)
+                return false;
-        // tetrahedon test (see 4)
-        if(k == 2 && (link0.size() + link1.size() == 6))
+            // remove the vertex shared by the two rings from the intersection
-            return false;
+            vector<VertexID>::iterator iter;
+            iter = find(lisect.begin(), lisect.end(), v);
-        // test that we do not merge holes (see 6)
+            assert(iter != lisect.end());
-        if(boundary(m, v0) && boundary(m, v1) && hew.face() != InvalidFaceID && hew.opp().face() != InvalidFaceID)
+            lisect.erase(iter);
-            return false;
+            ++k;
+        }
-        return true;
+        // if the adjacent face is a triangle (see 2)
-    }
+        if(hew.opp().next().next().next().halfedge() == hew.opp().halfedge()){
+            VertexID v = hew.opp().next().vertex();
-    bool precond_flip_edge(const Manifold& m, HalfEdgeID h)
-    {
+            // valency test (see 5)
-        Walker j = m.walker(h);
+            if(valency(m, v) < 3)
+                return false;
-        FaceID hf = j.face();
-        FaceID hof = j.opp().face();
+            // remove the vertex shared by the two rings from the intersection
+            vector<VertexID>::iterator iter;
-        // boundary case
+            iter = find(lisect.begin(), lisect.end(), v);
-        if(hf == InvalidFaceID || hof == InvalidFaceID)
+            assert(iter != lisect.end());
-           return false;
+            lisect.erase(iter);
+            ++k;
+        }
-        // We can only flip an edge if both incident polygons are triangles.
+        // double edge test (see 3)
-        if(no_edges(m, hf) != 3 || no_edges(m, hof) !=3)
+        if(lisect.size() != 0)
             return false;
+        // tetrahedon test (see 4)
-        // non boundary vertices with a valency of less than 4(less than 3 after operation) degenerates mesh.
+        if(k == 2 && (link0.size() + link1.size() == 6))
-        VertexID hv = j.vertex();
+            return false;
-        VertexID hov = j.opp().vertex();
-        if((valency(m, hv) < 4 && !boundary(m, hv)) || (valency(m, hov) < 4 && !boundary(m, hov))){
+        // test that we do not merge holes (see 6)
-            return false;
+        if(boundary(m, v0) && boundary(m, v1) && hew.face() != InvalidFaceID && hew.opp().face() != InvalidFaceID)
-        }
+            return false;
-        // Disallow flip if vertices being connected already are.
+        return true;
-        VertexID hnv = j.next().vertex();
+    }
-        VertexID honv = j.opp().next().vertex();
-        if(connected(m, hnv, honv)){
+    bool precond_flip_edge(const Manifold& m, HalfEdgeID h)
-           return false;
+    {
-        }
+        Walker j = m.walker(h);
-        return true;
+        FaceID hf = j.face();
-    }
+        FaceID hof = j.opp().face();
-    bool boundary(const Manifold& m, VertexID v)
+        // boundary case
-    {
+        if(hf == InvalidFaceID || hof == InvalidFaceID)
-        Walker j  = m.walker(v);
+           return false;
-        return boundary(m, j.halfedge());
-    }
+        // We can only flip an edge if both incident polygons are triangles.
-    int valency(const Manifold& m, VertexID v)
+        if(no_edges(m, hf) != 3 || no_edges(m, hof) !=3)
-    {
+            return false;
-        return circulate_vertex_ccw(m,v, (std::function<void(Walker&)>)[](Walker){});
-    }
+        // non boundary vertices with a valency of less than 4(less than 3 after operation) degenerates mesh.
-    Manifold::Vec normal(const Manifold& m, VertexID v)
+        VertexID hv = j.vertex();
-    {
+        VertexID hov = j.opp().vertex();
-        Manifold::Vec p0 = m.pos(v);
+        if((valency(m, hv) < 4 && !boundary(m, hv)) || (valency(m, hov) < 4 && !boundary(m, hov))){
-        vector<Manifold::Vec> one_ring;
+            return false;
+        }
-        // run through outgoing edges, and store them normalized
-        circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID vn) {
+        // Disallow flip if vertices being connected already are.
-            Manifold::Vec edge = m.pos(vn) - p0;
+        VertexID hnv = j.next().vertex();
-            double l = length(edge);
+        VertexID honv = j.opp().next().vertex();
-            if(l > 0.0)
+        if(connected(m, hnv, honv)){
-                one_ring.push_back(edge/l);
+           return false;
-        });
+        }
-        int N = one_ring.size();
-        if(N<2)
+        return true;
-            return Manifold::Vec(0);
+    }
-        size_t N_count = N;
+    bool boundary(const Manifold& m, VertexID v)
-        size_t N_start = 0;
+    {
-        if(boundary(m, v))
+        Walker j  = m.walker(v);
-            N_start = 1;
+        return boundary(m, j.halfedge());
+    }
-        // sum up the normals of each face surrounding the vertex
-        Manifold::Vec n(0);
+    int valency(const Manifold& m, VertexID v)
-        for(size_t i = N_start; i < N_count; ++i){
+    {
-            Manifold::Vec e0 = one_ring[i];
+        return circulate_vertex_ccw(m,v, (std::function<void(Walker&)>)[](Walker){});
-            Manifold::Vec e1 = one_ring[(i+1) % N];
+    }
-            Manifold::Vec n_part = normalize(cross(e0, e1));
+    Manifold::Vec normal(const Manifold& m, VertexID v)
-            n += n_part * acos(max(-1.0, min(1.0, dot(e0, e1))));
+    {
-        }
+        Manifold::Vec p0 = m.pos(v);
+        vector<Manifold::Vec> one_ring;
-        // normalize and return the normal
-        float sqr_l = sqr_length(n);
+        // run through outgoing edges, and store them normalized
-        if(sqr_l > 0.0f) return n / sqrt(sqr_l);
+        circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID vn) {
+            Manifold::Vec edge = m.pos(vn) - p0;
-        return n;
+            double l = length(edge);
-    }
+            if(l > 0.0)
+                one_ring.push_back(edge/l);
+        });
-    bool connected(const Manifold& m, VertexID v0, VertexID v1)
+        int N = one_ring.size();
-    {
+        if(N<2)
-        bool c=false;
+            return Manifold::Vec(0);
-        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID v){ c |= (v==v1);});
-        return c;
+        size_t N_count = N;
-    }
+        size_t N_start = 0;
+        if(boundary(m, v))
+            N_start = 1;
-    int no_edges(const Manifold& m, FaceID f)
-    {
+        // sum up the normals of each face surrounding the vertex
-        return circulate_face_ccw(m, f, (std::function<void(Walker&)>)[](Walker w){});
+        Manifold::Vec n(0);
-    }
+        for(size_t i = N_start; i < N_count; ++i){
+            Manifold::Vec e0 = one_ring[i];
-    Manifold::Vec normal(const Manifold& m, FaceID f)
+            Manifold::Vec e1 = one_ring[(i+1) % N];
-    {
-        vector<Manifold::Vec> v;
+            Manifold::Vec n_part = normalize(cross(e0, e1));
+            n += n_part * acos(max(-1.0, min(1.0, dot(e0, e1))));
-        int k= circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID vid) {
+        }
-            v.push_back(m.pos(vid));
-        });
+        // normalize and return the normal
+        float sqr_l = sqr_length(n);
-        Manifold::Vec norm(0);
+        if(sqr_l > 0.0f) return n / sqrt(sqr_l);
-        for(int i=0;i<k;++i)
-        {
+        return n;
-            norm[0] += (v[i][1]-v[(i+1)%k][1])*(v[i][2]+v[(i+1)%k][2]);
+    }
-            norm[1] += (v[i][2]-v[(i+1)%k][2])*(v[i][0]+v[(i+1)%k][0]);
-            norm[2] += (v[i][0]-v[(i+1)%k][0])*(v[i][1]+v[(i+1)%k][1]);
-        }
+    bool connected(const Manifold& m, VertexID v0, VertexID v1)
-        float l = sqr_length(norm);
+    {
-        if(l>0.0f)
+        bool c=false;
-            norm /= sqrt(l);
+        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID v){ c |= (v==v1);});
-        return norm;
+        return c;
     }
-    double area(const Manifold& m, FaceID fid)
+    int no_edges(const Manifold& m, FaceID f)
     {
-        // Get all projected vertices
+        return circulate_face_ccw(m, f, (std::function<void(Walker&)>)[](Walker w){});
-        vector<Manifold::Vec> vertices;
+    }
-        int N = circulate_face_ccw(m, fid, (std::function<void(VertexID)>)[&](VertexID vid) {
-            vertices.push_back(m.pos(vid));
+    Manifold::Vec normal(const Manifold& m, FaceID f)
-        });
+    {
+        vector<Manifold::Vec> v;
-        double area = 0;
+        int k= circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID vid) {
-        Manifold::Vec norm = normal(m,fid);
+            v.push_back(m.pos(vid));
-        for(int i = 1; i < N-1; ++i)
+        });
-            area += 0.5 * dot(norm,cross(vertices[i]-vertices[0], vertices[(i+1 )]-vertices[0]));
-        return area;
+        Manifold::Vec norm(0);
-    }
+        for(int i=0;i<k;++i)
+        {
-    Manifold::Vec centre(const Manifold& m, FaceID f)
+            norm[0] += (v[i][1]-v[(i+1)%k][1])*(v[i][2]+v[(i+1)%k][2]);
-    {
+            norm[1] += (v[i][2]-v[(i+1)%k][2])*(v[i][0]+v[(i+1)%k][0]);
-        Manifold::Vec c(0);
+            norm[2] += (v[i][0]-v[(i+1)%k][0])*(v[i][1]+v[(i+1)%k][1]);
-        int n = circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID v) {c+=m.pos(v);});
+        }
-        return c / n;
+        float l = sqr_length(norm);
-    }
+        if(l>0.0f)
+            norm /= sqrt(l);
-    double perimeter(const Manifold& m, FaceID f)
+        return norm;
-    {
+    }
-        double l=0.0;
-        circulate_face_ccw(m, f, (std::function<void(HalfEdgeID)>)[&](HalfEdgeID h) { l+= length(m, h);});
-        return l;
+    double area(const Manifold& m, FaceID fid)
-    }
+    {
+        // Get all projected vertices
-    bool boundary(const Manifold& m, HalfEdgeID h)
+        vector<Manifold::Vec> vertices;
-    {
+        int N = circulate_face_ccw(m, fid, (std::function<void(VertexID)>)[&](VertexID vid) {
-        Walker w = m.walker(h);
+            vertices.push_back(m.pos(vid));
-        return w.face() == InvalidFaceID || w.opp().face() == InvalidFaceID;
+        });
-    }
-    double length(const Manifold& m, HalfEdgeID h)
+        double area = 0;
-    {
+        Manifold::Vec norm = normal(m,fid);
-        Walker w = m.walker(h);
+        for(int i = 1; i < N-1; ++i)
-        return (m.pos(w.vertex()) - m.pos(w.opp().vertex())).length();
+            area += 0.5 * dot(norm,cross(vertices[i]-vertices[0], vertices[(i+1 )]-vertices[0]));
-    }
+        return area;
-}
+    }
+    Manifold::Vec centre(const Manifold& m, FaceID f)
+    {
+        Manifold::Vec c(0);
+        int n = circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID v) {c+=m.pos(v);});
+        return c / n;
+    }
+    double perimeter(const Manifold& m, FaceID f)
+    {
+        double l=0.0;
+        circulate_face_ccw(m, f, (std::function<void(HalfEdgeID)>)[&](HalfEdgeID h) { l+= length(m, h);});
+        return l;
+    }
+    bool boundary(const Manifold& m, HalfEdgeID h)
+    {
+        Walker w = m.walker(h);
+        return w.face() == InvalidFaceID || w.opp().face() == InvalidFaceID;
+    }
+    double length(const Manifold& m, HalfEdgeID h)
+    {
+        Walker w = m.walker(h);
+        return (m.pos(w.vertex()) - m.pos(w.opp().vertex())).length();
+    }
+}

Subversion Repositories gelsvn

(root)/branches/GEL_github/src/HMesh/Manifold.cpp @ 677 – Rev 667 → 677