Subversion Repositories gelsvn

Rev

Rev 506 | Rev 508 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

#include "Manifold.h"

#include <iostream>
#include <vector>
#include <map>

#include <CGLA/Vec3f.h>
#include <Geometry/TriMesh.h>

#include "VertexHandle.h"
#include "FaceHandle.h"
#include "HalfEdgeHandle.h"

#include "VertexCirculator.h"
#include "FaceCirculator.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(IndexType va, IndexType 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(v0 > k2.v0){
                    return false;
                }
                else{
                    return v1 < k2.v1;
                }
            }
        private:
            int v0;
            int v1;         
        };

        struct Edge
        {
            HalfEdgeIndex h0;
            HalfEdgeIndex h1;
            IndexType count;
            Edge() : h0(NULL_HIDX), h1(NULL_HIDX), 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<int>(mesh.geometry.no_vertices()),
                        reinterpret_cast<const float*>(&mesh.geometry.vertex(0)),
                        static_cast<int>(faces.size()),
                        static_cast<int*>(&faces[0]),
                        reinterpret_cast<const int*>(&mesh.geometry.face(0)));
    }
    void Manifold::build(   unsigned int no_vertices,
                            const float* vertvec,
                            unsigned int no_faces,
                            const unsigned int* facevec,
                            const unsigned int* indices,
                            const int* touch)
    {
        build_template(no_vertices, vertvec, no_faces, facevec, indices, touch);
    }
 
    void Manifold::bbox(Vec3f& pmin, Vec3f& pmax) 
    {
        VertexHandle v = vhandles_begin();
        pmin = pmax = v.get().pos;
        ++v;
        for(; v != vhandles_end(); ++v){
            pmin = v_min(v.get().pos, pmin);
            pmax = v_max(v.get().pos, pmax);
        }
    }
    void Manifold::bsphere(Vec3f& c, float& r) 
    {
        Vec3f p0, p7;
        bbox(p0, p7);
        Vec3f rad = (p7 - p0) * 0.5f;
        c = p0 + rad;
        r = rad.length();
    }
    void Manifold::collapse_halfedge(HalfEdgeHandle handle, bool avg_vertices)
    {
        HalfEdgeIndex h = handle.get_idx();

        HalfEdgeIndex ho = halfedges[h].opp_idx;
        VertexIndex hv = halfedges[h].vert_idx;
        VertexIndex hov = halfedges[ho].vert_idx;

        // average the vertex positions
        vertices[hv].pos = avg_vertices * 0.5f * (vertices[hov].pos - vertices[hv].pos) + vertices[hv].pos;

        HalfEdgeIndex hn = halfedges[h].next_idx;
        HalfEdgeIndex hp = halfedges[h].prev_idx;
        HalfEdgeIndex hon = halfedges[ho].next_idx;
        HalfEdgeIndex hop = halfedges[ho].prev_idx;

        // update all halfedges pointing to hov to point to hv, effectively removing hov from all loops
        HalfEdgeIndex he = vertices[hov].out_idx;
        HalfEdgeIndex last = he;
        unsigned int steps = 0;
        while(he != last || steps == 0){ 
            halfedges[halfedges[he].opp_idx].vert_idx = hv;
            ++steps;
            he = halfedges[halfedges[he].opp_idx].next_idx;
        }   
        vertices[hv].out_idx = hn;

        // link hp and hn, effectively removing h from opposite loop
        link(hp, hn);
        // make face owning h own hn instead
        if(halfedges[h].face_idx != NULL_FIDX)
            faces[halfedges[h].face_idx].last_idx = hn;
       
        // link hop and hon, effectively removing h from opposite face loop
        link(hop, hon);
        // make opposite face owning h own hon instead
        if(halfedges[ho].face_idx != NULL_FIDX)
            faces[halfedges[ho].face_idx].last_idx = hon;
      
        // remove the obsolete entities
        remove_vertex(hov);
        remove_halfedge(h);
        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);
    }
    FaceHandle Manifold::split_face_by_edge(FaceHandle fhandle, VertexHandle vhandle0, VertexHandle vhandle1)
    {
        FaceIndex f = fhandle.get_idx();
        VertexIndex v0 = vhandle0.get_idx();
        VertexIndex v1 = vhandle1.get_idx();

        HalfEdgeIndex he = faces[f].last_idx;
        HalfEdgeIndex last = he;
        unsigned int steps = 0;
        while(he != last || steps == 0){
            ++steps;
            he = halfedges[halfedges[he].opp_idx].next_idx;
        }

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

        // find the halfege emanating from v0 which belongs to face being split
        he = vertices[v0].out_idx;
        last = he;
        steps = 0;
        while(halfedges[he].face_idx != f && !(last != he || steps == 0)){ 
            ++steps;
            he = halfedges[halfedges[he].opp_idx].next_idx;
        }
        assert(!(last != he || steps == 0));

        HalfEdgeIndex h0 = he;
        assert(halfedges[h0].face_idx == 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.
        HalfEdgeIndex h = h0;
        while(halfedges[h].vert_idx != v1){
            h = halfedges[h].next_idx;
        }
        assert(h != h0);

        // create a new halfedge ha which connects v1 and v0 closing the first loop.
        HalfEdgeIndex h1 = halfedges[h].next_idx;
        HalfEdgeIndex ha = add_halfedge(HalfEdge());
        link(h, ha);
        link(ha, h0);
        halfedges[ha].face_idx = f;
        halfedges[ha].vert_idx = v0;
        faces[f].last_idx = ha;

        // create a new face, f2, and set all halfedges in the remaining part of the polygon to point to this face.
        h = h1;
        FaceIndex f2 = add_face(Face());
        while(halfedges[h].vert_idx != v0){
            halfedges[h].face_idx = f;
            h = halfedges[h].next_idx;
        }
        halfedges[h].face_idx = f2;
        assert(h != h1);

        // create a new halfedge hb to connect v0 and v1.
        HalfEdgeIndex hb = add_halfedge(HalfEdge());
        link(h, hb);
        link(hb, h1);
        halfedges[hb].face_idx = f2;
        halfedges[hb].vert_idx = v1;
        faces[f2].last_idx = hb;

        // complete the operation by gluing the two new halfedges
        glue(ha, hb);

        // assert sanity of operation
        assert(halfedges[halfedges[halfedges[h1].prev_idx].opp_idx].next_idx == h0);
        assert(halfedges[halfedges[halfedges[h0].prev_idx].opp_idx].next_idx == h1);
        assert(halfedges[halfedges[hb].next_idx].face_idx == f2);
        assert(halfedges[halfedges[halfedges[hb].next_idx].next_idx].face_idx == f2);
        assert(halfedges[hb].face_idx == f2);

        // return handle to newly created face
        return FaceHandle(*this, f2);
    }
    VertexHandle Manifold::split_face_by_vertex(FaceHandle fhandle)
    {
        FaceIndex f = fhandle.get_idx();

        //create the new vertex, with the barycenter of the face as position
        Vec3f p(0.0f);
        HalfEdgeIndex last = faces[f].last_idx;
        HalfEdgeIndex he = last;
        unsigned int steps = 0;
        
        while(he != last || steps == 0){
            p += vertices[halfedges[he].vert_idx].pos;
            ++steps;
            he = halfedges[halfedges[he].opp_idx].next_idx;
        }
        p /= steps;
        VertexIndex v = add_vertex(Vertex(p));

        //circulate the face, create halfedges and connect to vertex
        vector<HalfEdgeIndex> eout;
        vector<HalfEdgeIndex> ein;
        last = faces[f].last_idx;
        he = last;

        do{
            HalfEdgeIndex hn = halfedges[he].next_idx;

            HalfEdgeIndex ho = add_halfedge(HalfEdge());
            HalfEdgeIndex hi = add_halfedge(HalfEdge());

            FaceIndex fn = add_face(Face());
            halfedges[hi].face_idx = fn;
            halfedges[hi].vert_idx = v;
            halfedges[ho].face_idx = fn;
            halfedges[ho].vert_idx = halfedges[halfedges[he].opp_idx].vert_idx;
            halfedges[he].face_idx = fn;
            faces[fn].last_idx = ho;

            link(hi, ho);
            link(ho, he);
            link(he, hi);

            eout.push_back(ho);
            ein.push_back(hi);

            he = hn;
        }
        while(he != last);

        // glue new halfedges together
        IndexType N = eout.size();
        for(IndexType i = 0; i < N; ++i){
            glue(ein[i], eout[(i+1)%N]);
        }
        vertices[v].out_idx = eout[0];

        //remove the now replaced face
        remove_face(f);

        return VertexHandle(*this, v);
    }
    VertexHandle Manifold::split_edge(HalfEdgeHandle hhandle)
    {
        HalfEdgeIndex h = hhandle.get_idx();
      
        HalfEdgeIndex ho = halfedges[h].opp_idx;
        VertexIndex v = halfedges[h].vert_idx;
        VertexIndex vo = halfedges[ho].vert_idx;

        //create the new vertex with middle of edge as position and update connectivity
        Vec3f np = 0.5f * (vertices[v].pos + vertices[vo].pos);
        VertexIndex vn = add_vertex(Vertex(np));
        vertices[vn].out_idx = h;

        //create two necessary halfedges, and update connectivity
        HalfEdgeIndex hn = add_halfedge(HalfEdge());
        HalfEdgeIndex hno = add_halfedge(HalfEdge());
        vertices[vo].out_idx = hn;
        vertices[v].out_idx = ho;
        glue(hn, hno);
        link(halfedges[h].prev_idx, hn);
        link(hn, h);
        halfedges[hn].vert_idx = vn;
        link(hno, halfedges[ho].next_idx);
        link(ho, hno);
        halfedges[hno].vert_idx = vo;
        halfedges[ho].vert_idx = vn;

        // update faces in case of non boundary edge
        if(halfedges[h].face_idx != NULL_FIDX)
            faces[halfedges[h].face_idx].last_idx = hn;  
        if(halfedges[ho].face_idx != NULL_FIDX)
            faces[halfedges[ho].face_idx].last_idx = ho;
        
        // update new edge with faces from existing edge
        halfedges[hn].face_idx = halfedges[h].face_idx;
        halfedges[hno].face_idx = halfedges[ho].face_idx;

        //if split occurs on a boundary, consistency must be ensured.
        ensure_boundary_consistency(vn);
        ensure_boundary_consistency(v);
        ensure_boundary_consistency(vo);

        return VertexHandle(*this, vn);
    }
    bool Manifold::merge_faces(FaceHandle fhandle, HalfEdgeHandle hhandle)
    {
        FaceIndex f = fhandle.get_idx();
        HalfEdgeIndex h = hhandle.get_idx();

        //assert that we're merging a valid face with the corresponding halfedge
        assert(halfedges[h].face_idx == f);

        HalfEdgeIndex ho = halfedges[h].opp_idx;
        FaceIndex fo = halfedges[ho].face_idx;
        HalfEdgeIndex hn = halfedges[h].next_idx;
        HalfEdgeIndex hon = halfedges[ho].next_idx;

        //boundary case
        if(halfedges[hn].vert_idx == halfedges[hon].vert_idx)
            return false;
        
        HalfEdgeIndex hp = halfedges[h].prev_idx;
        HalfEdgeIndex hop = halfedges[ho].prev_idx;
        VertexIndex v = halfedges[h].vert_idx;
        VertexIndex vo = halfedges[ho].vert_idx;

        link(hop, hn);
        vertices[v].out_idx = hn;

        link(hp, hon);
        vertices[vo].out_idx = hn;

        faces[f].last_idx = hn;
        HalfEdgeIndex hx = hon;
        assert(faces[hx].face_idx == fo);
        while(halfedges[h].face_idx != f){
            halfedges[hx].face_idx = f;
            hx = halfedges[hx].next_idx;
        }

        ensure_boundary_consistency(v);
        ensure_boundary_consistency(vo);

        remove_halfedge(h);
        remove_halfedge(ho);
        remove_face(fo);

        return true;
    }
    void Manifold::close_hole(HalfEdgeHandle hhandle)
    {
        HalfEdgeIndex h = hhandle.get_idx();
        
        if(halfedges[h].face_idx == NULL_FIDX){
            FaceIndex f = add_face(Face());
            faces[f].last_idx  = h;
            do{
                halfedges[h].face_idx = f;
                h = halfedges[h].next_idx;
            }
            while(halfedges[h].face_idx != f);
        }
    }
    bool Manifold::flip(HalfEdgeHandle h)
    {      
        FaceHandle f = h.face();
        HalfEdgeHandle ho = h.opp();
        FaceHandle fo = ho.face();

        // boundary case
        if(f == NULL_FHANDLE || fo == NULL_FHANDLE)
            return false;
        
                // We can only flip an edge if both incident polygons are triangles. 
        if(no_edges(f) != 3 || no_edges(fo) !=3){
                        return false;
        }
                // non boundary vertices with a valency of less than 4(less than 3 after operation) degenerates mesh.
        VertexHandle v = h.vert();
        ConstVertexHandle cv(v);
        VertexHandle vo = ho.vert();
        ConstVertexHandle cvo(vo);
        if((valency(v) < 4 && !boundary(v)) || (valency(vo) < 4 && !boundary(vo))){
            return false;
        }
                // Disallow flip if vertices being connected already are.
        VertexHandle vn = h.next().vert();
        VertexHandle von = ho.next().vert();
        if(connected(vn, von)){
            return false;
        }
        HalfEdgeHandle hn = h.next();
        HalfEdgeHandle hp = h.prev();
        HalfEdgeHandle hon = ho.next();
        HalfEdgeHandle hop = ho.prev();

        hop.get().face_idx = f.get_idx();
        hp.get().face_idx = fo.get_idx();
        f.get().last_idx = h.get_idx();
        fo.get().last_idx = ho.get_idx();

        link(hn.get_idx(), h.get_idx());
        link(h.get_idx(), hop.get_idx());
        link(hop.get_idx(), hn.get_idx());

        link(hon.get_idx(), ho.get_idx());
        link(ho.get_idx(), hp.get_idx());
        link(hp.get_idx(), hon.get_idx());

        h.get().vert_idx = von.get_idx();
        ho.get().vert_idx = vn.get_idx();

        v.get().out_idx = hn.get_idx();
        vo.get().out_idx = hon.get_idx();

        ensure_boundary_consistency(v.get_idx());
        ensure_boundary_consistency(vo.get_idx());

        return true;
    }
    void Manifold::cleanup()
    {
        //1. store new location of all entities in seperate arrays for future lookups
        vector<VertexIndex> vert_idx(vertices.full_size());
        vector<FaceIndex> face_idx(faces.full_size());
        vector<HalfEdgeIndex> he_idx(halfedges.full_size());

        IndexType k = 0;
        for(VertexIndex v = 0; v < vertices.full_size(); ++v){
            if(vertices.in_use(v)){
                vert_idx[v] = k;
                ++k;
            }
        }
        k = 0;
        for(FaceIndex f = 0; f < faces.full_size(); ++f){
            if(faces.in_use(f)){
                face_idx[f] = k;
                ++k;
            }
        }
        k = 0;
        for(HalfEdgeIndex h = 0; h < halfedges.full_size(); ++h){
            if(halfedges.in_use(h)){
                he_idx[h] = k;
                ++k;
            }
        }
        //2. update the remainding attributes of all entities with the indices stored in the touched attributes of their known entities
        for(VertexIndex v = 0; v < vertices.full_size(); ++v){
            if(vertices.in_use(v))
                vertices[v].out_idx = he_idx[vertices[v].out_idx];
        }
        for(FaceIndex f = 0; f < faces.full_size(); ++f){
            if(faces.in_use(f))
                faces[f].last_idx = he_idx[faces[f].last_idx];
        }
        for(HalfEdgeIndex h = 0; h < halfedges.full_size(); ++h){ 
            if(halfedges.in_use(h)){
                        if(halfedges[h].face_idx != NULL_FIDX)
                                halfedges[h].face_idx = face_idx[halfedges[h].face_idx];
                halfedges[h].next_idx = he_idx[halfedges[h].next_idx];
                halfedges[h].prev_idx = he_idx[halfedges[h].prev_idx];
                halfedges[h].opp_idx = he_idx[halfedges[h].opp_idx];
                halfedges[h].vert_idx = vert_idx[halfedges[h].vert_idx];
            }
        }
        vertices.cleanup();
        halfedges.cleanup();
        faces.cleanup();
    }
    void Manifold::enumerate_vertices()
    {
        IndexType k = 0;
        for(VertexHandle v = vhandles_begin(); v != vhandles_end(); ++v, ++k){
            v.get().touched = k;
        }
    }
    void Manifold::enumerate_halfedges()
    {
        IndexType k = 0;
        for(HalfEdgeHandle h = hhandles_begin(); h != hhandles_end(); ++h, ++k){
            h.get().touched = k;
        }
    }
    void Manifold::enumerate_faces()
    {
        IndexType k = 0;
        for(FaceHandle f = fhandles_begin(); f != fhandles_end(); ++f, ++k){
            f.get().touched = k;
        }
    }
    void Manifold::clear()
    {
        vertices.clear();
        halfedges.clear();
        faces.clear();
    }

    VertexHandle Manifold::vhandles_begin(){                return VertexHandle(*this, vertices.index_begin());}  
    HalfEdgeHandle Manifold::hhandles_begin(){              return HalfEdgeHandle(*this, halfedges.index_begin()); }
    FaceHandle Manifold::fhandles_begin(){                  return FaceHandle(*this, faces.index_begin());}
    ConstVertexHandle Manifold::vhandles_begin()const{      return ConstVertexHandle(*this, vertices.index_begin());  }
    ConstHalfEdgeHandle Manifold::hhandles_begin() const{   return ConstHalfEdgeHandle(*this, halfedges.index_begin()); }
    ConstFaceHandle Manifold::fhandles_begin() const{       return ConstFaceHandle(*this, faces.index_begin()); }

    VertexHandle Manifold::vhandles_end(){                  return VertexHandle(*this, vertices.index_end()); }
    HalfEdgeHandle Manifold::hhandles_end(){                return HalfEdgeHandle(*this, halfedges.index_end()); }
    FaceHandle Manifold::fhandles_end(){                    return FaceHandle(*this, faces.index_end()); }
    ConstVertexHandle Manifold::vhandles_end() const{       return ConstVertexHandle(*this, vertices.index_end()); }
    ConstHalfEdgeHandle Manifold::hhandles_end()const{      return ConstHalfEdgeHandle(*this, halfedges.index_end()); }
    ConstFaceHandle Manifold::fhandles_end()const{          return ConstFaceHandle(*this, faces.index_end()); }

    /**********************************************
    * Private functions
    **********************************************/

    template<typename float_type, typename int_type>
    void Manifold::build_template(  int_type no_vertices,
                                    const float_type* vertvec,
                                    int_type no_faces,
                                    const int_type* facevec,
                                    const int_type* indices,
                                    const int* touch)
    {
        // clean up the manifold for the build
        clear();

        // allocate space
        vertices.resize(no_vertices);
        faces.resize(no_faces);

        // create vertices corresponding to positions stored in vertvec
        for(int_type i = 0; i < no_vertices; ++i){
            Vec3f v;
            v[0] = vertvec[i*3];
            v[1] = vertvec[i*3+1];
            v[2] = vertvec[i*3+2];              
            vertices[i].pos = v;
        }

        // update touched values if the optional touch array was passed to build
        if(touch){
            for(int_type i = 0; i < no_vertices; ++i){
                vertices[i].touched = touch[i];
            }
        }

        //map over the created edges - needed to preserve edge uniqueness
        typedef map<EdgeKey, Edge> EdgeMap;
        EdgeMap edge_map;

        // counter that jumps between faces in indices vector
        int_type n  = 0;

        // create faces correspponding to faces stored in facevec
        for(int_type i = 0; i < no_faces; ++i){
            //amount of vertices in current face
            int_type N = facevec[i];
            vector<HalfEdgeIndex> fh;

            //each face indice results corresponds to 1 edge, 2 halfedges
            for(int_type j = 0; j < N; ++j){
                // each iteration uses two indices from the face
                VertexIndex v_idx0 = indices[j+n];
                VertexIndex v_idx1 = indices[(j+1)%N+n];

                // ensure indice integrity
                assert(v_idx0 < no_vertices);
                assert(v_idx1 < no_vertices);

                // create key and search map for edge
                EdgeKey ek(v_idx0, v_idx1);
                EdgeMap::iterator em_iter = edge_map.find(ek);

                // current edge has not been created
                if(em_iter == edge_map.end()){ 
                    // create edge for map
                    Edge e;
                    e.h0 = add_halfedge(HalfEdge());
                    e.h1 = add_halfedge(HalfEdge());
                    e.count = 1;

                    // glue operation: 1 edge = 2 glued halfedges
                    glue(e.h0, e.h1);

                    // update vertices with their outgoing halfedges
                    vertices[v_idx0].out_idx = e.h0;
                    vertices[v_idx1].out_idx = e.h1;

                    // update halfedges with the vertices they point to
                    halfedges[e.h0].vert_idx = v_idx1;
                    halfedges[e.h1].vert_idx = v_idx0;

                    // update map
                    edge_map[ek] = e;

                    // update vector of halfedges belonging to current face
                    fh.push_back(e.h0);
                }
                else{
                    // update current face with halfedge from edge
                    fh.push_back(em_iter->second.h1);

                    // asserting that a halfedge is visited exactly twice;
                    // once for each face on either side of the edge. 
                    em_iter->second.count++;
                    assert(em_iter->second.count == 2);
                }
            }
        
            for(int_type j = 0; j < N; ++j){
                // update halfedge with face
                halfedges[fh[j]].face_idx = i;

                // link operation: link two halfedges in the same face
                link(fh[j], fh[(j+1)%N]);
            }
            //update face with the first halfedge created
            faces[i].last_idx = fh[0];

            // step to first indice of next face
            n += N;
        }

        // test for unused vertices
        for(int_type i = 0; i < no_vertices; ++i){
            if(vertices[i].out_idx == NULL_HIDX)
                vertices.remove(i);
        }

        // boundary check while avoiding unused vertices
        for(int_type i = 0; i < no_vertices; ++i){
            if(vertices[i].out_idx != NULL_HIDX)
                ensure_boundary_consistency(i);
        }     
    }
    void Manifold::link(HalfEdgeIndex h0, HalfEdgeIndex h1)
    {
        halfedges[h0].next_idx = h1;
        halfedges[h1].prev_idx = h0;
    }
    void Manifold::glue(HalfEdgeIndex h0, HalfEdgeIndex h1)
    {
        halfedges[h0].opp_idx = h1;
        halfedges[h1].opp_idx = h0;
    }
    void Manifold::ensure_boundary_consistency(VertexIndex v)
    {
        // boundary consistency check by performing two vertex circulations
        HalfEdgeIndex h = vertices[v].out_idx;
        HalfEdgeIndex last = h;

        // step 1: circle through edges pointing away from vertex until reaching a null face
        while(halfedges[h].face_idx != NULL_FIDX){
            h = halfedges[halfedges[h].prev_idx].opp_idx;
            if(h == last) 
                return;
        }
        // null face encountered, we update our vertex with half edge index and prepare for step 2
        vertices[v].out_idx = h;
        HalfEdgeIndex ho = halfedges[h].opp_idx;
        
        // step 2: circle through edges pointing towards vertex until reaching a null face
        while(halfedges[ho].face_idx != NULL_FIDX){
            ho = halfedges[halfedges[ho].next_idx].opp_idx;
        }
        // null face encountered again, we update our edge with vertex index
        halfedges[ho].vert_idx = v;

        // remaining step is to make the in and out going edges link to each other
        link(ho, h);
    }
    void Manifold::remove_face_if_degenerate(HalfEdgeIndex h)
    {
        // face is degenerate if there is only two halfedges in face loop
        if(halfedges[halfedges[h].next_idx].next_idx == h)
        {
            HalfEdgeIndex hn = halfedges[h].next_idx;
            HalfEdgeIndex ho = halfedges[h].opp_idx;
            HalfEdgeIndex hno = halfedges[hn].opp_idx;

            // glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
            glue(ho, hno);

            VertexIndex hv = halfedges[h].vert_idx;
            VertexIndex hnv = halfedges[hn].vert_idx;

            // make v and vn point to opposite and next opposite halfedge, obsoleting h and hn from loop
            vertices[hnv].out_idx = hno;
            vertices[hv].out_idx = ho;

            FaceIndex f = halfedges[h].face_idx;

            // if face owning h is valid, remove face
            if(f != NULL_FIDX)
                remove_face(f);
            // remove the two invalid halfedges and the invalid face loop
            remove_halfedge(h);
            remove_halfedge(hn);

            // verify that v and vn is sane after removing the degenerate face
            ensure_boundary_consistency(hv);
            ensure_boundary_consistency(hnv);
        }
    }

    /***************************************************
    * Namespace functions
    ***************************************************/
    bool valid(const Manifold& m)
    {
        // Verify components of halfedges
        for(ConstHalfEdgeHandle h = m.hhandles_begin(); h != m.hhandles_end(); ++h){
                /* NOTE: Should not return false for boundary edges - thus commented out
                if(halfedges[h].face_idx == NULL_FIDX){
                    cout << "Halfedge [" << h << "] lacks face" << endl;
                   return false;
                }
                */
            if(h.vert() == NULL_CVHANDLE){
                cout << "Halfedge [" << h.get_idx() << "] lacks vert" << endl;
                return false;
            }
            if(h.next() == NULL_CHHANDLE){
                cout << "Halfedge [" << h.get_idx() << "] lacks next" << endl;
                return false;
            }
            if(h.prev() == NULL_CHHANDLE){
                cout << "Halfedge [" << h.get_idx() << "] lacks prev" << endl;
                return false;
            }
            if(h.opp() == NULL_CHHANDLE){
                cout << "Halfedge [" << h.get_idx() << "] lacks opp" << endl;
                return false;
            }

        }
        // Verify components of vertices
        for(ConstVertexHandle v = m.vhandles_begin(); v != m.vhandles_end(); ++v){
            vector<VertexIndex> link;

            // circulate the halfedges of vertex
            for(ConstVertexCirculator vc(v); !vc.end(); ++vc){
                // test halfedges around v
                if(vc.halfedge() == NULL_CHHANDLE){
                    cout << "Vertex [" << v.get_idx() << "] circulation produced null halfedge at step " << vc.no_steps() << endl;
                    return false;
                }
                VertexIndex ring_v = vc.halfedge().vert().get_idx();

                // test one-ring for multiple occurences of vertex
                if(find(link.begin(), link.end(), ring_v) != link.end()){
                    cout << "Vertex [" << ring_v << "] appears two times in the one-ring of vertex [" << v.get_idx() << "]" << endl;
                    return false;
                }
                link.push_back(ring_v);

                // test for infinite loop around vertex
                if(vc.no_steps() > m.no_vertices()){
                    cout << "Vertex [" << v.get_idx() << "] loop contains more vertices than manifold" << endl;
                    return false;
                }
            }

            // test one_ring size for boundary consistency
            if(link.size() == 2){
                ConstHalfEdgeHandle he = v.out();
                ConstHalfEdgeHandle opp = he.opp();

                if(he.face() != NULL_CFHANDLE && opp.face() != NULL_CFHANDLE)
                    cout << "Vertex [" << v .get_idx() << "] contains only two or less incident edges" << endl;
                else
                    cout << "Boundary vertex [" << v.get_idx() << "]  contains only two or less incident edges" << endl;
            }
            assert(link.size() >= 2);
        }
        // verify components of faces
        for(ConstFaceHandle f = m.fhandles_begin(); f != m.fhandles_end(); ++f){
            // count edges on face
            ConstFaceCirculator fc(f);
            for(; f != m.fhandles_end(); ++f){ 
                // test that all halfedges in faces bind properly to their face     
                if(fc.next_halfedge().face() != f){
                    cout << "Face [" << f.get_idx() << "] is inconsistent, halfedge [" << fc.next_halfedge().get_idx() << "] is not bound to face" << endl;
                    return false;
                }
            }
            // test faces for valid geometrical properties
            if(fc.no_steps() < 3){
                cout << "Face [" << f.get_idx() << "] contains only " << fc.no_steps() << " edge(s)" << endl;
                return false;
            }
            // test for infinite loop around face
            if(fc.no_steps() > m.no_halfedges() * 0.5f){
                cout << "Face [" << f.get_idx() << "] loop contains more halfedges than manifold" << endl;
                return false;
            }
        }
        return true;
    }
}