diff --git a/Surface_mesher/include/CGAL/builder.h b/Surface_mesher/include/CGAL/builder.h new file mode 100644 index 00000000000..5690874eade --- /dev/null +++ b/Surface_mesher/include/CGAL/builder.h @@ -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 +// 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 +#else // CGAL_USE_POLYHEDRON_DESIGN_ONE // +#define CGAL_USE_POLYHEDRON_DESIGN_TWO 1 + +#include +#include +#include +#include +#include +#include + +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 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 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::" + << std::endl; + verr << "add_vertex(): capacity error: more than " << new_vertices + << " vertices added." << std::endl; + m_error = true; + return Vertex_handle(); + } + HalfedgeDS_decorator 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::" + << 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 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 decorator; + Halfedge_handle h = get_vertex_to_edge_map(w1); + decorator.set_face_halfedge( current_face, h); + ++new_faces; + return h; + } + + template + 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 + 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 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::" + << 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" + "::" << 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" + "::" << 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::" + << 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 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::" << 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:: +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:: +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:: +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::" + << 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 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:: +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:: +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:: +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::\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:: +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 // diff --git a/Surface_mesher/include/CGAL/enriched_polyhedron.h b/Surface_mesher/include/CGAL/enriched_polyhedron.h new file mode 100644 index 00000000000..11b2866228c --- /dev/null +++ b/Surface_mesher/include/CGAL/enriched_polyhedron.h @@ -0,0 +1,898 @@ +// -*- tab-width: 2 -*- + +/////////////////////////////////////////////////////////////////////////// +// // +// Class: Enriched_polyhedron // +// // +/////////////////////////////////////////////////////////////////////////// + +#ifndef _POLYGON_MESH_ +#define _POLYGON_MESH_ + +// CGAL stuff +#include +#include +#include + +#include + +// a refined facet with a normal and a tag +template +class Enriched_facet : public CGAL::HalfedgeDS_face_base +{ + // 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 Enriched_halfedge : public CGAL::HalfedgeDS_halfedge_base +{ +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 Enriched_vertex : public CGAL::HalfedgeDS_vertex_base +{ + // 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(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 struct Vertex_wrapper + { + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Normal; + typedef Enriched_vertex Vertex; + }; + + // wrap face + template struct Face_wrapper + { + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Normal; + typedef Enriched_facet Face; + }; + + // wrap halfedge + template struct Halfedge_wrapper + { + typedef typename Traits::Vector_3 Normal; + typedef Enriched_halfedge Halfedge; + }; +}; + +static const double PI = 3.1415926535897932384626; + +// compute facet normal +struct Facet_normal // (functor) +{ + template 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 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 Enriched_polyhedron : public CGAL::Polyhedron_3 +{ +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 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 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_ diff --git a/Surface_mesher/include/CGAL/pws_loop_subdivision.h b/Surface_mesher/include/CGAL/pws_loop_subdivision.h new file mode 100644 index 00000000000..ff0ba9158af --- /dev/null +++ b/Surface_mesher/include/CGAL/pws_loop_subdivision.h @@ -0,0 +1,386 @@ +#ifndef PWS_LOOP_SUBDIVISION +#define PWS_LOOP_SUBDIVISION + +#include "enriched_polyhedron.h" +#include "builder.h" + +template +class CModifierLoop : public CGAL::Modifier_base +{ +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 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 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 CSubdivider_loop +{ + public: + typedef typename Polyhedron::HalfedgeDS HalfedgeDS; + CSubdivider_loop() {} + ~CSubdivider_loop() {} + + public: + void subdivide(Polyhedron &input, + Polyhedron &output) + { + // subdivide, then smooth + CModifierLoop builder(&input); + output.delegate(builder); + } +}; + + +#endif // PWS_LOOP_SUBDIVISION