Rev 535 | Rev 537 | 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);
// average the vertex positions
pos(hv) = avg_vertices * 0.5f * (pos(hov) - pos(hv)) + pos(hv);
HalfEdgeID hn = ck.next(h);
HalfEdgeID hp = ck.prev(h);
HalfEdgeID hon = ck.next(ho);
HalfEdgeID hop = ck.prev(ho);
// update all halfedges pointing to hov to point to hv, effectively removing hov from all loops
HalfEdgeID he = ck.out(hov);
HalfEdgeID last = he;
int steps = 0;
while(he != last || steps == 0){
ck.set_vert(ck.opp(he), hv);
++steps;
he = ck.next(ck.opp(he));
}
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(ck.face(h), 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(ck.face(ho), 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;
}
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);
//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);
link(hop, hn);
ck.set_out(v, hn);
link(hp, hon);
ck.set_out(vo, hn);
ck.set_last(f, hn);
HalfEdgeID hx = hon;
assert(ck.face(hx) == fo);
while(ck.face(h) != 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)
{
// clean up the manifold for the build
clear();
// allocate space
ck.resize_vertices(no_vertices);
positions.resize(no_vertices);
ck.resize_faces(no_faces);
vector<VertexID> vids;
vids.reserve(no_vertices);
copy(vertices_begin(), vertices_end(), back_inserter(vids));
vector<FaceID> fids;
fids.reserve(no_faces);
copy(faces_begin(), faces_end(), back_inserter(fids));
// create vertices corresponding to positions stored in vertvec
int i = 0;
for(VertexID v = ck.vertices_begin(); v != ck.vertices_end(); v = ck.vertices_next(v), ++i){
Vec3f p;
p[0] = vertvec[i*3];
p[1] = vertvec[i*3+1];
p[2] = vertvec[i*3+2];
positions[v] = p;
}
//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);
}
}
for(size_type j = 0; j < N; ++j){
// update halfedge with face
ck.set_face(fh[j], fids[i]);
// 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(fids[i], fh[0]);
// step to first indice 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)
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);
// glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
glue(ho, hno);
VertexID hv = ck.vert(h);
VertexID hnv = ck.vert(hn);
// 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);
FaceID f = ck.face(h);
// 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)
cout << "Vertex contains only two or less incident edges" << endl;
else
cout << "Boundary vertex contains only two or less incident edges" << endl;
}
assert(link.size() >= 2);
}
// 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 j = m.halfedgewalker(h);
// get the two vertices of the edge being tested
VertexID v0 = j.opp().vertex();
VertexID v1 = j.vertex();
// get the one-ring of v0
vector<VertexID> link0;
for(HalfEdgeWalker vj = m.halfedgewalker(v0); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
link0.push_back(vj.vertex());
}
assert(link0.size() > 1);
// get the one-ring of v1
vector<VertexID> link1;
for(HalfEdgeWalker vj = m.halfedgewalker(v1); !vj.full_circle(); vj = vj.circulate_vertex_cw()){
link1.push_back(vj.vertex());
}
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);
int k = 0;
// if the adjacent face is a triangle (see 2)
if(j.next().next().next().halfedge() == h){
VertexID v = j.next().vertex();
// valency test (see 5)
if(valency(m, v) < 4)
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(j.opp().next().next().next().halfedge() == j.opp().halfedge()){
VertexID v = j.opp().next().vertex();
// valency test (see 5)
if(valency(m, v) < 4)
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) && j.face() != InvalidFaceID && j.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_cw()){
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 e1 = one_ring[i];
Vec3f e0 = 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> vertices;
HalfEdgeWalker w = m.halfedgewalker(f);
// Assuming that a face is valid, i.e. consisting of AT LEAST 3 vertices...
for(int i = 0; i < 3; ++i, w = w.circulate_face_cw())
vertices.push_back(m.pos(w.vertex()));
Vec3f n = cross(vertices[0] - vertices[1], vertices[2] - vertices[1]);
float l = sqr_length(n);
if(l > 0.0f) return n/sqrt(l);
return Vec3f(0);
}
float area(const Manifold& m, FaceID f)
{
HalfEdgeWalker w = m.halfedgewalker(f);
int k = 0;
// M is a matrix that projects a vector onto the orthogonal
// complement of the face normal
Vec3f n = normal(m, f);
if(sqr_length(n)<0.1) // actually this number is arbitrary - the normal is unit length or zero length
return 0; // threshold .1 just to avoid floating point issues.
Vec3f a,b;
orthogonal(n, a, b);
Mat3x3f M(a, b, n);
// Get all projected vertices
vector<Vec2f> v;
for(k = 0; !w.full_circle(); w = w.circulate_face_cw(), ++k){
Vec3f p = M * m.pos(w.vertex());
v.push_back(Vec2f(p[0], p[1]));
}
float area = 0;
for(int i = 0; i < k; ++i) area += 0.5 * cross(v[i], v[(i+1 )% k]);
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;
}
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();
}
}