commit Pierre files, maybe temporarily

This commit is contained in:
Laurent Rineau 2007-11-14 11:02:04 +00:00
parent a61decfce6
commit 5aab4f8ab7
3 changed files with 2107 additions and 0 deletions

View File

@ -0,0 +1,823 @@
/***************************************************************************
builder.h
----------------------------------------------------------------------------
begin : june 2003
copyright : (C) 2003 by Pierre Alliez - INRIA
email : pierre.alliez@sophia.inria.fr
Enriched incremental builder for CGAL Polyhedral Surfaces.
***************************************************************************/
#ifndef Enriched_polyhedron_incremental_builder_3_H
#define Enriched_polyhedron_incremental_builder_3_H 1
#include <CGAL/basic.h>
// MS Visual C++ 6.0 does not work with the new design.
#if defined( _MSC_VER) && (_MSC_VER <= 1200)
#ifndef CGAL_USE_POLYHEDRON_DESIGN_TWO
#define CGAL_USE_POLYHEDRON_DESIGN_ONE 1
#endif
#endif
#ifdef CGAL_USE_POLYHEDRON_DESIGN_ONE
#include <CGAL/Polyhedron_old/Polyhedron_incremental_builder_3.h>
#else // CGAL_USE_POLYHEDRON_DESIGN_ONE //
#define CGAL_USE_POLYHEDRON_DESIGN_TWO 1
#include <CGAL/Random_access_adaptor.h>
#include <CGAL/HalfedgeDS_decorator.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/IO/Verbose_ostream.h>
#include <vector>
#include <cstddef>
CGAL_BEGIN_NAMESPACE
template < class HalfedgeDS_>
class Enriched_polyhedron_incremental_builder_3 {
public:
typedef HalfedgeDS_ HDS; // internal
typedef HalfedgeDS_ HalfedgeDS;
typedef typename HDS::Vertex Vertex;
typedef typename HDS::Halfedge Halfedge;
typedef typename HDS::Face Face;
typedef typename HDS::Vertex_handle Vertex_handle;
typedef typename HDS::Halfedge_handle Halfedge_handle;
typedef typename HDS::Face_handle Face_handle;
typedef typename HDS::Face_handle Facet_handle;
typedef typename Vertex::Base VBase;
typedef typename Halfedge::Base HBase;
typedef typename Vertex::Point Point_3;
typedef typename HDS::size_type size_type;
protected:
typedef typename HDS::Supports_vertex_halfedge Supports_vertex_halfedge;
typedef typename HDS::Supports_removal Supports_removal;
typedef typename HDS::Vertex_iterator Vertex_iterator;
typedef typename HDS::Halfedge_iterator Halfedge_iterator;
typedef Random_access_adaptor<Vertex_iterator> Random_access_index;
bool m_error;
bool m_verbose;
HDS& hds;
size_type rollback_v;
size_type rollback_f;
size_type rollback_h;
size_type new_vertices;
size_type new_faces;
size_type new_halfedges;
Face_handle current_face;
Random_access_index index_to_vertex_map;
std::vector< Halfedge_handle> vertex_to_edge_map;
Halfedge_handle g1; // first halfedge, 0 denotes none.
Halfedge_handle gprime;
Halfedge_handle h1; // current halfedge
size_type w1; // first vertex.
size_type w2; // second vertex.
size_type v1; // current vertex
bool first_vertex;
bool last_vertex;
CGAL_assertion_code( int check_protocoll;) // use to check protocoll.
// states for checking: 0 = created, 1 = constructing, 2 = make face.
// Implement the vertex_to_edge_map either with an array or
// the halfedge pointer in the vertices (if supported).
// ----------------------------------------------------
void initialize_vertex_to_edge_map( size_type , bool, Tag_true) {}
void initialize_vertex_to_edge_map( size_type n, bool mode, Tag_false){
vertex_to_edge_map = std::vector< Halfedge_handle>();
vertex_to_edge_map.reserve(n);
if ( mode) {
// go through all halfedges and keep a halfedge for each
// vertex found in a hashmap.
typedef Unique_hash_map< Vertex_iterator, Halfedge_handle> V_map;
Halfedge_handle hh;
V_map v_map( hh, hds.size_of_vertices());
for ( Halfedge_iterator hi = hds.halfedges_begin();
hi != hds.halfedges_end();
++hi) {
v_map[ hi->vertex()] = hi;
}
size_type i = 0;
for ( Vertex_iterator vi = hds.vertices_begin();
vi != hds.vertices_end();
++vi) {
set_vertex_to_edge_map( i, v_map[ index_to_vertex_map[i]]);
++i;
}
}
}
void initialize_vertex_to_edge_map( size_type n, bool mode) {
initialize_vertex_to_edge_map(n, mode, Supports_vertex_halfedge());
}
void push_back_vertex_to_edge_map( Halfedge_handle , Tag_true) {}
void push_back_vertex_to_edge_map( Halfedge_handle h, Tag_false) {
vertex_to_edge_map.push_back(h);
}
void push_back_vertex_to_edge_map( Halfedge_handle h) {
push_back_vertex_to_edge_map( h, Supports_vertex_halfedge());
}
Halfedge_handle get_vertex_to_edge_map( int i, Tag_true) {
// Use the halfedge pointer within the vertex.
return index_to_vertex_map[i]->halfedge();
}
Halfedge_handle get_vertex_to_edge_map( int i, Tag_false) {
// Use the self-managed array vertex_to_edge_map.
return vertex_to_edge_map[i];
}
Halfedge_handle get_vertex_to_edge_map( int i) {
return get_vertex_to_edge_map( i, Supports_vertex_halfedge());
}
void set_vertex_to_edge_map( int i, Halfedge_handle h, Tag_true) {
// Use the halfedge pointer within the vertex.
index_to_vertex_map[i]->VBase::set_halfedge(h);
}
void set_vertex_to_edge_map( int i, Halfedge_handle h, Tag_false) {
// Use the self-managed array vertex_to_edge_map.
vertex_to_edge_map[i] = h;
}
void set_vertex_to_edge_map( int i, Halfedge_handle h) {
set_vertex_to_edge_map( i, h, Supports_vertex_halfedge());
}
// An Incremental Builder for Polyhedral Surfaces
// ----------------------------------------------
// DEFINITION
//
// Enriched_polyhedron_incremental_builder_3<HDS> is an auxiliary class that
// supports the incremental construction of polyhedral surfaces. This is
// for example convinient when constructing polyhedral surfaces from
// files. The incremental construction starts with a list of all point
// coordinates and concludes with a list of all facet polygons. Edges are
// not explicitly specified. They are derived from the incidence
// information provided from the facet polygons. These are given as a
// sequence of vertex indices. The correct protocol of method calls to
// build a polyhedral surface can be stated as regular expression:
//
// `begin_surface (add_vertex | (begin_facet add_vertex_to_facet*
// end_facet))* end_surface '
//
// PARAMETERS
//
// `HDS' is the halfedge data structure used to represent the
// polyhedral surface that is to be constructed.
//
// CREATION
public:
bool error() const { return m_error; }
Enriched_polyhedron_incremental_builder_3( HDS& h, bool verbose = false)
// stores a reference to the halfedge data structure `h' in the
// internal state. The previous polyhedral surface in `h'
// remains unchanged. The incremental builder adds the new
// polyhedral surface to the old one.
: m_error( false), m_verbose( verbose), hds(h) {
CGAL_assertion_code(check_protocoll = 0;)
}
~Enriched_polyhedron_incremental_builder_3() {
CGAL_assertion( check_protocoll == 0);
}
// OPERATIONS
enum { CGAL_RELATIVE = 0, CGAL_ABSOLUTE = 1};
void begin_surface( std::size_t v, std::size_t f, std::size_t h = 0,
int mode = CGAL_RELATIVE);
// starts the construction. v is the number of new
// vertices to expect, f the number of new facets, and h the number of
// new halfedges. If h is unspecified (`== 0') it is estimated using
// Euler equations (plus 5% for the so far unkown holes and genus
// of the object). These values are used to reserve space in the
// polyhedron representation `HDS'. If the representation
// supports insertion these values do not restrict the class of
// readable polyhedrons. If the representation does not support
// insertion the object must fit in the reserved sizes.
// If `mode' is set to CGAL_ABSOLUTE the incremental builder uses
// CGAL_ABSOLUTE indexing and the vertices of the old polyhedral surface
// can be used in new facets. Otherwise CGAL_RELATIVE indexing is used
// starting with new indices for the new construction.
Vertex_handle add_vertex( const Point_3& p) {
// adds p to the vertex list.
CGAL_assertion( check_protocoll == 1);
if ( hds.size_of_vertices() >= hds.capacity_of_vertices()) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::"
<< std::endl;
verr << "add_vertex(): capacity error: more than " << new_vertices
<< " vertices added." << std::endl;
m_error = true;
return Vertex_handle();
}
HalfedgeDS_decorator<HDS> decorator(hds);
Vertex_handle v = decorator.vertices_push_back( Vertex(p));
index_to_vertex_map.push_back( v);
decorator.set_vertex_halfedge( v, Halfedge_handle());
push_back_vertex_to_edge_map( Halfedge_handle());
++new_vertices;
return v;
}
// returns handle for the vertex of index i
Vertex_handle vertex( std::size_t i) {
if ( i < new_vertices)
return index_to_vertex_map[i];
return Vertex_handle();
}
Facet_handle begin_facet() {
// starts a facet.
if ( m_error)
return Facet_handle();
CGAL_assertion( check_protocoll == 1);
CGAL_assertion_code( check_protocoll = 2;)
if ( hds.size_of_faces() >= hds.capacity_of_faces()) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::"
<< std::endl;
verr << "begin_facet(): capacity error: more than " << new_vertices
<< " facets added." << std::endl;
m_error = true;
return Facet_handle();
}
// initialize all status variables.
first_vertex = true; // denotes 'no vertex yet'
g1 = Halfedge_handle(); // denotes 'no halfedge yet'
last_vertex = false;
HalfedgeDS_decorator<HDS> decorator(hds);
current_face = decorator.faces_push_back( Face());
return current_face;
}
void add_vertex_to_facet( std::size_t i);
// adds a vertex with index i to the current facet. The first
// point added with `add_vertex()' has the index 0.
Halfedge_handle end_facet() {
// ends a facet.
if ( m_error)
return Halfedge_handle();
CGAL_assertion( check_protocoll == 2);
CGAL_assertion( ! first_vertex);
// cleanup all static status variables
add_vertex_to_facet( w1);
if ( m_error)
return Halfedge_handle();
last_vertex = true;
add_vertex_to_facet( w2);
if ( m_error)
return Halfedge_handle();
CGAL_assertion( check_protocoll == 2);
CGAL_assertion_code( check_protocoll = 1;)
HalfedgeDS_items_decorator<HDS> decorator;
Halfedge_handle h = get_vertex_to_edge_map(w1);
decorator.set_face_halfedge( current_face, h);
++new_faces;
return h;
}
template <class InputIterator>
Halfedge_handle add_facet( InputIterator first, InputIterator beyond) {
// synonym for begin_facet(), a call to add_facet() for each iterator
// value type, and end_facet().
begin_facet();
for ( ; ! m_error && first != beyond; ++first)
add_vertex_to_facet( *first);
if ( m_error)
return Halfedge_handle();
return end_facet();
}
template <class InputIterator>
bool test_facet( InputIterator first, InputIterator beyond) {
// tests if the facet described by the vertex indices in the
// range [first,beyond) can be inserted without creating a
// a non-manifold (and therefore invalid) situation.
// First, create a copy of the indices and close it cyclically
std::vector< std::size_t> indices( first, beyond);
if ( indices.size() < 3)
return false;
indices.push_back( indices[0]);
return test_facet_indices( indices);
}
bool test_facet_indices( std::vector< std::size_t> indices);
void end_surface();
// ends the construction.
bool check_unconnected_vertices();
// returns `true' if unconnected vertices are detected. If `verb'
// is set to `true' debug information about the unconnected
// vertices is printed.
bool remove_unconnected_vertices( Tag_true);
bool remove_unconnected_vertices( Tag_false) {
return ! check_unconnected_vertices();
}
bool remove_unconnected_vertices() {
// returns `true' if all unconnected vertices could be removed
// succesfully.
return remove_unconnected_vertices( Supports_removal());
}
void rollback();
protected:
Halfedge_handle lookup_hole( std::size_t w) {
CGAL_assertion( w < new_vertices);
return lookup_hole( get_vertex_to_edge_map( w));
}
size_type find_vertex( Vertex_handle v) {
// Returns 0 if v == NULL.
if ( v == Vertex_handle() )
return 0;
size_type n = 0;
typename HDS::Vertex_iterator it = hds.vertices_begin();
while ( it != v) {
CGAL_assertion( it != hds.vertices_end());
++n;
++it;
}
n = n - rollback_v;
return n;
}
size_type find_facet( Face_handle f) {
// Returns 0 if f == NULL.
if ( f == Face_handle())
return 0;
size_type n = 0;
typename HDS::Face_iterator it = hds.faces_begin();
while ( it != f) {
CGAL_assertion( it != hds.faces_end());
++n;
++it;
}
n = n - rollback_f;
return n;
}
Halfedge_handle lookup_halfedge( size_type w, size_type v) {
// Pre: 0 <= w,v < new_vertices
// Case a: It exists an halfedge g from w to v:
// g must be a border halfedge and the facet of g->opposite()
// must be set and different from the current facet.
// Set the facet of g to the current facet. Return the
// halfedge pointing to g.
// Case b: It exists no halfedge from w to v:
// Create a new pair of halfedges g and g->opposite().
// Set the facet of g to the current facet and g->opposite()
// to a border halfedge. Assign the vertex references.
// Set g->opposite()->next() to g. Return g->opposite().
typedef typename HDS::Supports_halfedge_vertex
Supports_halfedge_vertex;
Assert_compile_time_tag( Supports_halfedge_vertex(), Tag_true());
CGAL_assertion( w < new_vertices);
CGAL_assertion( v < new_vertices);
CGAL_assertion( ! last_vertex);
HalfedgeDS_items_decorator<HDS> decorator;
Halfedge_handle e = get_vertex_to_edge_map( w);
if ( e != Halfedge_handle()) {
CGAL_assertion( e->vertex() == index_to_vertex_map[w]);
// check that the facet has no self intersections
if ( current_face != Face_handle()
&& current_face == decorator.get_face(e)) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::"
<< std::endl;
verr << "lookup_halfedge(): input error: facet "
<< new_faces << " has a self intersection at vertex "
<< w << "." << std::endl;
m_error = true;
return Halfedge_handle();
}
Halfedge_handle start_edge( e);
do {
if ( e->next()->vertex() == index_to_vertex_map[v]) {
if ( ! e->next()->is_border()) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3"
"<HDS>::" << std::endl;
verr << "lookup_halfedge(): input error: facet "
<< new_faces << " shares a halfedge from "
"vertex " << w << " to vertex " << v
<< " with";
if ( m_verbose && current_face != Face_handle())
verr << " facet "
<< find_facet( decorator.get_face(e->next()))
<< '.' << std::endl;
else
verr << " another facet." << std::endl;
m_error = true;
return Halfedge_handle();
}
CGAL_assertion( ! e->next()->opposite()->is_border());
if ( current_face != Face_handle() && current_face ==
decorator.get_face( e->next()->opposite())) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3"
"<HDS>::" << std::endl;
verr << "lookup_halfedge(): input error: facet "
<< new_faces << " has a self intersection "
"at the halfedge from vertex " << w
<< " to vertex " << v << "." << std::endl;
m_error = true;
return Halfedge_handle();
}
decorator.set_face( e->next(), current_face);
return e;
}
e = e->next()->opposite();
} while ( e != start_edge);
}
// create a new halfedge
if ( hds.size_of_halfedges() >= hds.capacity_of_halfedges()) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::"
<< std::endl;
verr << "lookup_halfedge(): capacity error: more than "
<< new_halfedges << " halfedges added while creating facet"
<< new_faces << '.' << std::endl;
m_error = true;
return Halfedge_handle();
}
e = hds.edges_push_back( Halfedge(), Halfedge());
new_halfedges++;
new_halfedges++;
decorator.set_face( e, current_face);
e->HBase::set_vertex( index_to_vertex_map[v]);
e->HBase::set_next( Halfedge_handle());
decorator.set_prev( e, e->opposite());
e = e->opposite();
e->HBase::set_vertex( index_to_vertex_map[w]);
e->HBase::set_next( e->opposite());
return e;
}
Halfedge_handle lookup_hole( Halfedge_handle e) {
// Halfedge e points to a vertex w. Walk around w to find a hole
// in the facet structure. Report an error if none exist. Return
// the halfedge at this hole that points to the vertex w.
CGAL_assertion( e != Halfedge_handle());
HalfedgeDS_items_decorator<HDS> decorator;
Halfedge_handle start_edge( e);
do {
if ( e->next()->is_border()) {
return e;
}
e = e->next()->opposite();
} while ( e != start_edge);
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::" << std::endl;
verr << "lookup_hole(): input error: at vertex "
<< find_vertex( e->vertex())
<< " a closed surface already exists and facet "
<< new_faces << " is nonetheless adjacent." << std::endl;
if ( m_verbose && current_face != Face_handle()) {
verr << " The closed cycle of facets is:";
do {
if ( ! e->is_border())
verr << " " << find_facet( decorator.get_face(e));
e = e->next()->opposite();
} while ( e != start_edge);
verr << '.' << std::endl;
}
m_error = true;
return Halfedge_handle();
}
};
template < class HDS>
void
Enriched_polyhedron_incremental_builder_3<HDS>::
rollback() {
CGAL_assertion( rollback_v <= hds.size_of_vertices());
CGAL_assertion( rollback_h <= hds.size_of_halfedges());
CGAL_assertion( rollback_f <= hds.size_of_faces());
if ( rollback_v == 0 && rollback_h == 0 && rollback_f == 0) {
hds.clear();
} else {
while ( rollback_v != hds.size_of_vertices())
hds.vertices_pop_back();
CGAL_assertion((( hds.size_of_halfedges() - rollback_h) & 1) == 0);
while ( rollback_h != hds.size_of_halfedges())
hds.edges_pop_back();
while ( rollback_f != hds.size_of_faces())
hds.faces_pop_back();
}
m_error = false;
CGAL_assertion_code( check_protocoll = 0;)
}
template < class HDS> CGAL_MEDIUM_INLINE
void
Enriched_polyhedron_incremental_builder_3<HDS>::
begin_surface( std::size_t v, std::size_t f, std::size_t h, int mode) {
CGAL_assertion( check_protocoll == 0);
CGAL_assertion_code( check_protocoll = 1;)
CGAL_assertion( ! m_error);
if ( mode == CGAL_RELATIVE) {
new_vertices = 0;
new_faces = 0;
new_halfedges = 0;
rollback_v = hds.size_of_vertices();
rollback_f = hds.size_of_faces();
rollback_h = hds.size_of_halfedges();
} else {
new_vertices = hds.size_of_vertices();
new_faces = hds.size_of_faces();
new_halfedges = hds.size_of_halfedges();
rollback_v = 0;
rollback_f = 0;
rollback_h = 0;
}
if ( h == 0) {
// Use the Eulerian equation for connected planar graphs. We do
// not know the number of facets that are holes and we do not
// know the genus of the surface. So we add 12 and a factor of
// 5 percent.
h = int((v + f - 2 + 12) * 2.1);
}
hds.reserve( hds.size_of_vertices() + v,
hds.size_of_halfedges() + h,
hds.size_of_faces() + f);
if ( mode == CGAL_RELATIVE) {
index_to_vertex_map = Random_access_index( hds.vertices_end());
index_to_vertex_map.reserve(v);
initialize_vertex_to_edge_map( v, false);
} else {
index_to_vertex_map = Random_access_index( hds.vertices_begin(),
hds.vertices_end());
index_to_vertex_map.reserve( hds.size_of_vertices() + v);
initialize_vertex_to_edge_map( hds.size_of_vertices() + v, true);
}
}
template < class HDS>
void
Enriched_polyhedron_incremental_builder_3<HDS>::
add_vertex_to_facet( std::size_t v2) {
if ( m_error)
return;
CGAL_assertion( check_protocoll == 2);
if ( v2 >= new_vertices) {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::"
<< std::endl;
verr << "add_vertex_to_facet(): vertex index " << v2
<< " is out-of-range [0," << new_vertices-1 << "]."
<< std::endl;
m_error = true;
return;
}
HalfedgeDS_items_decorator<HDS> decorator;
if ( first_vertex) {
CGAL_assertion( ! last_vertex);
w1 = v2;
first_vertex = false;
return;
}
if ( g1 == Halfedge_handle()) {
CGAL_assertion( ! last_vertex);
gprime = lookup_halfedge( w1, v2);
if ( m_error)
return;
h1 = g1 = gprime->next();
v1 = w2 = v2;
return;
}
// g1, h1, v1, w1, w2 are set. Insert halfedge.
// Lookup v1-->v2
Halfedge_handle hprime;
if ( last_vertex)
hprime = gprime;
else {
hprime = lookup_halfedge( v1, v2);
if ( m_error)
return;
}
Halfedge_handle h2 = hprime->next();
CGAL_assertion( ! last_vertex || h2 == g1);
Halfedge_handle prev = h1->next();
h1->HBase::set_next( h2);
decorator.set_prev( h2, h1);
if ( get_vertex_to_edge_map( v1) == Halfedge_handle()) { // case 1:
h2->opposite()->HBase::set_next( h1->opposite());
decorator.set_prev( h1->opposite(), h2->opposite());
} else { // case 2:
bool b1 = h1->opposite()->is_border();
bool b2 = h2->opposite()->is_border();
if ( b1 && b2) {
Halfedge_handle hole = lookup_hole( v1);
if ( m_error)
return;
CGAL_assertion( hole != Halfedge_handle());
h2->opposite()->HBase::set_next( hole->next());
decorator.set_prev( hole->next(), h2->opposite());
hole->HBase::set_next( h1->opposite());
decorator.set_prev( h1->opposite(), hole);
} else if ( b2) { // case 2.b:
CGAL_assertion( prev->is_border());
h2->opposite()->HBase::set_next( prev);
decorator.set_prev( prev, h2->opposite());
} else if ( b1) { // case 2.c:
CGAL_assertion( hprime->is_border());
hprime->HBase::set_next( h1->opposite());
decorator.set_prev( h1->opposite(), hprime);
} else if ( h2->opposite()->next() == h1->opposite()) {// case 2.d:
// f1 == f2
CGAL_assertion( decorator.get_face( h1->opposite()) ==
decorator.get_face( h2->opposite()));
} else { // case 2.e:
if ( prev == h2) { // case _i:
// nothing to be done, hole is closed.
} else { // case _ii:
CGAL_assertion( prev->is_border());
CGAL_assertion( hprime->is_border());
hprime->HBase::set_next( prev);
decorator.set_prev( prev, hprime);
// Check whether the halfedges around v1 are connected.
// It is sufficient to check it for h1 to prev.
// Assert loop termination:
CGAL_assertion_code( std::size_t k = 0;)
// Look for a hole in the facet complex starting at h1.
Halfedge_handle hole;
Halfedge_handle e = h1;
do {
if ( e->is_border())
hole = e;
e = e->next()->opposite();
CGAL_assertion( k++ < hds.size_of_halfedges());
} while ( e->next() != prev && e != h1);
if ( e == h1) {
// disconnected facet complexes
if ( hole != Halfedge_handle()) {
// The complex can be connected with
// the hole at hprime.
hprime->HBase::set_next( hole->next());
decorator.set_prev( hole->next(), hprime);
hole->HBase::set_next( prev);
decorator.set_prev( prev, hole);
} else {
Verbose_ostream verr( m_verbose);
verr << " " << std::endl;
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<"
"HDS>::" << std::endl;
verr << "add_vertex_to_facet(): input error: "
"disconnected facet complexes at vertex "
<< v1 << ":" << std::endl;
if ( m_verbose && current_face != Face_handle()) {
verr << " involved facets are:";
do {
if ( ! e->is_border())
verr << " " << find_facet(
decorator.get_face(e));
e = e->next()->opposite();
} while ( e != h1);
verr << " (closed cycle) and";
e = hprime;
do {
if ( ! e->is_border())
verr << " " << find_facet(
decorator.get_face(e));
} while ( e != hprime);
verr << "." << std::endl;
}
m_error = true;
return;
}
}
}
}
}
if ( h1->vertex() == index_to_vertex_map[v1])
set_vertex_to_edge_map( v1, h1);
CGAL_assertion( h1->vertex() == index_to_vertex_map[v1]);
h1 = h2;
v1 = v2;
}
template < class HDS>
bool
Enriched_polyhedron_incremental_builder_3<HDS>::
test_facet_indices( std::vector< std::size_t> indices) {
typedef typename HDS::Supports_halfedge_vertex Supports_halfedge_vertex;
Assert_compile_time_tag( Supports_halfedge_vertex(), Tag_true());
// tests if the facet described by the vertex indices can be inserted
// without creating a a non-manifold (and therefore invalid) situation.
// indices are cyclically closed once.
std::size_t n = indices.size() - 1;
// Test if a vertex is not twice in the indices
for ( std::size_t i = 0; i < n; ++i) {
CGAL_precondition( indices[i] < new_vertices);
// check if vertex indices[i] is already in the sequence [0..i-1]
for ( std::size_t k = 0; k+1 < i; ++k) {
if ( indices[k] == indices[i])
return false;
}
}
// Test non-manifold edges
for ( std::size_t i = 0; i < n; ++i) {
// edge goes from vertex indices[i] to indices[i+1]
// we know already that the edge is only once in the sequence
// (otherwise the end-vertices would be twice in the sequence too)
// check if edge is already in the HDS and is not border edge
Halfedge_handle v = get_vertex_to_edge_map(indices[i]);
Vertex_handle w = index_to_vertex_map[indices[i+1]];
if ( v != Halfedge_handle()
&& get_vertex_to_edge_map(indices[i+1]) != Halfedge_handle()) {
// cycle through halfedge-loop and find edge to indices[i+1]
Halfedge_handle vstart = v;
do {
v = v->next()->opposite();
} while ( v->next()->vertex() != w && v != vstart);
if ( v->next()->vertex() == w && ! v->next()->is_border())
return false;
}
}
// test non-manifold vertices
for ( std::size_t i = 0; i < n; ++i) {
// since we don't allow duplicates in indices[..] and we
// tested for non-manifold edges already, we just need to check
// if the vertex indices[i] is not a closed manifold yet.
Halfedge_handle v = get_vertex_to_edge_map(indices[i]);
if ( v != Halfedge_handle()) {
Halfedge_handle vstart = v;
do {
v = v->next()->opposite();
} while ( ! v->is_border() && v != vstart);
if ( ! v->is_border())
return false;
}
}
return true;
}
template < class HDS> CGAL_MEDIUM_INLINE
void
Enriched_polyhedron_incremental_builder_3<HDS>::
end_surface() {
if ( m_error)
return;
CGAL_assertion( check_protocoll == 1);
CGAL_assertion_code( check_protocoll = 0;)
}
template < class HDS>
bool
Enriched_polyhedron_incremental_builder_3<HDS>::
check_unconnected_vertices() {
if ( m_error)
return false;
bool unconnected = false;
Verbose_ostream verr( m_verbose);
for ( std::size_t i = 0; i < new_vertices; i++) {
if ( get_vertex_to_edge_map( i) == Halfedge_handle()) {
verr << "CGAL::Enriched_polyhedron_incremental_builder_3<HDS>::\n"
<< "check_unconnected_vertices( verb = true): "
<< "vertex " << i << " is unconnected." << std::endl;
unconnected = true;
}
}
return unconnected;
}
template < class HDS>
bool
Enriched_polyhedron_incremental_builder_3<HDS>::
remove_unconnected_vertices( Tag_true) {
if ( m_error)
return true;
for( std::size_t i = 0; i < new_vertices; i++) {
if( get_vertex_to_edge_map( i) == Halfedge_handle()) {
hds.vertices_erase( index_to_vertex_map[i]);
}
}
return true;
}
CGAL_END_NAMESPACE
#endif // CGAL_USE_POLYHEDRON_DESIGN_ONE //
#endif // CGAL_Enriched_polyhedron_incremental_builder_3_H //
// EOF //

