Subversion Repositories gelsvn

/* ----------------------------------------------------------------------- *

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

                }

            }

            bool operator==(const EdgeKey& k2) const {return k2.v0 == v0 && k2.v1==v1;}

            size_t hash() const {return v0.get_index()*3125*49+v1.get_index()*3125+7;}

        private:

            VertexID v0;

            VertexID v1;

        };

        struct Edge

        {

            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

        vector<int> faces(mesh.geometry.no_faces(), 3);

        build_template( static_cast<size_t>(mesh.geometry.no_vertices()),

					   reinterpret_cast<const float*>(&mesh.geometry.vertex(0)),

					   static_cast<size_t>(faces.size()),

					   static_cast<const int*>(&faces[0]),

					   reinterpret_cast<const int*>(&mesh.geometry.face(0)));

    }

    void Manifold::build(   size_t no_vertices,

						 const float* vertvec,

						 size_t no_faces,

						 const int* facevec,

						 const int* indices)

    {

        build_template(no_vertices, vertvec, no_faces, facevec, indices);

    }

    void Manifold::build(   size_t no_vertices,

						 const double* vertvec,

						 size_t no_faces,

						 const int* facevec,

						 const int* indices)

    {

        build_template(no_vertices, vertvec, no_faces, facevec, indices);

    }

    FaceID Manifold::add_face(std::vector<Manifold::Vec> points)

    {

        int F = points.size();

        vector<int> indices;

        for(size_t i=0;i<points.size(); ++i)

            indices.push_back(i);

        FaceID fid = *faces_end();

        build(points.size(), reinterpret_cast<double*>(&points[0]), 1, &F, &indices[0]);

        return fid;

    }

    bool Manifold::remove_face(FaceID fid)

    {

        if(!in_use(fid))

            return false;

        HalfEdgeID he = kernel.last(fid);

        HalfEdgeID last = he;

        vector<HalfEdgeID> halfedges;

        do

        {

            halfedges.push_back(he);

            kernel.set_face(he, InvalidFaceID);

            he = kernel.next(he);

        }

        while(he != last);

        vector<HalfEdgeID> halfedges_garbage;

        vector<VertexID> vertices_garbage;

        for(size_t i=0;i<halfedges.size(); ++i)

        {

            HalfEdgeID h = halfedges[i];

            HalfEdgeID hopp = kernel.opp(h);

            if(kernel.face(hopp) == InvalidFaceID)

            {

                halfedges_garbage.push_back(h);

                VertexID v0 = kernel.vert(hopp);

                if(valency(*this, v0) <= 2 )

                    vertices_garbage.push_back(v0);

                else {

                    link(kernel.prev(h), kernel.next(hopp));

                    kernel.set_out(v0, kernel.opp(kernel.prev(h)));

                }

                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)

            kernel.remove_vertex(vertices_garbage[i]);

        kernel.remove_face(fid);

        return true;

    }

    bool Manifold::remove_edge(HalfEdgeID hid)

    {

        if(!in_use(hid))

            return false;

        FaceID f0 = kernel.face(hid);

        FaceID f1 = kernel.face(kernel.opp(hid));

        remove_face(f0);

        remove_face(f1);

        return true;

    }

    bool Manifold::remove_vertex(VertexID vid)

    {

        if(!in_use(vid))

            return false;

        vector<FaceID> faces;

        int N = circulate_vertex_ccw(*this, vid, (std::function<void(FaceID)>)[&](FaceID f) {

            faces.push_back(f);

        });

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

        HalfEdgeID hn = kernel.next(h);

        HalfEdgeID hp = kernel.prev(h);

        HalfEdgeID hon = kernel.next(ho);

        HalfEdgeID hop = kernel.prev(ho);

		FaceID f = kernel.face(h);

		FaceID fo = kernel.face(ho);

        // average the vertex positions

        pos(hv) = avg_vertices ? (0.5f * (pos(hov) + pos(hv))) : pos(hv);

        // update all halfedges pointing to hov to point to hv, effectively removing hov from all loops

        HalfEdgeID he = kernel.out(hov);

        HalfEdgeID last = he;

		do {

			assert(kernel.vert(kernel.opp(he)) == hov);

            kernel.set_vert(kernel.opp(he), hv);

            he = kernel.next(kernel.opp(he));

        }

		while(he != last);

        kernel.set_out(hv, hn);

        // link hp and hn, effectively removing h from opposite loop

        link(hp, hn);

        // make face owning h own hn instead

        if(kernel.face(h) != InvalidFaceID)

            kernel.set_last(f, hn);

        // link hop and hon, effectively removing h from opposite face loop

        link(hop, hon);

        // make opposite face owning h own hon instead

        if(kernel.face(ho) != InvalidFaceID)

            kernel.set_last(fo, hon);

        // remove the obsolete entities

        kernel.remove_vertex(hov);

        kernel.remove_halfedge(h);

        kernel.remove_halfedge(ho);

        // verify that remaining faces haven't become degenerate because of collapse

        remove_face_if_degenerate(hn);

        remove_face_if_degenerate(hon);

        // verify that v is sane after collapse

        ensure_boundary_consistency(hv);

    }

    FaceID Manifold::split_face_by_edge(FaceID f, VertexID v0, VertexID v1)

    {

		if(connected(*this, v0, v1))

            return InvalidFaceID;

        HalfEdgeID he = kernel.last(f);

        HalfEdgeID last = he;

        int steps = 0;

        while(he != last || steps == 0){

            ++steps;

            he = kernel.next(he);

        }

        // make sure this is not a triangle

        assert(steps > 3);

        // make sure we are not trying to connect a vertex to itself

        assert(v0 != v1);

		HalfEdgeID h0 = kernel.out(v0);

		for(Walker w = walker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){

			if(w.face() == f){

				h0 = w.halfedge();

				break;

			}

		}

		assert(kernel.face(h0) != InvalidFaceID);

        assert(kernel.face(h0) == f);

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

        HalfEdgeID h = h0;

        while(kernel.vert(h) != v1){

            h = kernel.next(h);

        }

        assert(h != h0);

        // create a new halfedge ha which connects v1 and v0 closing the first loop.

        HalfEdgeID h1 = kernel.next(h);

        HalfEdgeID ha = kernel.add_halfedge();

        link(h, ha);

        link(ha, h0);

        kernel.set_face(ha, f);

        kernel.set_vert(ha, v0);

        kernel.set_last(f, ha);

        // create a new face, f2, and set all halfedges in the remaining part of the polygon to point to this face.

        h = h1;

        FaceID f2 = kernel.add_face();

        while(kernel.vert(h) != v0){

            kernel.set_face(h, f2);

            h = kernel.next(h);

        }

        kernel.set_face(h, f2);

        assert(h != h1);

        // create a new halfedge hb to connect v0 and v1.

        HalfEdgeID hb = kernel.add_halfedge();

        link(h, hb);

        link(hb, h1);

        kernel.set_face(hb, f2);

        kernel.set_vert(hb, v1);

        kernel.set_last(f2, hb);

        // complete the operation by gluing the two new halfedges

        glue(ha, hb);

        // assert sanity of operation

        assert(kernel.next(kernel.opp(kernel.prev(h1))) == h0);

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

        assert(kernel.face(hb) == f2);

        // return handle to newly created face

        return f2;

    }

    VertexID Manifold::split_face_by_vertex(FaceID f)

    {

        //create the new vertex, with the barycenter of the face as position

        Manifold::Vec p(0.0f);

        HalfEdgeID last_he = kernel.last(f);

        HalfEdgeID he = last_he;

        int steps = 0;

        while(he != last_he || steps == 0){

            p += positions[kernel.vert(he)];

            ++steps;

            he = kernel.next(he);

        }

        p /= steps;

        VertexID v = kernel.add_vertex();

        positions[v]  = p;

        //circulate the face, create halfedges and connect to vertex

        vector<HalfEdgeID> eout;

        vector<HalfEdgeID> ein;

        last_he = kernel.last(f);

        he = last_he;

        do{

            HalfEdgeID hn = kernel.next(he);

            HalfEdgeID ho = kernel.add_halfedge();

            HalfEdgeID hi = kernel.add_halfedge();

            FaceID fn = kernel.add_face();

            kernel.set_face(hi, fn);

            kernel.set_vert(hi, v);

            kernel.set_face(ho, fn);

            kernel.set_vert(ho, kernel.vert(kernel.opp(he)));

            kernel.set_face(he, fn);

            kernel.set_last(fn, ho);

            link(hi, ho);

            link(ho, he);

            link(he, hi);

            eout.push_back(ho);

            ein.push_back(hi);

            he = hn;

        }

        while(he != last_he);

        // glue new halfedges together

        size_t N = eout.size();

        for(size_t i = 0; i < N; ++i){

            glue(ein[i], eout[(i+1)%N]);

        }

        kernel.set_out(v, eout[0]);

        //remove the now replaced face

        kernel.remove_face(f);

        return v;

    }

    VertexID Manifold::split_edge(HalfEdgeID h)

    {

        HalfEdgeID ho = kernel.opp(h);

        VertexID v = kernel.vert(h);

        VertexID vo = kernel.vert(ho);

        //create the new vertex with middle of edge as position and update connectivity

        VertexID vn = kernel.add_vertex();

        positions[vn] = .5f * (positions[v] + positions[vo]);

        kernel.set_out(vn, h);

        //create two necessary halfedges, and update connectivity

        HalfEdgeID hn = kernel.add_halfedge();

        HalfEdgeID hno = kernel.add_halfedge();

        kernel.set_out(vo, hn);

        kernel.set_out(v, hno);

        glue(h, hno);

        glue(hn, ho);

        link(kernel.prev(h), hn);

        link(hn, h);

        kernel.set_vert(hn, vn);

        link(kernel.prev(ho), hno);

        link(hno, ho);

        kernel.set_vert(hno, vn);

        // update faces in case of non boundary edge

        if(kernel.face(h) != InvalidFaceID)

            kernel.set_last(kernel.face(h), hn);

        if(kernel.face(ho) != InvalidFaceID)

            kernel.set_last(kernel.face(ho), ho);

        // update new edge with faces from existing edge

        kernel.set_face(hn, kernel.face(h));

        kernel.set_face(hno, kernel.face(ho));

        //if split occurs on a boundary, consistency must be ensured.

        ensure_boundary_consistency(vn);

        ensure_boundary_consistency(v);

        ensure_boundary_consistency(vo);

        return vn;

    }

    size_t link_intersection(const Manifold& m, VertexID v0, VertexID v1, vector<VertexID>& lisect)

    {

        // get the one-ring of v0

        vector<VertexID> link0;

        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID vn) {

            link0.push_back(vn);

        });

        // get the one-ring of v1

        vector<VertexID> link1;

        circulate_vertex_ccw(m, v1, (std::function<void(VertexID)>)[&](VertexID vn) {

            link1.push_back(vn);

        });

        // sort the vertices of the two rings

        sort(link0.begin(), link0.end());

        sort(link1.begin(), link1.end());

        // get the intersection of the shared vertices from both rings

        std::back_insert_iterator<vector<VertexID> > lii(lisect);

        set_intersection(link0.begin(), link0.end(),

						 link1.begin(), link1.end(),

						 lii);

        return lisect.size();

    }

    bool Manifold::stitch_boundary_edges(HalfEdgeID h0, HalfEdgeID h1)

    {

        // Cannot stitch an edge with itself

        if(h0 == h1)

            return false;

        // Only stitch a pair of boundary edges.

        if(kernel.face(h0) == InvalidFaceID && kernel.face(h1) == InvalidFaceID)

        {

            HalfEdgeID h0o = kernel.opp(h0);

            HalfEdgeID h1o = kernel.opp(h1);

            VertexID v0a = kernel.vert(h0);

            VertexID v0b = kernel.vert(kernel.opp(h1));

            VertexID v1b = kernel.vert(h1);

            VertexID v1a = kernel.vert(kernel.opp(h0));

            // Case below implies that h0 and h1 are the same edge with different ID

            // That should not happen.

            assert(!(v0a == v0b && v1a == v1b));

            if(connected(*this, v0a, v0b))

                return false;

            if(connected(*this, v1a, v1b))

                return false;

            if(v0b != v0a)

            {

                // Check the link intersection v0a and v0b are welded together

                // 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(link_intersection(*this, v0a, v0b, lisect))

                {

                    vector<VertexID>::iterator iter;

                    if(v1a == v1b)

                    {

                        iter = find(lisect.begin(), lisect.end(), v1a);

                        if(iter != lisect.end())

                            lisect.erase(iter);

                    }

                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h0)));

                    if(iter != lisect.end())

                        lisect.erase(iter);

                    if(lisect.size() > 0)

                        return false;

                }

            }

            if(v1a != v1b)

            {

                // Check the same for the other endpoints.

                vector<VertexID> lisect;

                if(link_intersection(*this, v1a, v1b, lisect))

                {

                    vector<VertexID>::iterator iter;

                    if(v0a == v0b)

                    {

                        iter = find(lisect.begin(), lisect.end(), v0a);

                        if(iter != lisect.end())

                            lisect.erase(iter);

                    }

                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h1)));

                    if(iter != lisect.end())

                        lisect.erase(iter);

                    if(lisect.size() > 0)

                        return false;

                }

            }

            if(v0b != v0a)

                circulate_vertex_ccw(*this, v0b, (std::function<void(Walker&)>)[&](Walker hew) {

                    kernel.set_vert(hew.opp().halfedge(), v0a);

                });

            if(v1b != v1a)

                circulate_vertex_ccw(*this, v1b, (std::function<void(Walker&)>)[&](Walker hew) {

                    kernel.set_vert(hew.opp().halfedge(), v1a);

                });

            if(v0a != v0b)

            {

                HalfEdgeID h1p = kernel.prev(h1);

                HalfEdgeID h0n = kernel.next(h0);

                if(kernel.next(h0n) == h1p)

                {

                    glue(kernel.opp(h0n), kernel.opp(h1p));

                    kernel.set_out(kernel.vert(h0n),kernel.opp(h0n));

                    kernel.remove_halfedge(h0n);

                    kernel.remove_halfedge(h1p);

                }

                else

                    link(h1p, h0n);

                kernel.remove_vertex(v0b);

            }

            if(v1a != v1b)

            {

                HalfEdgeID h0p = kernel.prev(h0);

                HalfEdgeID h1n = kernel.next(h1);

                if(kernel.next(h1n) == h0p)

                {

                    glue(kernel.opp(h0p), kernel.opp(h1n));

                    kernel.set_out(kernel.vert(h1n), kernel.opp(h1n));

                    kernel.remove_halfedge(h0p);

                    kernel.remove_halfedge(h1n);

                }

                else

                    link(h0p, h1n);

                kernel.remove_vertex(v1b);

            }

            glue(h0o, h1o);

            kernel.remove_halfedge(h0);

            kernel.remove_halfedge(h1);

            kernel.set_out(v1a, h1o);

            kernel.set_out(v0a, h0o);

            ensure_boundary_consistency(v1a);

            ensure_boundary_consistency(v0a);

            return true;

        }

        return false;

    }

    FaceID Manifold::merge_one_ring(VertexID v, float max_loop_length)

    {

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

        // If the vertex is  a boundary vertex, we preparte the walker so that the

        // first face visited is not the invalid face outside the boundary. If the boundary

        // vertex is adjacent to only one vertex, there is little to do and we bail out.

        bool vertex_is_boundary = false;

        Walker hew = walker(v);

        if(boundary(*this, v))

        {

            if(val==2) return InvalidFaceID;

            vertex_is_boundary = true;

            hew = hew.circulate_vertex_ccw();

        }

        // Prepare some vectors for taking out the trash: We remove all old faces and all orphaned edges

        // and vertices

        vector<HalfEdgeID> garbage_halfedges;

        vector<FaceID> garbage_faces;

        vector<VertexID> garbage_vertices;

        vector<HalfEdgeID> loop; // The halfedges which form the outer loop of the merged one ring.

        // Below we loop over all faces adjacent to the vertex and add their halfedges to a big loop

        // which will form the loop of the new merged face. Below we remove from the loop edges

        // that appear twice (as each other's opposite).

        for(;!hew.full_circle(); hew = hew.circulate_vertex_ccw())

        {

            garbage_faces.push_back(hew.face());

            for(Walker hew2 = walker(hew.halfedge());

                !hew2.full_circle(); hew2 = hew2.circulate_face_ccw())

                loop.push_back(hew2.halfedge());

        }

        // Now merge the loops. We iteratively remove pairs of adjacent halfedges from the loop

        // if we find that the second is the opposite of the first since this is a degenerate

        // situation. However, we stop at four remaining halfedges since otherwise the loop degenerates

        // to zero after the next step if these four are also pairwise each other's opposites.

        int did_work;

        do

        {

            did_work = 0;

            vector<HalfEdgeID> loop_tmp(0);

            for(size_t i=0;i<loop.size();++i)

                if(walker(loop[i]).opp().halfedge() == loop[(i+1)%loop.size()])

                {

                    VertexID vid = walker(loop[i]).vertex();

                    if(vid != v)

                        garbage_vertices.push_back(walker(loop[i]).vertex());

                    garbage_halfedges.push_back(loop[i]);

                    garbage_halfedges.push_back(loop[(i+1)%loop.size()]);

                    loop[i] = InvalidHalfEdgeID;

                    loop[(i+1)%loop.size()] = InvalidHalfEdgeID;

                    ++did_work;

                    ++i;

                }

            for(size_t i=0;i<loop.size();++i)

                if(loop[i] != InvalidHalfEdgeID)

                    loop_tmp.push_back(loop[i]);

            loop = loop_tmp;

        } while(did_work > 0 && loop.size() > 4);

        // Check whether the loop is too long

        float loop_len=0.0;

        for(size_t i=0;i<loop.size();++i)

            loop_len += length(*this, loop[i]);

        if(loop_len > max_loop_length)

            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.

        vector<HalfEdgeID> loop_tmp = loop;

        sort(loop_tmp.begin(), loop_tmp.end());

        vector<HalfEdgeID>::iterator end_iter = unique(loop_tmp.begin(), loop_tmp.end());

        if(end_iter != loop_tmp.end())

            return InvalidFaceID;

        // Remove all faces and connected halfedges and the original vertex v.

        for(size_t i=0;i<garbage_vertices.size(); ++i)

            kernel.remove_vertex(garbage_vertices[i]);

        for(size_t i=0;i<garbage_faces.size(); ++i)

            kernel.remove_face(garbage_faces[i]);

        for(size_t i=0;i<garbage_halfedges.size(); ++i)

            kernel.remove_halfedge(garbage_halfedges[i]);

        if(!vertex_is_boundary)

            kernel.remove_vertex(v);

        // Create a new face for the merged one ring and link up all the halfedges

        // in the loop

        FaceID f = kernel.add_face();

        kernel.set_last(f,loop[0]);

        for(size_t i=0;i<loop.size(); ++i)

        {

            kernel.set_face(loop[i], f);

            Walker hw = walker(loop[i]);

            kernel.set_out(hw.vertex(),hw.opp().halfedge());

            link(loop[i],loop[(i+1)%loop.size()]);

            assert(hw.vertex() == walker(loop[(i+1)%loop.size()]).opp().vertex());

        }

        // Finally, we ensure boundary consitency for all vertices in the loop.

        for(size_t i=0;i<loop.size(); ++i)

        {

            Walker hw = walker(loop[i]);

            ensure_boundary_consistency(hw.vertex());

        }

        // Return the brand new merged face.

        return f;

    }

    bool Manifold::merge_faces(FaceID f, HalfEdgeID h)

    {

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

        HalfEdgeID hn = kernel.next(h);

        HalfEdgeID hon = kernel.next(ho);

        if(fo == f)

            return false;

        //boundary case

        if(kernel.vert(hn) == kernel.vert(hon))

            return false;

        HalfEdgeID hp = kernel.prev(h);

        HalfEdgeID hop = kernel.prev(ho);

        VertexID v = kernel.vert(h);

        VertexID vo = kernel.vert(ho);

        if(valency(*this, v) < 3 || valency(*this, vo) < 3)

            return false;

        link(hop, hn);

        link(hp, hon);

        kernel.set_out(vo, hon);

        kernel.set_out(v, hn);

        kernel.set_last(f, hn);

        HalfEdgeID hx = hon;

        assert(kernel.face(hx) == fo);

        while(kernel.face(hx) != f){

            kernel.set_face(hx, f);

            hx = kernel.next(hx);

        }

        ensure_boundary_consistency(v);

        ensure_boundary_consistency(vo);

        kernel.remove_halfedge(h);

        kernel.remove_halfedge(ho);

        kernel.remove_face(fo);

        return true;

    }

    FaceID Manifold::close_hole(HalfEdgeID h)

    {

        // invalid face is a hole

        if(kernel.face(h) == InvalidFaceID){

            FaceID f = kernel.add_face();

            kernel.set_last(f, h);

            do{

                kernel.set_face(h, f);

                h = kernel.next(h);

            }

            while(kernel.face(h) != f);

            return f;

        }

        return kernel.face(h);

    }

    VertexID Manifold::slit_vertex(VertexID v, HalfEdgeID h_in, HalfEdgeID h_out)

    {

        assert(kernel.face(h_in) != InvalidFaceID);

        assert(kernel.face(h_out) != InvalidFaceID);

        assert(kernel.opp(h_out) != h_in);

        VertexID v_new = kernel.add_vertex();

        pos(v_new) = pos(v);

        HalfEdgeID h = kernel.prev(h_out);

        kernel.set_vert(h, v_new);

        while ( h != h_in) {

            h = kernel.prev(kernel.opp(h));

            kernel.set_vert(h, v_new);

        }

        HalfEdgeID h_in_opp = kernel.opp(h_in);

        HalfEdgeID hn_in, hn_in_opp;

        if(kernel.face(h_in_opp) != InvalidFaceID)

        {

            hn_in = kernel.add_halfedge();

            kernel.set_face(hn_in, InvalidFaceID);

            glue(h_in_opp, hn_in);

            hn_in_opp = kernel.add_halfedge();

            kernel.set_face(hn_in_opp, InvalidFaceID);

            glue(h_in, hn_in_opp);

            link(hn_in_opp, hn_in);

            VertexID v_i = kernel.vert(h_in_opp);

            kernel.set_vert(hn_in_opp, v_i);

            kernel.set_out(v_i, hn_in);

        }

        else

        {

            hn_in_opp = h_in_opp;

            hn_in = kernel.prev(hn_in_opp);

            h_in_opp = kernel.opp(hn_in);

        }

        HalfEdgeID h_out_opp = kernel.opp(h_out);

        HalfEdgeID hn_out,hn_out_opp;

        if(kernel.face(h_out_opp) != InvalidFaceID)

        {

            hn_out = kernel.add_halfedge();

            kernel.set_face(hn_out, InvalidFaceID);

            glue(h_out_opp, hn_out);

            hn_out_opp = kernel.add_halfedge();

            kernel.set_face(hn_out_opp, InvalidFaceID);

            glue(h_out, hn_out_opp);

            link(hn_out, hn_out_opp);

            VertexID v_o = kernel.vert(h_out);

            kernel.set_vert(hn_out, v_o);

            kernel.set_out(v_o, hn_out_opp);

        }

        else

        {

            hn_out_opp = h_out_opp;

            hn_out = kernel.next(hn_out_opp);

            h_out_opp = kernel.opp(hn_out);

        }

        link(hn_out_opp, hn_in_opp);

        link(hn_in, hn_out);

        kernel.set_vert(hn_in, v);

        kernel.set_vert(hn_out_opp, v_new);

        kernel.set_out(v, hn_out);

        kernel.set_out(v_new, hn_in_opp);

        return v_new;

    }

    HalfEdgeID Manifold::slit_edges(VertexAttributeVector<int>& insel)

    {

        HalfEdgeID h;

        for(auto vid : vertices())

        {

            if(insel[vid])

            {

                HalfEdgeID h_in = InvalidHalfEdgeID, h_out = InvalidHalfEdgeID;

                Walker w = walker(vid);

                while(!w.full_circle())

                {

                    if(insel[w.vertex()]) {

                        if(h_in == InvalidHalfEdgeID) {

                            if(w.opp().face() == InvalidFaceID)

                                h_in = w.opp().next().opp().halfedge();

                            else

                                h_in = w.opp().halfedge();

                        }

                        else {

                            if(w.face() == InvalidFaceID)

                                h_out = w.prev().opp().halfedge();

                            else

                                h_out = w.halfedge();

                            break;

                        }

                    }

                    w = w.circulate_vertex_ccw();

                }