Rev 550 | Rev 553 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed
/* ----------------------------------------------------------------------- *
* This file is part of GEL, www.imm.dtu.dk/GEL
* Copyright (C) the authors (see AUTHORS.txt) and DTU Informatics
*
* Principal authors:
* Christian Thode Larsen (thode2d@gmail.com)
* J. Andreas Baerentzen (jab@imm.dtu.dk)
*
* See LICENSE.txt for licensing information
* ----------------------------------------------------------------------- */
#include "Manifold.h"
#include <iostream>
#include <vector>
#include <map>
#include <iterator>
#include <CGLA/Vec2f.h>
#include <CGLA/Vec3f.h>
#include <CGLA/Mat3x3f.h>
#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;
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::collapse_edge(HalfEdgeID h, bool avg_vertices)
{
HalfEdgeID ho = ck.opp(h);
VertexID hv = ck.vert(h);
VertexID hov = ck.vert(ho);
HalfEdgeID hn = ck.next(h);
HalfEdgeID hp = ck.prev(h);
HalfEdgeID hon = ck.next(ho);
HalfEdgeID hop = ck.prev(ho);
FaceID f = ck.face(h);
FaceID fo = ck.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 = ck.out(hov);
HalfEdgeID last = he;
do {
assert(ck.vert(ck.opp(he)) == hov);
ck.set_vert(ck.opp(he), hv);
he = ck.next(ck.opp(he));
}
while(he != last);
ck.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(ck.face(h) != InvalidFaceID)
ck.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(ck.face(ho) != InvalidFaceID)
ck.set_last(fo, hon);
// remove the obsolete entities
ck.remove_vertex(hov);
ck.remove_halfedge(h);
ck.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)
{
assert(!connected(*this, v0, v1));
HalfEdgeID he = ck.last(f);
HalfEdgeID last = he;
int steps = 0;
while(he != last || steps == 0){
++steps;
he = ck.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 = ck.out(v0);
for(HalfEdgeWalker w = halfedgewalker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){
if(w.face() == f){
h0 = w.halfedge();
break;
}
}
assert(ck.face(h0) != InvalidFaceID);
assert(ck.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(ck.vert(h) != v1){
h = ck.next(h);
}
assert(h != h0);
// create a new halfedge ha which connects v1 and v0 closing the first loop.
HalfEdgeID h1 = ck.next(h);
HalfEdgeID ha = ck.add_halfedge();
link(h, ha);
link(ha, h0);
ck.set_face(ha, f);
ck.set_vert(ha, v0);
ck.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 = ck.add_face();
while(ck.vert(h) != v0){
ck.set_face(h, f2);
h = ck.next(h);
}
ck.set_face(h, f2);
assert(h != h1);
// create a new halfedge hb to connect v0 and v1.
HalfEdgeID hb = ck.add_halfedge();
link(h, hb);
link(hb, h1);
ck.set_face(hb, f2);
ck.set_vert(hb, v1);
ck.set_last(f2, hb);
// complete the operation by gluing the two new halfedges
glue(ha, hb);
// assert sanity of operation
assert(ck.next(ck.opp(ck.prev(h1))) == h0);
assert(ck.next(ck.opp(ck.prev(h0))) == h1);
assert(ck.face(ck.next(hb)) == f2);
assert(ck.face(ck.next(ck.next(hb))) == f2);
assert(ck.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
Vec3f p(0.0f);
HalfEdgeID last = ck.last(f);
HalfEdgeID he = last;
int steps = 0;
while(he != last || steps == 0){
p += positions[ck.vert(he)];
++steps;
he = ck.next(he);
}
p /= steps;
VertexID v = ck.add_vertex();
positions[v] = p;
//circulate the face, create halfedges and connect to vertex
vector<HalfEdgeID> eout;
vector<HalfEdgeID> ein;
last = ck.last(f);
he = last;
do{
HalfEdgeID hn = ck.next(he);
HalfEdgeID ho = ck.add_halfedge();
HalfEdgeID hi = ck.add_halfedge();
FaceID fn = ck.add_face();
ck.set_face(hi, fn);
ck.set_vert(hi, v);
ck.set_face(ho, fn);
ck.set_vert(ho, ck.vert(ck.opp(he)));
ck.set_face(he, fn);
ck.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);
// glue new halfedges together
size_t N = eout.size();
for(size_t i = 0; i < N; ++i){
glue(ein[i], eout[(i+1)%N]);
}
ck.set_out(v, eout[0]);
//remove the now replaced face
ck.remove_face(f);
return v;
}
VertexID Manifold::split_edge(HalfEdgeID h)
{
HalfEdgeID ho = ck.opp(h);
VertexID v = ck.vert(h);
VertexID vo = ck.vert(ho);
//create the new vertex with middle of edge as position and update connectivity
VertexID vn = ck.add_vertex();
positions[vn] = .5f * (positions[v] + positions[vo]);
ck.set_out(vn, h);
//create two necessary halfedges, and update connectivity
HalfEdgeID hn = ck.add_halfedge();
HalfEdgeID hno = ck.add_halfedge();
ck.set_out(vo, hn);
ck.set_out(v, ho);
glue(hn, hno);
link(ck.prev(h), hn);
link(hn, h);
ck.set_vert(hn, vn);
link(hno, ck.next(ho));
link(ho, hno);
ck.set_vert(hno, vo);
ck.set_vert(ho, vn);
// update faces in case of non boundary edge
if(ck.face(h) != InvalidFaceID)
ck.set_last(ck.face(h), hn);
if(ck.face(ho) != InvalidFaceID)
ck.set_last(ck.face(ho), ho);
// update new edge with faces from existing edge
ck.set_face(hn, ck.face(h));
ck.set_face(hno, ck.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;
}
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 halfedgewalker 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;
HalfEdgeWalker hew = halfedgewalker(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(HalfEdgeWalker hew2 = halfedgewalker(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(int i=0;i<loop.size();++i)
if(halfedgewalker(loop[i]).opp().halfedge() == loop[(i+1)%loop.size()])
{
VertexID vid = halfedgewalker(loop[i]).vertex();
if(vid != v)
garbage_vertices.push_back(halfedgewalker(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(int 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(int 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(int i=0;i<garbage_vertices.size(); ++i)
ck.remove_vertex(garbage_vertices[i]);
for(int i=0;i<garbage_faces.size(); ++i)
ck.remove_face(garbage_faces[i]);
for(int i=0;i<garbage_halfedges.size(); ++i)
ck.remove_halfedge(garbage_halfedges[i]);
if(!vertex_is_boundary)
ck.remove_vertex(v);
// Create a new face for the merged one ring and link up all the halfedges
// in the loop
FaceID f = ck.add_face();
ck.set_last(f,loop[0]);
for(int i=0;i<loop.size(); ++i)
{
ck.set_face(loop[i], f);
HalfEdgeWalker hw = halfedgewalker(loop[i]);
ck.set_out(hw.vertex(),hw.opp().halfedge());
link(loop[i],loop[(i+1)%loop.size()]);
assert(hw.vertex() == halfedgewalker(loop[(i+1)%loop.size()]).opp().vertex());
}
// Finally, we ensure boundary consitency for all vertices in the loop.
for(int i=0;i<loop.size(); ++i)
{
HalfEdgeWalker hw = halfedgewalker(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(ck.face(h) == f);
HalfEdgeID ho = ck.opp(h);
FaceID fo = ck.face(ho);
HalfEdgeID hn = ck.next(h);
HalfEdgeID hon = ck.next(ho);
if(fo == f)
return false;
//boundary case
if(ck.vert(hn) == ck.vert(hon))
return false;
HalfEdgeID hp = ck.prev(h);
HalfEdgeID hop = ck.prev(ho);
VertexID v = ck.vert(h);
VertexID vo = ck.vert(ho);
if(valency(*this, v) < 3 || valency(*this, vo) < 3)
return false;
link(hop, hn);
link(hp, hon);
ck.set_out(vo, hon);
ck.set_out(v, hn);
ck.set_last(f, hn);
HalfEdgeID hx = hon;
assert(ck.face(hx) == fo);
while(ck.face(hx) != f){
ck.set_face(hx, f);
hx = ck.next(hx);
}
ensure_boundary_consistency(v);
ensure_boundary_consistency(vo);
ck.remove_halfedge(h);
ck.remove_halfedge(ho);
ck.remove_face(fo);
return true;
}
void Manifold::close_hole(HalfEdgeID h)
{
// invalid face is a hole
if(ck.face(h) == InvalidFaceID){
FaceID f = ck.add_face();
ck.set_last(f, h);
do{
ck.set_face(h, f);
h = ck.next(h);
}
while(ck.face(h) != f);
}
}
void Manifold::flip_edge(HalfEdgeID h)
{
HalfEdgeID hn = ck.next(h);
HalfEdgeID hp = ck.prev(h);
HalfEdgeID ho = ck.opp(h);
HalfEdgeID hon = ck.next(ho);
HalfEdgeID hop = ck.prev(ho);
FaceID hf = ck.face(h);
FaceID hof = ck.face(ho);
// update face connectivity of halfedges changing face
ck.set_face(hop, hf);
ck.set_face(hp, hof);
// connectivity of faces with halfedges of flipped edge
ck.set_last(hf, h);
ck.set_last(hof, ho);
// connectivity of halfedges of first face after flip
link(hn, h);
link(h, hop);
link(hop, hn);
// connectivity of halfedges of second face after flip
link(hon, ho);
link(ho, hp);
link(hp, hon);
VertexID hv = ck.vert(h);
VertexID hnv = ck.vert(hn);
VertexID hov = ck.vert(ho);
VertexID honv = ck.vert(hon);
// connectivity of flipped edge and corresponding vertices
ck.set_vert(h, honv);
ck.set_vert(ho, hnv);
ck.set_out(hv, hn);
ck.set_out(hov, hon);
// if the flip occurs next to a boundary, ensure the boundary consistency
ensure_boundary_consistency(hv);
ensure_boundary_consistency(hov);
}
/**********************************************
* Private functions
**********************************************/
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,
const int_type* facevec,
const int_type* indices)
{
vector<VertexID> vids(no_vertices);
// create vertices corresponding to positions stored in vertvec
for(int i=0;i<no_vertices;++i)
{
const float* v0 = &vertvec[i*3];
pos(vids[i] = ck.add_vertex()) = Vec3f(v0[0], v0[1], v0[2]);
}
//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(size_type i = 0; i < no_faces; ++i){
//amount of vertices in current face
size_type N = facevec[i];
vector<HalfEdgeID> fh;
//each face indice results corresponds to 1 edge, 2 halfedges
for(size_type j = 0; j < N; ++j){
// ensure indice integrity
assert(static_cast<size_type>(indices[j + n]) < no_vertices);
assert(static_cast<size_type>(indices[(j + 1) % N + n]) < no_vertices);
// each iteration uses two indices from the face
VertexID v0 = vids[static_cast<size_type>(indices[j + n])];
VertexID v1 = vids[static_cast<size_type>(indices[(j + 1) % N + n])];
// create key and search map for edge
EdgeKey ek(v0, v1);
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 = ck.add_halfedge();
e.h1 = ck.add_halfedge();
e.count = 1;
// glue operation: 1 edge = 2 glued halfedges
glue(e.h0, e.h1);
// update vertices with their outgoing halfedges
ck.set_out(v0, e.h0);
ck.set_out(v1, e.h1);
// update halfedges with the vertices they point to
ck.set_vert(e.h0, v1);
ck.set_vert(e.h1, v0);
// 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);
}
}
FaceID fid = ck.add_face();
for(size_type j = 0; j < N; ++j){
// update halfedge with face
ck.set_face(fh[j], fid);
// link operation: link two halfedges in the same face
link(fh[j], fh[(j + 1) % N]);
}
//update face with the first halfedge created
ck.set_last(fid, fh[0]);
// step to first index of next face
n += N;
}
// test for unused vertices
for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
assert( (*v) != InvalidVertexID);
if(ck.out(*v) == InvalidHalfEdgeID)
ck.remove_vertex(*v);
}
// boundary check while avoiding unused vertices
for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
if((*v) != InvalidVertexID)
ensure_boundary_consistency(*v);
}
}
void Manifold::link(HalfEdgeID h0, HalfEdgeID h1)
{
ck.set_next(h0, h1);
ck.set_prev(h1, h0);
}
void Manifold::glue(HalfEdgeID h0, HalfEdgeID h1)
{
ck.set_opp(h0, h1);
ck.set_opp(h1, h0);
}
void Manifold::ensure_boundary_consistency(VertexID v)
{
// boundary consistency check by performing two vertex circulations
HalfEdgeID h = ck.out(v);
HalfEdgeID last = h;
// step 1: circle through edges pointing away from vertex until reaching a null face
while(ck.face(h) != InvalidFaceID){
h = ck.opp(ck.prev(h));
if(h == last) // We came full circle - vertex not boundary - return.
return;
}
// null face encountered, we update our vertex with half edge index and prepare for step 2
ck.set_out(v, h);
HalfEdgeID ho = ck.opp(h);
// step 2: circle through edges pointing towards vertex until reaching a null face
while(ck.face(ho) != InvalidFaceID){
ho = ck.opp(ck.next(ho));
}
// null face encountered again, we update our edge with vertex index
ck.set_vert(ho, 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(HalfEdgeID h)
{
// face is degenerate if there is only two halfedges in face loop
if(ck.next(ck.next(h)) == h)
{
HalfEdgeID hn = ck.next(h);
HalfEdgeID ho = ck.opp(h);
HalfEdgeID hno = ck.opp(hn);
VertexID hv = ck.vert(h);
VertexID hnv = ck.vert(hn);
FaceID f = ck.face(h);
assert(ho != hn);
assert(h != hno);
assert(hv != hnv);
// glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
glue(ho, hno);
// make v and vn point to opposite and next opposite halfedge, obsoleting h and hn from loop
ck.set_out(hnv, hno);
ck.set_out(hv, ho);
// if face owning h is valid, remove face
if(f != InvalidFaceID)
ck.remove_face(f);
// remove the two invalid halfedges and the invalid face loop
ck.remove_halfedge(h);
ck.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(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
HalfEdgeWalker j = m.halfedgewalker(*h);
if(j.vertex() == InvalidVertexID){
cout << "Halfedge lacks vert" << endl;
return false;
}
if(j.next().halfedge() == InvalidHalfEdgeID){
cout << "Halfedge lacks next" << endl;
return false;
}
if(j.prev().halfedge() == InvalidHalfEdgeID){
cout << "Halfedge lacks prev" << endl;
return false;
}
if(j.opp().halfedge() == InvalidHalfEdgeID){
cout << "Halfedge lacks opp" << endl;
return false;
}
}
// Verify components of vertices
for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
vector<VertexID> link;
// circulate the halfedges of vertex
for(HalfEdgeWalker j = m.halfedgewalker(*v); !j.full_circle(); j = j.circulate_vertex_cw()){
// test halfedges around v
if(j.halfedge() == InvalidHalfEdgeID){
cout << "Vertex circulation produced invalid halfedge" << endl;
return false;
}
VertexID ring_v = j.vertex();
// test one-ring for multiple occurences of vertex
if(find(link.begin(), link.end(), ring_v) != link.end()){
cout << "Vertex appears two times in one-ring of vertex" << endl;
return false;
}
link.push_back(ring_v);
// test for infinite loop around vertex
if(static_cast<size_t>(j.no_steps()) > m.active_vertices()){
cout << "Vertex loop contains more vertices than manifold" << endl;
return false;
}
}
// test one_ring size for boundary consistency
if(link.size() <= 2){
HalfEdgeWalker j = m.halfedgewalker(*v);
if(j.face() != InvalidFaceID && j.opp().face() != InvalidFaceID)
{
if(link.size()==1)
cout << "Vertex contains only a single incident edge" << endl;
//
// cout << "Vertex contains only " << link.size() <<" incident edges" << endl;
}
else
cout << "Boundary vertex contains only " << link.size() <<" incident edges" << endl;
}
}
// verify components of faces
for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
// count edges on face
HalfEdgeWalker j = m.halfedgewalker(*f);
for(; !j.full_circle(); j = j.circulate_face_cw()){
// test that all halfedges in faces bind properly to their face
if(j.face() != *f){
cout << "Face is inconsistent, halfedge is not bound to face" << endl;
return false;
}
}
// test faces for valid geometrical properties
if(j.no_steps() < 3){
cout << "Face contains less than 3 edges" << endl;
return false;
}
// test for infinite loop around face
if(j.no_steps() > m.active_halfedges() * 0.5f){
cout << "Face loop contains more halfedges than manifold" << endl;
return false;
}
}
return true;
}
void bbox(const Manifold& m, Vec3f& pmin, Vec3f& pmax)
{
VertexIDIterator v = m.vertices_begin();
pmin = pmax = m.pos(*v);
++v;
for(; v != m.vertices_end(); ++v){
pmin = v_min(m.pos(*v), pmin);
pmax = v_max(m.pos(*v), pmax);
}
}
void bsphere(const Manifold& m, Vec3f& c, float& r)
{
Vec3f p0, p7;
bbox(m, p0, p7);
Vec3f rad = (p7 - p0) * 0.5f;
c = p0 + rad;
r = rad.length();
}
bool precond_collapse_edge(const Manifold& m, HalfEdgeID h)
{
HalfEdgeWalker hew = m.halfedgewalker(h);
HalfEdgeID ho = hew.opp().halfedge();
VertexID v0 = hew.opp().vertex();
VertexID v1 = hew.vertex();
// get the one-ring of v0
vector<VertexID> link0;
vector<FaceID> faces0;
for(HalfEdgeWalker vj = m.halfedgewalker(h);
!vj.full_circle(); vj = vj.circulate_vertex_ccw()){
link0.push_back(vj.vertex());
if(vj.halfedge() != h)
faces0.push_back(vj.face());
}
assert(link0.size() > 1);
// get the one-ring of v1
vector<VertexID> link1;
vector<FaceID> faces1;
for(HalfEdgeWalker vj = m.halfedgewalker(ho);
!vj.full_circle(); vj = vj.circulate_vertex_ccw()){
link1.push_back(vj.vertex());
if(vj.halfedge() != ho)
faces1.push_back(vj.face());
}
assert(link1.size() > 1);
// 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
vector<VertexID> lisect;
std::back_insert_iterator<vector<VertexID> > lii(lisect);
set_intersection(link0.begin(), link0.end(),
link1.begin(), link1.end(),
lii);
// sort the faces of the two rings
sort(faces0.begin(), faces0.end());
sort(faces1.begin(), faces1.end());
// get the intersection of the shared faces from both rings
vector<FaceID> fisect;
std::back_insert_iterator<vector<FaceID> > fii(fisect);
set_intersection(faces0.begin(), faces0.end(),
faces1.begin(), faces1.end(),
fii);
if(fisect.size() > 0)
return false;
int k = 0;
// if the adjacent face is a triangle (see 2)
if(hew.next().next().next().halfedge() == h){
VertexID v = hew.next().vertex();
// valency test (see 5)
if(valency(m, v) < 3)
return false;
// remove the vertex shared by the two rings from the intersection
vector<VertexID>::iterator iter;
iter = find(lisect.begin(), lisect.end(), v);
assert(iter != lisect.end());
lisect.erase(iter);
++k;
}
// 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();
// valency test (see 5)
if(valency(m, v) < 3)
return false;
// remove the vertex shared by the two rings from the intersection
vector<VertexID>::iterator iter;
iter = find(lisect.begin(), lisect.end(), v);
assert(iter != lisect.end());
lisect.erase(iter);
++k;
}
// double edge test (see 3)
if(lisect.size() != 0)
return false;
// tetrahedon test (see 4)
if(k == 2 && (link0.size() + link1.size() == 6))
return false;
// test that we do not merge holes (see 6)
if(boundary(m, v0) && boundary(m, v1) && hew.face() != InvalidFaceID && hew.opp().face() != InvalidFaceID)
return false;
return true;
}
bool precond_flip_edge(const Manifold& m, HalfEdgeID h)
{
HalfEdgeWalker j = m.halfedgewalker(h);
FaceID hf = j.face();
HalfEdgeID ho = j.opp().halfedge();
FaceID hof = j.opp().face();
// boundary case
if(hf == InvalidFaceID || hof == InvalidFaceID)
return false;
// We can only flip an edge if both incident polygons are triangles.
if(no_edges(m, hf) != 3 || no_edges(m, hof) !=3)
return false;
// non boundary vertices with a valency of less than 4(less than 3 after operation) degenerates mesh.
VertexID hv = j.vertex();
VertexID hov = j.opp().vertex();
if((valency(m, hv) < 4 && !boundary(m, hv)) || (valency(m, hov) < 4 && !boundary(m, hov)))
return false;
// Disallow flip if vertices being connected already are.
VertexID hnv = j.next().vertex();
VertexID honv = j.opp().next().vertex();
if(connected(m, hnv, honv))
return false;
return true;
}
bool boundary(const Manifold& m, VertexID v)
{
HalfEdgeWalker j = m.halfedgewalker(v);
return boundary(m, j.halfedge());
}
int valency(const Manifold& m, VertexID v)
{
// perform full circulation to get valency
HalfEdgeWalker vj = m.halfedgewalker(v);
while(!vj.full_circle())
vj = vj.circulate_vertex_cw();
return vj.no_steps();
}
CGLA::Vec3f normal(const Manifold& m, VertexID v)
{
Vec3f p0 = m.pos(v);
vector<Vec3f> one_ring;
// run through outgoing edges, and store them normalized
for(HalfEdgeWalker vj = m.halfedgewalker(v); !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
Vec3f edge = p0 - m.pos(vj.vertex());
float l = length(edge);
if(l > 0.0f) one_ring.push_back(edge/l);
}
Vec3f n(0);
size_t N = one_ring.size();
// sum up the normals of each face surrounding the vertex
for(size_t i = 0; i < N; ++i){
Vec3f e0 = one_ring[i];
Vec3f e1 = one_ring[(i+1) % N];
Vec3f n_part = normalize(cross(e0, e1));
n += n_part * acos(max(-1.0f, min(1.0f, dot(e0, e1))));
}
// normalize and return the normal
float sqr_l = sqr_length(n);
if(sqr_l > 0.0f) return n / sqrt(sqr_l);
return n;
}
bool connected(const Manifold& m, VertexID v0, VertexID v1)
{
for(HalfEdgeWalker vj = m.halfedgewalker(v0); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
if(vj.vertex() == v1)
return true;
}
return false;
}
int no_edges(const Manifold& m, FaceID f)
{
// perform full circulation to get valency
HalfEdgeWalker w = m.halfedgewalker(f);
for(; !w.full_circle(); w = w.circulate_face_cw());
return w.no_steps();
}
CGLA::Vec3f normal(const Manifold& m, FaceID f)
{
vector<Vec3f> v;
int k=0;
for(HalfEdgeWalker w = m.halfedgewalker(f); !w.full_circle(); w = w.circulate_face_ccw(),++k)
v.push_back(m.pos(w.vertex()));
Vec3f norm(0);
for(int i=0;i<k;++i)
{
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]);
}
float l = sqr_length(norm);
if(l>0.0f)
norm /= sqrt(l);
return norm;
}
void Manifold::bridge_faces(FaceID f0, FaceID f1, const vector<pair<VertexID, VertexID> >& pairs)
{
// Let N be the number of vertex pairs to connect by edges
int N = pairs.size();
// 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);
for(int i=0;i<N; ++i)
{
new_halfedges[i] = ck.add_halfedge();
new_halfedges_opp[i] = ck.add_halfedge();
glue(new_halfedges[i], new_halfedges_opp[i]);
ck.set_vert(new_halfedges[i], pairs[i].second);
ck.set_vert(new_halfedges_opp[i], pairs[i].first);
}
// We now maintain some halfedge indices to right before
// and right after the point we are trying to connect on
// each of the two loops.
HalfEdgeWalker hw0 = halfedgewalker(f0);
while (hw0.vertex() != pairs[0].first)
hw0 = hw0.circulate_face_ccw();
HalfEdgeWalker hw1 = halfedgewalker(f1);
while (hw1.vertex() != pairs[0].second)
hw1 = hw1.circulate_face_ccw();
HalfEdgeID h0p = hw0.halfedge();
HalfEdgeID h1n = hw1.next().halfedge();
// loop over all connections and link the new halfedges with the old
// ones.
for(int i=0;i<N; ++i)
{
HalfEdgeID h0n = ck.next(h0p);
HalfEdgeID h1p = ck.prev(h1n);
link(h0p, new_halfedges[i]);
link(new_halfedges[i],h1n);
link(h1p, new_halfedges_opp[i]);
link(new_halfedges_opp[i],h0n);
h0p = new_halfedges_opp[i];
h1n = new_halfedges_opp[i];
while(ck.vert(h0p) != pairs[(i+1)%N].first)
h0p = ck.next(h0p);
while(ck.vert(ck.prev(h1n)) != pairs[(i+1)%N].second)
h1n = ck.prev(h1n);
}
// Create the faces and their connections
for(int i=0;i<N; ++i)
{
HalfEdgeID last = new_halfedges[i];
FaceID f = ck.add_face();
ck.set_last(f, last);
HalfEdgeID h = last;
do
{
ck.set_face(h, f);
h = ck.next(h);
} while(h != last);
}
// Delete the old faces that were bridged.
ck.remove_face(f0);
ck.remove_face(f1);
}
float area(const Manifold& m, FaceID f)
{
// Get all projected vertices
vector<Vec3f> vertices;
HalfEdgeWalker w = m.halfedgewalker(f);
for(; !w.full_circle(); w = w.circulate_face_ccw())
vertices.push_back(m.pos(w.vertex()));
float area = 0;
int k = w.no_steps();
Vec3f n = normal(m,f);
for(int i = 0; i < k; ++i)
{
Vec3f A = cross(vertices[i], vertices[(i+1 )%k]);
float s = sign(dot(n, A));
area += 0.5 * s * length(A);
}
return abs(area);
}
CGLA::Vec3f centre(const Manifold& m, FaceID f)
{
Vec3f c(0);
HalfEdgeWalker w = m.halfedgewalker(f);
for(; !w.full_circle(); w = w.circulate_face_cw()) c += m.pos(w.vertex());
c /= w.no_steps();
return c;
}
float perimeter(const Manifold& m, FaceID f)
{
float p=0.0;
HalfEdgeWalker w = m.halfedgewalker(f);
for(; !w.full_circle(); w = w.circulate_face_cw()) p += length(m, w.halfedge());
return p;
}
bool boundary(const Manifold& m, HalfEdgeID h)
{
HalfEdgeWalker w = m.halfedgewalker(h);
return w.face() == InvalidFaceID || w.opp().face() == InvalidFaceID;
}
float length(const Manifold& m, HalfEdgeID h)
{
HalfEdgeWalker w = m.halfedgewalker(h);
return (m.pos(w.vertex()) - m.pos(w.opp().vertex())).length();
}
}