View File

@ -0,0 +1,898 @@
// -*- tab-width: 2 -*-
///////////////////////////////////////////////////////////////////////////
// //
// Class: Enriched_polyhedron //
// //
///////////////////////////////////////////////////////////////////////////
#ifndef _POLYGON_MESH_
#define _POLYGON_MESH_
// CGAL stuff
#include <CGAL/Cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <list>
#include <GL/gl.h>
// a refined facet with a normal and a tag
template <class Refs, class T, class P, class Norm>
class Enriched_facet : public CGAL::HalfedgeDS_face_base<Refs, T>
{
// tag
int m_tag;
// normal
Norm m_normal;
public:
// life cycle
// no constructors to repeat, since only
// default constructor mandatory
Enriched_facet()
{
}
// tag
const int& tag() const { return m_tag; }
int& tag() { return m_tag; }
// normal
typedef Norm Normal_3;
Normal_3& normal() { return m_normal; }
const Normal_3& normal() const { return m_normal; }
};
// a refined halfedge with a general tag and
// a binary tag to indicate wether it belongs
// to the control mesh or not
template <class Refs, class Tprev, class Tvertex, class Tface, class Norm>
class Enriched_halfedge : public CGAL::HalfedgeDS_halfedge_base<Refs,Tprev,Tvertex,Tface>
{
private:
// general purpose tag
int m_tag;
// option for edge superimposing
bool m_control_edge;
bool m_sharp;
public:
// life cycle
Enriched_halfedge()
{
m_control_edge = true;
m_sharp = false;
}
// tag
const int& tag() const { return m_tag; }
int& tag() { return m_tag; }
void tag(const int& t) { m_tag = t; }
// control edge
bool& control_edge() { return m_control_edge; }
const bool& control_edge() const { return m_control_edge; }
// sharp
bool& sharp() { return m_sharp; }
const bool& sharp() const { return m_sharp; }
};
// a refined vertex with a normal and a tag
template <class Refs, class T, class P, class Norm>
class Enriched_vertex : public CGAL::HalfedgeDS_vertex_base<Refs, T, P>
{
// tag
int m_tag;
// normal
Norm m_normal;
public:
// life cycle
Enriched_vertex() {}
// repeat mandatory constructors
Enriched_vertex(const P& pt)
: CGAL::HalfedgeDS_vertex_base<Refs, T, P>(pt)
{
}
// normal
typedef Norm Normal_3;
Normal_3& normal() { return m_normal; }
const Normal_3& normal() const { return m_normal; }
// tag
int& tag() { return m_tag; }
const int& tag() const { return m_tag; }
void tag(const int& t) { m_tag = t; }
};
// A redefined items class for the Polyhedron_3
// with a refined vertex class that contains a
// member for the normal vector and a refined
// facet with a normal vector instead of the
// plane equation (this is an alternative
// solution instead of using
// Polyhedron_traits_with_normals_3).
struct Enriched_items : public CGAL::Polyhedron_items_3
{
// wrap vertex
template<class Refs, class Traits> struct Vertex_wrapper
{
typedef typename Traits::Point_3 Point;
typedef typename Traits::Vector_3 Normal;
typedef Enriched_vertex<Refs,
CGAL::Tag_true,
Point,
Normal> Vertex;
};
// wrap face
template<class Refs, class Traits> struct Face_wrapper
{
typedef typename Traits::Point_3 Point;
typedef typename Traits::Vector_3 Normal;
typedef Enriched_facet<Refs,
CGAL::Tag_true,
Point,
Normal> Face;
};
// wrap halfedge
template<class Refs, class Traits> struct Halfedge_wrapper
{
typedef typename Traits::Vector_3 Normal;
typedef Enriched_halfedge<Refs,
CGAL::Tag_true,
CGAL::Tag_true,
CGAL::Tag_true,
Normal> Halfedge;
};
};
static const double PI = 3.1415926535897932384626;
// compute facet normal
struct Facet_normal // (functor)
{
template<class Facet> void operator()(Facet& f)
{
typename Facet::Normal_3 sum = CGAL::NULL_VECTOR;
typename Facet::Halfedge_around_facet_circulator h = f.facet_begin();
do
{
typename Facet::Normal_3 normal = CGAL::cross_product(h->next()->vertex()->point() - h->vertex()->point(), h->next()->next()->vertex()->point() - h->next()->vertex()->point());
double sqnorm = normal * normal;
if (sqnorm != 0)
normal = normal / (float)std::sqrt(sqnorm);
sum = sum + normal;
} while (++h != f.facet_begin());
float sqnorm = sum * sum;
if (sqnorm != 0.0)
f.normal() = sum / std::sqrt(sqnorm);
else
{
f.normal() = CGAL::NULL_VECTOR;
// TRACE("degenerate face\n");
}
}
};
// compute vertex normal
struct Vertex_normal // (functor)
{
template<class Vertex> void operator()(Vertex& v)
{
typename Vertex::Normal_3 normal = CGAL::NULL_VECTOR;
typename Vertex::Halfedge_around_vertex_const_circulator pHalfedge =
v.vertex_begin();
typename Vertex::Halfedge_around_vertex_const_circulator begin =
pHalfedge;
CGAL_For_all(pHalfedge,begin)
if(!pHalfedge->is_border())
normal = normal + pHalfedge->facet()->normal();
float sqnorm = normal * normal;
if (sqnorm != 0.0f)
v.normal() = normal / (float)std::sqrt(sqnorm);
else
v.normal() = CGAL::NULL_VECTOR;
}
};
//*********************************************************
template <class kernel, class items>
class Enriched_polyhedron : public CGAL::Polyhedron_3<kernel,items>
{
public :
typedef typename kernel::FT FT;
typedef typename kernel::Point_3 Point;
typedef typename kernel::Vector_3 Vector;
typedef typename kernel::Iso_cuboid_3 Iso_cuboid;
typedef CGAL::Polyhedron_3<kernel,items> Base;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Vertex_iterator Vertex_iterator;
typedef typename Base::Halfedge_handle Halfedge_handle;
typedef typename Base::Halfedge_iterator Halfedge_iterator;
typedef typename Base::Halfedge_around_facet_circulator Halfedge_around_facet_circulator;
typedef typename Base::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator;
typedef typename Base::Edge_iterator Edge_iterator;
typedef typename Base::Facet_iterator Facet_iterator;
typedef typename Base::Facet_handle Facet_handle;
typedef typename Base::Facet Facet;
using Base::vertices_begin;
using Base::vertices_end;
using Base::edges_begin;
using Base::edges_end;
using Base::halfedges_begin;
using Base::halfedges_end;
using Base::facets_begin;
using Base::facets_end;
using Base::size_of_halfedges;
using Base::size_of_facets;
using Base::size_of_vertices;
private :
// bounding box
Iso_cuboid m_bbox;
// type
bool m_pure_quad;
bool m_pure_triangle;
public :
enum Vertex_type { SMOOTH, // 0 sharp edge
DART, // 1
CREASE_REGULAR, // 2 - with two non-sharp edges on each side
CREASE_IRREGULAR, // 2
CORNER }; // 3 and more
public :
// life cycle
Enriched_polyhedron()
{
m_pure_quad = false;
m_pure_triangle = false;
}
virtual ~Enriched_polyhedron()
{
}
// type
bool is_pure_triangle() { return m_pure_triangle; }
bool is_pure_quad() { return m_pure_quad; }
// normals (per facet, then per vertex)
void compute_normals_per_facet()
{
std::for_each(facets_begin(),facets_end(),::Facet_normal());
}
void compute_normals_per_vertex()
{
std::for_each(vertices_begin(),vertices_end(),::Vertex_normal());
}
void compute_normals()
{
compute_normals_per_facet();
compute_normals_per_vertex();
}
// bounding box
Iso_cuboid& bbox() { return m_bbox; }
const Iso_cuboid bbox() const { return m_bbox; }
// compute bounding box
void compute_bounding_box()
{
CGAL_assertion(size_of_vertices() != 0);
FT xmin,xmax,ymin,ymax,zmin,zmax;
Vertex_iterator pVertex = vertices_begin();
xmin = xmax = pVertex->point().x();
ymin = ymax = pVertex->point().y();
zmin = zmax = pVertex->point().z();
for(;
pVertex != vertices_end();
pVertex++)
{
const Point& p = pVertex->point();
xmin = std::min(xmin,p.x());
ymin = std::min(ymin,p.y());
zmin = std::min(zmin,p.z());
xmax = std::max(xmax,p.x());
ymax = std::max(ymax,p.y());
zmax = std::max(zmax,p.z());
}
m_bbox = Iso_cuboid(xmin,ymin,zmin,
xmax,ymax,zmax);
}
// bounding box
FT xmin() { return m_bbox.xmin(); }
FT xmax() { return m_bbox.xmax(); }
FT ymin() { return m_bbox.ymin(); }
FT ymax() { return m_bbox.ymax(); }
FT zmin() { return m_bbox.zmin(); }
FT zmax() { return m_bbox.zmax(); }
Point center()
{
FT cx = (FT)0.5 * (xmin() + xmax());
FT cy = (FT)0.5 * (ymin() + ymax());
FT cz = (FT)0.5 * (zmin() + zmax());
return Point(cx,cy,cz);
}
FT size()
{
FT dx = xmax() - xmin();
FT dy = ymax() - ymin();
FT dz = zmax() - zmin();
return std::max(dx,std::max(dy,dz));
}
// degree of a face
static unsigned int degree(Facet_handle pFace)
{
return CGAL::circulator_size(pFace->facet_begin());
}
// valence of a vertex
static unsigned int valence(Vertex_handle pVertex)
{
return CGAL::circulator_size(pVertex->vertex_begin());
}
// check wether a vertex is on a boundary or not
static bool is_border(Vertex_handle pVertex)
{
Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
if(pHalfEdge == NULL) // isolated vertex
return true;
Halfedge_around_vertex_circulator d = pHalfEdge;
CGAL_For_all(pHalfEdge,d)
if(pHalfEdge->is_border())
return true;
return false;
}
// get any border halfedge attached to a vertex
Halfedge_handle get_border_halfedge(Vertex_handle pVertex)
{
Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
Halfedge_around_vertex_circulator d = pHalfEdge;
CGAL_For_all(pHalfEdge,d)
if(pHalfEdge->is_border())
return pHalfEdge;
return NULL;
}
// tag all halfedges
void tag_halfedges(const int tag)
{
for(Halfedge_iterator pHalfedge = halfedges_begin();
pHalfedge != halfedges_end();
pHalfedge++)
pHalfedge->tag(tag);
}
// tag all facets
void tag_facets(const int tag)
{
for(Facet_iterator pFacet = facets_begin();
pFacet != facets_end();
pFacet++)
pFacet->tag() = tag;
}
// set index for all vertices
void set_index_vertices()
{
int index = 0;
for(Vertex_iterator pVertex = vertices_begin();
pVertex != vertices_end();
pVertex++)
pVertex->tag(index++);
}
// set index for all edges
void set_index_edges()
{
int index = 0;
for(Edge_iterator he = edges_begin();
he != edges_end();
he++)
{
he->tag(index);
he->opposite()->tag(index);
index++;
}
}
// is pure degree ?
bool is_pure_degree(unsigned int d)
{
for(Facet_iterator pFace = facets_begin();
pFace != facets_end();
pFace++)
if(degree(pFace) != d)
return false;
return true;
}
// compute type
void compute_type()
{
m_pure_quad = is_pure_degree(4);
m_pure_triangle = is_pure_degree(3);
}
// compute facet center
void compute_facet_center(Facet_handle pFace,
Point& center)
{
Halfedge_around_facet_circulator pHalfEdge = pFace->facet_begin();
Halfedge_around_facet_circulator end = pHalfEdge;
Vector vec(0.0,0.0,0.0);
int degree = 0;
CGAL_For_all(pHalfEdge,end)
{
vec = vec + (pHalfEdge->vertex()->point()-CGAL::ORIGIN);
degree++;
}
center = CGAL::ORIGIN + (vec/FT(degree));
}
void gl_draw_direct_triangles(bool smooth_shading,
bool use_normals)
{
// draw triangles
::glBegin(GL_TRIANGLES);
Facet_iterator pFacet = facets_begin();
for(;pFacet != facets_end();pFacet++)
gl_draw_facet(pFacet,smooth_shading,use_normals);
::glEnd(); // end polygon assembly
}
void gl_draw_direct(bool smooth_shading,
bool use_normals)
{
// draw polygons
Facet_iterator pFacet = facets_begin();
for(;pFacet != facets_end();pFacet++)
{
// begin polygon assembly
::glBegin(GL_POLYGON);
gl_draw_facet(pFacet,smooth_shading,use_normals);
::glEnd(); // end polygon assembly
}
}
void gl_draw_facet(Facet_handle pFacet,
bool smooth_shading,
bool use_normals)
{
// one normal per face
if(use_normals && !smooth_shading)
{
const typename Facet::Normal_3& normal = pFacet->normal();
::glNormal3f(normal[0],normal[1],normal[2]);
}
// revolve around current face to get vertices
Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
do
{
// one normal per vertex
if(use_normals && smooth_shading)
{
const typename Facet::Normal_3& normal = pHalfedge->vertex()->normal();
::glNormal3f(normal[0],normal[1],normal[2]);
}
// polygon assembly is performed per vertex
const Point& point = pHalfedge->vertex()->point();
::glVertex3d(point[0],point[1],point[2]);
}
while(++pHalfedge != pFacet->facet_begin());
}
// superimpose edges
void superimpose_edges(bool skip_ordinary_edges = true,
bool skip_control_edges = false)
{
::glBegin(GL_LINES);
for(Edge_iterator h = edges_begin();
h != edges_end();
h++)
{
if(h->sharp())
continue;
// ignore this edges
if(skip_ordinary_edges && !h->control_edge())
continue;
// ignore control edges
if(skip_control_edges && h->control_edge())
continue;
// assembly and draw line segment
const Point& p1 = h->prev()->vertex()->point();
const Point& p2 = h->vertex()->point();
::glVertex3f(p1[0],p1[1],p1[2]);
::glVertex3f(p2[0],p2[1],p2[2]);
}
::glEnd();
}
bool is_sharp(Halfedge_handle he,
const double angle_sharp)
{
Facet_handle f1 = he->facet();
Facet_handle f2 = he->opposite()->facet();
if(f1 == NULL || f2 == NULL)
return false;
const Vector& n1 = f1->normal();
const Vector& n2 = f2->normal();
if(angle_deg(n1,n2) > angle_sharp)
return true;
else
return false;
}
Vertex_type type(Vertex_handle v)
{
unsigned int nb = nb_sharp_edges(v);
switch(nb)
{
case 0:
return SMOOTH;
case 1:
return DART;
case 2: // crease vertex - may be regular or not
return crease_type(v);
default: // 3 and more
return CORNER;
}
}
// regular crease vertex must have valence 6, with
// exactly two smooth edges on each side.
Vertex_type crease_type(Vertex_handle v)
{
if(valence(v) != 6)
return CREASE_IRREGULAR;
// valence = 6 - let us check regularity
// pick first sharp edge
Halfedge_around_vertex_circulator he = v->vertex_begin();
Halfedge_around_vertex_circulator end = he;
CGAL_For_all(he,end)
if(he->sharp())
break;
// next two must be smooth
for(int i=0;i<2;i++)
if(++he->sharp())
return CREASE_IRREGULAR;
// next one must be sharp
if(!++he->sharp())
return CREASE_IRREGULAR;
// next two must be smooth
for(int i=0;i<2;i++)
if(++he->sharp())
return CREASE_IRREGULAR;
return CREASE_REGULAR;
}
// return true if succeds
void incident_points_on_crease(Vertex_handle v,
Point& a,
Point& b)
{
#ifdef _DEBUG
Vertex_type vertex_type = type(v,angle_sharp);
ASSERT(vertex_type == CREASE_IRREGULAR || vertex_type == CREASE_REGULAR);
#endif // _DEBUG
// pick first sharp edge
Halfedge_around_vertex_circulator he = v->vertex_begin();
Halfedge_around_vertex_circulator end = he;
CGAL_For_all(he,end)
if(he->sharp())
break;
a = he->opposite()->vertex()->point();
// pick next sharp edge
he++;
CGAL_For_all(he,end)
if(he->sharp())
break;
b = he->opposite()->vertex()->point();
}
// return nb of sharp edges incident to v
unsigned int nb_sharp_edges(Vertex_handle v)
{
Halfedge_around_vertex_circulator he = v->vertex_begin();
Halfedge_around_vertex_circulator end = he;
unsigned int nb_sharp_edges = 0;
CGAL_For_all(he,end)
if(he->sharp())
nb_sharp_edges++;
return nb_sharp_edges;
}
// Angle between two vectors (in degrees)
// we use this formula
// uv = |u||v| cos(u,v)
// u ^ v = w
// |w| = |u||v| |sin(u,v)|
//**************************************************
static double angle_deg(const Vector &u,
const Vector &v)
{
static const double conv = 1.0/PI*180.0;
return conv * angle_rad(u,v);
}
static FT len(const Vector &v)
{
return (FT)std::sqrt(CGAL_NTS to_double(v*v));
}
// Angle between two vectors (in rad)
// uv = |u||v| cos(u,v)
// u ^ v = w
// |w| = |u||v| |sin(u,v)|
//**************************************************
static double angle_rad(const Vector &u,
const Vector &v)
{
// check
double product = len(u)*len(v);
if(product == 0)
return 0.0;
// cosine
double dot = (u*v);
double cosine = dot / product;
// sine
Vector w = CGAL::cross_product(u,v);
double AbsSine = len(w) / product;
if(cosine >= 0)
return std::asin(fix_sine(AbsSine));
else
return PI-std::asin(fix_sine(AbsSine));
}
//**********************************************
// fix sine
//**********************************************
static double fix_sine(double sine)
{
if(sine >= 1)
return 1;
else
if(sine <= -1)
return -1;
else
return sine;
}
unsigned int tag_sharp_edges(const double angle_sharp)
{
unsigned int nb = 0;
for(Halfedge_iterator he = edges_begin();
he != edges_end();
he++)
{
const bool tag = is_sharp(he,angle_sharp);
he->sharp() = tag;
he->opposite()->sharp() = tag;
nb += tag ? 1 : 0;
}
return nb;
}
// draw edges
void gl_draw_sharp_edges(const float line_width,
unsigned char r,
unsigned char g,
unsigned char b)
{
::glLineWidth(line_width);
::glColor3ub(r,g,b);
::glBegin(GL_LINES);
for(Halfedge_iterator he = edges_begin();
he != edges_end();
he++)
{
if(he->sharp())
{
const Point& a = he->opposite()->vertex()->point();
const Point& b = he->vertex()->point();
::glVertex3d(a[0],a[1],a[2]);
::glVertex3d(b[0],b[1],b[2]);
}
}
::glEnd();
}
void gl_draw_boundaries()
{
::glBegin(GL_LINES);
for(Halfedge_iterator he = halfedges_begin();
he != halfedges_end();
he++)
{
if(he->is_border())
{
const Point& a = he->vertex()->point();
const Point& b = he->opposite()->vertex()->point();
::glVertex3d(a.x(),a.y(),a.z());
::glVertex3d(b.x(),b.y(),b.z());
}
}
::glEnd();
}
// draw bounding box
void gl_draw_bounding_box()
{
::glBegin(GL_LINES);
// along x axis
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax());
// along y axis
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax());
// along z axis
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmin(),m_bbox.ymax(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymin(),m_bbox.zmax());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmin());
::glVertex3f(m_bbox.xmax(),m_bbox.ymax(),m_bbox.zmax());
::glEnd();
}
// count #boundaries
unsigned int nb_boundaries()
{
unsigned int nb = 0;
tag_halfedges(0);
for(Halfedge_iterator he = halfedges_begin();
he != halfedges_end();
he++)
{
if(he->is_border() && he->tag() == 0)
{
nb++;
Halfedge_handle curr = he;
do
{
curr = curr->next();
curr->tag(1);
}
while(curr != he);
}
}
return nb;
}
// tag component
void tag_component(Facet_handle pSeedFacet,
const int tag_free,
const int tag_done)
{
pSeedFacet->tag() = tag_done;
std::list<Facet_handle> facets;
facets.push_front(pSeedFacet);
while(!facets.empty())
{
Facet_handle pFacet = facets.front();
facets.pop_front();
pFacet->tag() = tag_done;
Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
Halfedge_around_facet_circulator end = pHalfedge;
CGAL_For_all(pHalfedge,end)
{
Facet_handle pNFacet = pHalfedge->opposite()->facet();
if(pNFacet != NULL && pNFacet->tag() == tag_free)
{
facets.push_front(pNFacet);
pNFacet->tag() = tag_done;
}
}
}
}
// count #components
unsigned int nb_components()
{
unsigned int nb = 0;
tag_facets(0);
for(Facet_iterator pFacet = facets_begin();
pFacet != facets_end();
pFacet++)
{
if(pFacet->tag() == 0)
{
nb++;
tag_component(pFacet,0,1);
}
}
return nb;
}
// compute the genus
// V - E + F + B = 2 (C - G)
// C -> #connected components
// G : genus
// B : #boundaries
int genus()
{
int c = nb_components();
int b = nb_boundaries();
int v = size_of_vertices();
int e = size_of_halfedges()/2;
int f = size_of_facets();
return genus(c,v,f,e,b);
}
int genus(int c,
int v,
int f,
int e,
int b)
{
return (2*c+e-b-f-v)/2;
}
};
#endif // _POLYGON_MESH_

View File

@ -0,0 +1,386 @@
#ifndef PWS_LOOP_SUBDIVISION
#define PWS_LOOP_SUBDIVISION
#include "enriched_polyhedron.h"
#include "builder.h"
template <class HDS,class Polyhedron,class Kernel>
class CModifierLoop : public CGAL::Modifier_base<HDS>
{
private:
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;
typedef typename HDS::Face_handle Face_handle;
typedef typename HDS::Vertex_handle Vertex_handle;
typedef typename HDS::Halfedge_handle Halfedge_handle;
typedef typename CGAL::Enriched_polyhedron_incremental_builder_3<HDS> Builder;
typedef typename Polyhedron::Halfedge_around_vertex_circulator HV_circulator;
typedef typename Polyhedron::Vertex_type Vertex_type;
Polyhedron *m_pMesh;
bool m_piecewise_smooth;
public:
// life cycle
CModifierLoop(Polyhedron *pMesh,
const bool piecewise_smooth = true)
{
CGAL_assertion(pMesh != NULL);
m_pMesh = pMesh;
m_piecewise_smooth = piecewise_smooth;
}
~CModifierLoop() {}
// subdivision
void operator()( HDS& hds)
{
Builder builder(hds,true);
builder.begin_surface(3,1,6);
add_vertices(builder);
add_facets(builder);
builder.end_surface();
}
// add vertices
void add_vertices(Builder &builder)
{
// put original vertices
int index = 0;
typename Polyhedron::Vertex_iterator v;
for(v = m_pMesh->vertices_begin();
v != m_pMesh->vertices_end();
v++)
{
v->tag(index);
// border vertices
if(Polyhedron::is_border(v))
{
const Point& curr = v->point();
const typename Polyhedron::Halfedge_handle& he = m_pMesh->get_border_halfedge(v);
CGAL_assertion(he != NULL);
const Point& next = he->next()->vertex()->point();
const Point& prev = he->prev()->vertex()->point();
Point p = CGAL::ORIGIN + 0.25 *(prev - CGAL::ORIGIN) +
0.50 *(curr - CGAL::ORIGIN) +
0.25 *(next - CGAL::ORIGIN);
builder.add_vertex(p);
} // end is border
else
{
Vertex_type vertex_type = m_pMesh->type(v);
switch(vertex_type)
{
// corner - do not move
case Polyhedron::CORNER:
builder.add_vertex(v->point());
break;
// crease type (regular or irregular)
case Polyhedron::CREASE_REGULAR:
case Polyhedron::CREASE_IRREGULAR:
{
const Point& curr = v->point();
// get two other points along crease
// prev ==== curr ==== next
Point prev,next;
m_pMesh->incident_points_on_crease(v,prev,next);
Point p = CGAL::ORIGIN + 0.125 *(prev - CGAL::ORIGIN) +
0.750 *(curr - CGAL::ORIGIN) +
0.125 *(next - CGAL::ORIGIN);
builder.add_vertex(p);
break;
}
// smooth or dart
case Polyhedron::SMOOTH:
case Polyhedron::DART:
{
unsigned int n = Polyhedron::valence(v);
FT w = alpha((FT)n);
FT sum_weights = w + (FT)n;
// new position
Vector pos = w/sum_weights * (v->point() - CGAL::ORIGIN);
HV_circulator he = v->vertex_begin();
HV_circulator end = he;
CGAL_For_all(he,end)
{
// add edge-vertex contribution
const Point& point = he->opposite()->vertex()->point();
pos = pos + (point - CGAL::ORIGIN) / sum_weights;
}
builder.add_vertex(CGAL::ORIGIN + pos);
} // end smooth or dart
} // end switch
} // end !is border
index++;
}
// as many as #edges
m_pMesh->tag_halfedges(-1);
typename Polyhedron::Halfedge_iterator he;
for(he = m_pMesh->halfedges_begin();
he != m_pMesh->halfedges_end();
he++)
{
if(he->tag() != -1)
continue;
if(he->is_border() || he->opposite()->is_border())
{
const Point& p1 = he->vertex()->point();
const Point& p2 = he->opposite()->vertex()->point();
Point point = CGAL::ORIGIN + 0.5 * ((p1 - CGAL::ORIGIN) +
(p2 - CGAL::ORIGIN));
builder.add_vertex(point);
}
else
{
int mask_type = 1; // smooth by default
Vertex_type type_v1;
Vertex_type type_v2;
// if sharp, get type of end vertices to deduce mask type
if(he->sharp())
{
Vertex_handle v1 = he->vertex();
Vertex_handle v2 = he->opposite()->vertex();
type_v1 = m_pMesh->type(v1);
type_v2 = m_pMesh->type(v2);
if(type_v1 == Polyhedron::DART || type_v2 == Polyhedron::DART)
mask_type = 1; // smooth
else
{
if(type_v1 == type_v2) // all diagonals of minor in table 1
mask_type = 2; // regular crease case
else
if(type_v1 == Polyhedron::CREASE_REGULAR ||
type_v2 == Polyhedron::CREASE_REGULAR)
mask_type = 3; // irregular crease case
else
{
mask_type = 2; // regular crease case
// this is all must remain
CGAL_assertion( (type_v1 == Polyhedron::CREASE_IRREGULAR && type_v2 == Polyhedron::CORNER) ||
(type_v2 == Polyhedron::CREASE_IRREGULAR && type_v1 == Polyhedron::CORNER));
}
}
}
const Point& p1 = he->vertex()->point();
const Point& p2 = he->opposite()->vertex()->point();
switch(mask_type)
{
case 1:
{
const Point& p3 = he->next()->vertex()->point();
const Point& p4 = he->opposite()->next()->vertex()->point();
Point p = CGAL::ORIGIN + 0.375 *(p1 - CGAL::ORIGIN) +
0.375 *(p2 - CGAL::ORIGIN) +
0.125 *(p3 - CGAL::ORIGIN) +
0.125 *(p4 - CGAL::ORIGIN);
builder.add_vertex(p);
break;
}
case 2: // midpoint
{
builder.add_vertex(CGAL::midpoint(p1,p2));
break;
}
case 3: // 5/3 (regular vertex gets weight 5)
{
// 0.625
Point regular,irregular;
if(type_v1 == Polyhedron::CREASE_REGULAR)
{
regular = p1;
irregular = p2;
CGAL_assertion(type_v2 == Polyhedron::CREASE_IRREGULAR);
}
else
{
regular = p2;
irregular = p1;
CGAL_assertion(type_v1 == Polyhedron::CREASE_IRREGULAR);
CGAL_assertion(type_v2 == Polyhedron::CREASE_REGULAR);
}
Point p = CGAL::ORIGIN + 0.625 *(regular - CGAL::ORIGIN) +
0.375 *(irregular - CGAL::ORIGIN);
builder.add_vertex(p);
}
}
}
// put vertex indices on both halfedges
he->tag(index);
he->opposite()->tag(index);
index++;
}
}
// add facets
void add_facets(Builder &builder)
{
typename Polyhedron::Facet_iterator f;
for(f = m_pMesh->facets_begin();
f != m_pMesh->facets_end();
f++)
{
unsigned int degree = Polyhedron::degree(f);
CGAL_assertion(degree == 3);
typename Polyhedron::Halfedge_handle he = f->halfedge();
int i0 = he->tag();
int i1 = he->vertex()->tag();
int i2 = he->next()->tag();
int i3 = he->next()->vertex()->tag();
int i4 = he->next()->next()->tag();
int i5 = he->next()->next()->vertex()->tag();
// create 4 triangles
builder.begin_facet();
builder.add_vertex_to_facet(i1);
builder.add_vertex_to_facet(i2);
builder.add_vertex_to_facet(i0);
const Halfedge_handle& h1 = builder.end_facet();
builder.begin_facet();
builder.add_vertex_to_facet(i3);
builder.add_vertex_to_facet(i4);
builder.add_vertex_to_facet(i2);
const Halfedge_handle& h2 = builder.end_facet();
builder.begin_facet();
builder.add_vertex_to_facet(i5);
builder.add_vertex_to_facet(i0);
builder.add_vertex_to_facet(i4);
const Halfedge_handle& h3 = builder.end_facet();
// center face
builder.begin_facet();
builder.add_vertex_to_facet(i0);
builder.add_vertex_to_facet(i2);
builder.add_vertex_to_facet(i4);
Halfedge_handle h4 = builder.end_facet();
h1->next()->next()->control_edge() = false;
h2->next()->next()->control_edge() = false;
h3->next()->next()->control_edge() = false;
h4->control_edge() = false;
h4->next()->control_edge() = false;
h4->next()->next()->control_edge() = false;
h1->next()->next()->sharp() = false;
h2->next()->next()->sharp() = false;
h3->next()->next()->sharp() = false;
h4->sharp() = false;
h4->next()->sharp() = false;
h4->next()->next()->sharp() = false;
// propagate control edge tags
h1->control_edge() = he->control_edge();
h3->next()->control_edge() = he->control_edge();
h1->next()->control_edge() = he->next()->control_edge();
h2->control_edge() = he->next()->control_edge();
h2->next()->control_edge() = he->prev()->control_edge();
h3->control_edge() = he->prev()->control_edge();
// propagate sharp tags
h1->sharp() = he->sharp();
h3->next()->sharp() = he->sharp();
h1->next()->sharp() = he->next()->sharp();
h2->sharp() = he->next()->sharp();
h2->next()->sharp() = he->prev()->sharp();
h3->sharp() = he->prev()->sharp();
}
}
static FT alpha(const FT n)
{
FT t = (3.0 + 2.0 * cos(2 * 3.1415926535 / n));
FT a = 5.0 / 8.0 - t*t / 64.0;
return n * (1-a) / a;
}
// smooth vertex positions
static void smooth(Polyhedron *pMesh)
{
CGAL_assertion(pMesh != NULL);
// alloc position vectors
std::vector<Vector> pos(pMesh->size_of_vertices());
// compute new positions
unsigned int index = 0;
typename Polyhedron::Vertex_iterator v;
for(v = pMesh->vertices_begin();
v != pMesh->vertices_end();
v++)
{
// border vertices
if(Polyhedron::is_border(v))
{
const Point& curr = v->point();
const typename Polyhedron::Halfedge_handle& he = pMesh->get_border_halfedge(v);
CGAL_assertion(he != NULL);
const Point& next = he->next()->vertex()->point();
const Point& prev = he->prev()->vertex()->point();
pos[index] = 0.25 *(prev - CGAL::ORIGIN) +
0.50 *(curr - CGAL::ORIGIN) +
0.25 *(next - CGAL::ORIGIN);
} // end is border
else
{
unsigned int n = Polyhedron::valence(v);
FT w = alpha((FT)n);
FT sum_weights = w + (FT)n;
// new position
pos[index] = w/sum_weights * (v->point() - CGAL::ORIGIN);
// rotate around vertex to compute new position
HV_circulator he = v->vertex_begin();
HV_circulator end = he;
CGAL_For_all(he,end)
{
// add edge-vertex contribution
const Point& point = he->opposite()->vertex()->point();
pos[index] = pos[index] + (point - CGAL::ORIGIN) / sum_weights;
}
} // end !is border
index++;
}
// set new positions
for(v = pMesh->vertices_begin(), index = 0;
v != pMesh->vertices_end();
v++,index++)
v->point() = CGAL::ORIGIN + pos[index];
}
};
template <class Polyhedron,class Kernel>
class CSubdivider_loop
{
public:
typedef typename Polyhedron::HalfedgeDS HalfedgeDS;
CSubdivider_loop() {}
~CSubdivider_loop() {}
public:
void subdivide(Polyhedron &input,
Polyhedron &output)
{
// subdivide, then smooth
CModifierLoop<HalfedgeDS,Polyhedron,Kernel> builder(&input);
output.delegate(builder);
}
};
#endif // PWS_LOOP_SUBDIVISION