From 83546c53b102d470bf95dea4a261d1b95d1d3389 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 13 Feb 2014 13:11:49 +0100 Subject: [PATCH] Add a time stamp in vertices and cells of triangulations and polyhedra and use this for operator< of handles in order to make mesh generation deterministic --- .../Mesh_3/mesh_determinism_example.cpp | 7 +++- .../Mesh_3/mesh_polyhedral_domain.cpp | 24 +++++++++-- .../mesh_polyhedral_domain_with_features.cpp | 20 ++++++++- .../include/CGAL/Compact_mesh_cell_base_3.h | 5 ++- .../Mesh_3/Detect_features_in_polyhedra.h | 39 +++++++++--------- .../Mesh_3/Detect_polylines_in_polyhedra.h | 2 +- .../include/CGAL/Mesh_3/vertex_perturbation.h | 2 +- Mesh_3/include/CGAL/Mesh_polyhedron_3.h | 12 ++++++ Mesh_3/include/CGAL/Mesh_vertex_base_3.h | 4 ++ .../include/CGAL/Polyhedral_mesh_domain_3.h | 3 +- .../Polyhedral_mesh_domain_with_features_3.h | 28 +++++++++++++ Mesh_3/include/CGAL/perturb_mesh_3.h | 6 +-- .../include/CGAL/Compact_container.h | 4 ++ .../CGAL/Triangulation_data_structure_3.h | 41 +++++++++++++++---- .../CGAL/Triangulation_ds_cell_base_3.h | 2 + .../CGAL/Triangulation_ds_vertex_base_3.h | 2 + 16 files changed, 159 insertions(+), 42 deletions(-) diff --git a/Mesh_3/examples/Mesh_3/mesh_determinism_example.cpp b/Mesh_3/examples/Mesh_3/mesh_determinism_example.cpp index a1b80dd372c..99c04567588 100644 --- a/Mesh_3/examples/Mesh_3/mesh_determinism_example.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_determinism_example.cpp @@ -1,5 +1,7 @@ #include +std::size_t TS; + #include #include #include @@ -51,7 +53,7 @@ int main(int argc, char* argv[]) // Domain CGAL::Image_3 image; - image.read("data/liver.inr.gz"); + image.read("data/liver.inr"); Mesh_domain domain(image); // Mesh criteria @@ -63,6 +65,8 @@ int main(int argc, char* argv[]) for(std::size_t i = 0; i < nb_runs; ++i) { + TS = 0; + CGAL::default_random = CGAL::Random(0); std::cout << "------- Iteration " << (i+1) << " -------" << std::endl; @@ -106,6 +110,7 @@ int main(int argc, char* argv[]) std::ofstream medit_file4(num_str + std::string("out4-exude.mesh")); c3t3.output_to_medit(medit_file4); } + std::cerr << "TS = " << TS << std::endl; } return 0; diff --git a/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain.cpp b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain.cpp index 1cc6a6c5acf..510d227d04e 100644 --- a/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain.cpp @@ -1,5 +1,10 @@ + + + #include +std::size_t TS; + #include #include #include @@ -26,8 +31,10 @@ typedef CGAL::Mesh_criteria_3 Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; -int main() +int fct() { + CGAL::default_random = CGAL::Random(); + TS = 0; // Create input polyhedron Polyhedron polyhedron; std::ifstream input("data/elephant.off"); @@ -41,13 +48,14 @@ int main() cell_radius_edge_ratio=3); // Mesh generation - C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, no_perturb(), no_exude()); + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + // Output std::ofstream medit_file("out_1.mesh"); c3t3.output_to_medit(medit_file); medit_file.close(); - +#if 0 // Set tetrahedron size (keep cell_radius_edge_ratio), ignore facets Mesh_criteria new_criteria(cell_radius_edge_ratio=3, cell_size=0.03); @@ -57,6 +65,16 @@ int main() // Output medit_file.open("out_2.mesh"); c3t3.output_to_medit(medit_file); +#endif + std::cerr << "TS = " << TS << std::endl; return 0; } + + +int main() +{ + for(int i = 0; i < 3; i++){ + fct(); + } +} diff --git a/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp index 5d30be19e2d..7d94a7dfd46 100644 --- a/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp +++ b/Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp @@ -1,5 +1,6 @@ #include +std::size_t TS; #include #include #include @@ -22,8 +23,10 @@ typedef CGAL::Mesh_criteria_3 Mesh_criteria; // To avoid verbose function and named parameters call using namespace CGAL::parameters; -int main() +void fct() { + CGAL::default_random = CGAL::Random(); + TS = 0; // Create domain Mesh_domain domain("data/fandisk.off"); @@ -36,9 +39,22 @@ int main() cell_radius_edge_ratio = 3, cell_size = 0.05); // Mesh generation - C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, + odt(max_iteration_number =1), + no_perturb(), + no_exude()); // Output std::ofstream medit_file("out.mesh"); c3t3.output_to_medit(medit_file); + + + std::cerr << "TS = " << TS << std::endl; +} + +int main() +{ + for(int i = 0; i < 3; i++){ + fct(); + } } diff --git a/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h b/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h index b6bea2b9c40..2af1c3dad17 100644 --- a/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h +++ b/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h @@ -94,8 +94,7 @@ public: , subdomain_index_() , bits_(0) , sliver_cache_validity_(false) - { - } + {} Compact_mesh_cell_base_3(const Compact_mesh_cell_base_3& rhs) : circumcenter_(NULL) @@ -499,6 +498,8 @@ private: char bits_; mutable bool sliver_cache_validity_; +public: + std::size_t ts; }; // end class Compact_mesh_cell_base_3 diff --git a/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_polyhedra.h b/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_polyhedra.h index 3d496996379..1a22188d1ff 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_polyhedra.h +++ b/Mesh_3/include/CGAL/Mesh_3/Detect_features_in_polyhedra.h @@ -54,8 +54,8 @@ public: typedef typename Polyhedron::Halfedge Halfedge; typedef typename Polyhedron::Facet Facet; - typedef std::set Facet_handle_set; - typedef std::set He_handle_set; + typedef std::set Facet_handle_set; + typedef std::set He_handle_set; public: Detect_features_in_polyhedra() : current_surface_index_(1) {} @@ -68,7 +68,7 @@ public: private: Vector_3 facet_normal(const Facet_handle& f) const; bool is_sharp(const Halfedge_handle& he, FT cos_angle) const; - void flood(Facet& f, const int index, + void flood(Facet_handle f, const int index, Facet_handle_set& unsorted_faces) const; private: @@ -118,16 +118,17 @@ detect_surface_patches(Polyhedron& polyhedron) for ( typename Polyhedron::Facet_iterator fit = polyhedron.facets_begin(), end = polyhedron.facets_end() ; fit != end ; ++fit ) { - unsorted_faces.insert(&*fit); + Facet_handle fh = fit; + unsorted_faces.insert(fh); } // Flood while ( ! unsorted_faces.empty() ) { - Facet& f = **(unsorted_faces.begin()); + Facet_handle f = *(unsorted_faces.begin()); unsorted_faces.erase(unsorted_faces.begin()); - f.set_patch_id(current_surface_index_); + f->set_patch_id(current_surface_index_); flood(f,current_surface_index_,unsorted_faces); ++current_surface_index_; } @@ -218,52 +219,52 @@ is_sharp(const Halfedge_handle& he, FT cos_angle) const template void Detect_features_in_polyhedra:: -flood(Facet& f, const int index, Facet_handle_set& unsorted_faces) const +flood(Facet_handle f, const int index, Facet_handle_set& unsorted_faces) const { typedef typename Facet::Halfedge_around_facet_circulator Facet_he_circ; - Facet_he_circ begin = f.facet_begin(); + Facet_he_circ begin = f->facet_begin(); Facet_he_circ done = begin; // Initialize he_to_explore with halfedges of the starting facet He_handle_set he_to_explore; CGAL_For_all(begin,done) { - he_to_explore.insert(&*(begin->opposite())); + he_to_explore.insert(begin->opposite()); } // While there is something to explore while ( ! he_to_explore.empty() ) { // Get next halfedge to explore - Halfedge& he = **(he_to_explore.begin()); + Halfedge_handle he = *(he_to_explore.begin()); he_to_explore.erase(he_to_explore.begin()); // If we don't go through a border of the patch - if ( ! he.is_feature_edge() && ! he.is_border() ) + if ( ! he->is_feature_edge() && ! he->is_border() ) { - Facet& explored_facet = *(he.facet()); + Facet_handle explored_facet = he->facet(); // Mark facet and delete it from unsorted - explored_facet.set_patch_id(index); - unsorted_faces.erase(&explored_facet); + explored_facet->set_patch_id(index); + unsorted_faces.erase(explored_facet); // Add/Remove facet's halfedge to/from explore list - Facet_he_circ he_begin = explored_facet.facet_begin(); + Facet_he_circ he_begin = explored_facet->facet_begin(); Facet_he_circ he_done = he_begin; CGAL_For_all(he_begin,he_done) { - Halfedge& current_he = *he_begin; + Halfedge_handle current_he = he_begin; // do not explore heh again - if ( ¤t_he == &he ) { continue; } + if ( current_he == he ) { continue; } // if current_he is not in to_explore set, add it, otherwise remove it // (because we just explore the facet he_begin is pointing to) - if ( he_to_explore.erase(¤t_he) == 0 ) + if ( he_to_explore.erase(current_he) == 0 ) { - he_to_explore.insert(&*(current_he.opposite())); + he_to_explore.insert(current_he->opposite()); } } } diff --git a/Mesh_3/include/CGAL/Mesh_3/Detect_polylines_in_polyhedra.h b/Mesh_3/include/CGAL/Mesh_3/Detect_polylines_in_polyhedra.h index dac41b320b7..20a0a5433b9 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Detect_polylines_in_polyhedra.h +++ b/Mesh_3/include/CGAL/Mesh_3/Detect_polylines_in_polyhedra.h @@ -31,7 +31,7 @@ namespace CGAL { namespace Mesh_3 { struct Detect_polyline_less { template bool operator()(const Handle& va, const Handle& vb) const { - return &*va < &*vb; + return va->ts < vb->ts; } }; diff --git a/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h b/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h index 5ed80d48f63..e34bb60afff 100644 --- a/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h +++ b/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h @@ -731,7 +731,7 @@ private: const Vertex_handle& v) const { CGAL_assertion(cell->has_vertex(v)); - + std::cout << cell->ts << " " << v->ts << std::endl; const int i = cell->index(v); // fixed vertices: (the ones with index != i) diff --git a/Mesh_3/include/CGAL/Mesh_polyhedron_3.h b/Mesh_3/include/CGAL/Mesh_polyhedron_3.h index cf0bfbe0d48..677cf1e813f 100644 --- a/Mesh_3/include/CGAL/Mesh_polyhedron_3.h +++ b/Mesh_3/include/CGAL/Mesh_polyhedron_3.h @@ -44,6 +44,10 @@ private: typedef CGAL::HalfedgeDS_vertex_base Pdv_base; Set_of_indices indices; + +public: + std::size_t ts; + public: int nb_of_feature_edges; @@ -74,6 +78,10 @@ public CGAL::HalfedgeDS_halfedge_base { private: bool feature_edge; + +public: + std::size_t ts; + public: Mesh_polyhedron_halfedge() @@ -95,6 +103,10 @@ public CGAL::HalfedgeDS_face_base { private: Patch_id_ patch_id_; + +public: + std::size_t ts; + public: typedef Patch_id_ Patch_id; diff --git a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h index c802d377811..6765e9c4721 100644 --- a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h +++ b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h @@ -158,6 +158,10 @@ public: Get_io_signature()() + "+" + Get_io_signature()(); } + +public: + std::size_t ts; + private: int number_of_incident_facets_; diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index dceddaca71c..fe4386f438a 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -637,9 +637,8 @@ Construct_initial_points::operator()(OutputIterator pts, const Point_3 center( FT( (bbox.xmin() + bbox.xmax()) / 2), FT( (bbox.ymin() + bbox.ymax()) / 2), FT( (bbox.zmin() + bbox.zmax()) / 2) ); - - Random_points_on_sphere_3 random_point(1.); + Random_points_on_sphere_3 random_point(1.); int i = n; #ifdef CGAL_MESH_3_VERBOSE std::cerr << "construct initial points:" << std::endl; diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_with_features_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_with_features_3.h index 65ebf82fa4a..e43a80bf9df 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_with_features_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_with_features_3.h @@ -108,6 +108,8 @@ public: ~Polyhedral_mesh_domain_with_features_3() {} /// Detect features + void initialize_ts(Polyhedron& p); + void detect_features(FT angle_in_degree, Polyhedron& p); void detect_features(FT angle_in_degree = FT(60)) { detect_features(angle_in_degree, polyhedron_); } @@ -147,12 +149,38 @@ Polyhedral_mesh_domain_with_features_3(const std::string& filename) } +template < typename GT_, typename P_, typename TA_, + typename Tag_, typename E_tag_> +void +Polyhedral_mesh_domain_with_features_3:: +initialize_ts(Polyhedron& p) +{ + for(typename Polyhedron::Vertex_iterator v = p.vertices_begin(), + end = p.vertices_end() ; v != end ; ++v) + { + v->ts= TS++; + } + + for ( typename Polyhedron::Facet_iterator fit = p.facets_begin(), + end = p.facets_end() ; fit != end ; ++fit ) + { + fit->ts= TS++; + } + for ( typename Polyhedron::Halfedge_iterator hit = p.halfedges_begin(), + end = p.halfedges_end() ; hit != end ; ++hit ) + { + hit->ts= TS++; + } +} + + template < typename GT_, typename P_, typename TA_, typename Tag_, typename E_tag_> void Polyhedral_mesh_domain_with_features_3:: detect_features(FT angle_in_degree, Polyhedron& p) { + initialize_ts(p); // Get sharp features Mesh_3::detect_features(p,angle_in_degree); diff --git a/Mesh_3/include/CGAL/perturb_mesh_3.h b/Mesh_3/include/CGAL/perturb_mesh_3.h index 120fce89c0b..b8e3da75efa 100644 --- a/Mesh_3/include/CGAL/perturb_mesh_3.h +++ b/Mesh_3/include/CGAL/perturb_mesh_3.h @@ -77,9 +77,9 @@ default_perturbation_vector(const C3T3&, std::vector perturbation_vect; perturbation_vect.push_back(new Sq_radius(40,0.05)); - perturbation_vect.push_back(new Volume(40,0.05)); - perturbation_vect.push_back(new Dihedral_angle(40,0.05)); - perturbation_vect.push_back(new Li_random(100,0.15)); + //perturbation_vect.push_back(new Volume(40,0.05)); + //perturbation_vect.push_back(new Dihedral_angle(40,0.05)); + //perturbation_vect.push_back(new Li_random(100,0.15)); return perturbation_vect; } diff --git a/STL_Extension/include/CGAL/Compact_container.h b/STL_Extension/include/CGAL/Compact_container.h index 766576ff23b..55618a0ba09 100644 --- a/STL_Extension/include/CGAL/Compact_container.h +++ b/STL_Extension/include/CGAL/Compact_container.h @@ -864,21 +864,25 @@ namespace internal { // For std::less... bool operator<(const CC_iterator& other) const { + return m_ptr.p->ts < other.m_ptr.p->ts; return (m_ptr.p < other.m_ptr.p); } bool operator>(const CC_iterator& other) const { + std::cerr << ">"<< std::endl; return (m_ptr.p > other.m_ptr.p); } bool operator<=(const CC_iterator& other) const { + std::cerr << "<="<< std::endl; return (m_ptr.p <= other.m_ptr.p); } bool operator>=(const CC_iterator& other) const { + std::cerr << ">="<< std::endl; return (m_ptr.p >= other.m_ptr.p); } diff --git a/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h b/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h index 29520c89443..46e3af3eae9 100644 --- a/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h @@ -156,7 +156,7 @@ public: } }; //#endif - + public: Triangulation_data_structure_3() @@ -207,12 +207,16 @@ public: Vertex_handle create_vertex(const Vertex &v) { - return vertices().insert(v); + Vertex_handle vh = vertices().insert(v); + vh->ts = TS++; + return vh; } Vertex_handle create_vertex() { - return vertices().emplace(); + Vertex_handle vh = vertices().emplace(); + vh->ts = TS++; + return vh; } Vertex_handle create_vertex(Vertex_handle v) @@ -222,12 +226,16 @@ public: Cell_handle create_cell(const Cell &c) { - return cells().insert(c); + Cell_handle ch = cells().insert(c); + ch->ts = TS++; + return ch; } Cell_handle create_cell() { - return cells().emplace(); + Cell_handle ch = cells().emplace(); + ch->ts = TS++; + return ch; } Cell_handle create_cell(Cell_handle c) @@ -238,7 +246,9 @@ public: Cell_handle create_cell(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3) { - return cells().emplace(v0, v1, v2, v3); + Cell_handle ch = cells().emplace(v0, v1, v2, v3); + ch->ts = TS++; + return ch; } Cell_handle create_cell(Vertex_handle v0, Vertex_handle v1, @@ -246,7 +256,9 @@ public: Cell_handle n0, Cell_handle n1, Cell_handle n2, Cell_handle n3) { - return cells().emplace(v0, v1, v2, v3, n0, n1, n2, n3); + Cell_handle ch = cells().emplace(v0, v1, v2, v3, n0, n1, n2, n3); + ch->ts = TS++; + return ch; } Cell_handle create_face() @@ -259,7 +271,9 @@ public: Vertex_handle v2) { CGAL_triangulation_precondition(dimension()<3); - return cells().emplace(v0, v1, v2, Vertex_handle()); + Cell_handle ch = cells().emplace(v0, v1, v2, Vertex_handle()); + ch->ts = TS++; + return ch; } // The following functions come from TDS_2. @@ -3642,6 +3656,17 @@ count_cells(size_type & i, bool verbose, int level) const return true; } + + +template +bool + operator< (typename Triangulation_data_structure_3::Cell_handle a, typename Triangulation_data_structure_3::Cell_handle b) + { + std::cerr << "operator<" << std::endl; + return a->ts < b->ts; +} + + } //namespace CGAL #endif // CGAL_TRIANGULATION_DATA_STRUCTURE_3_H diff --git a/Triangulation_3/include/CGAL/Triangulation_ds_cell_base_3.h b/Triangulation_3/include/CGAL/Triangulation_ds_cell_base_3.h index 0c59e0eab57..dc63441d5d0 100644 --- a/Triangulation_3/include/CGAL/Triangulation_ds_cell_base_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_ds_cell_base_3.h @@ -207,6 +207,8 @@ private: Cell_handle N[4]; Vertex_handle V[4]; TDS_data _tds_data; +public: + std::size_t ts; }; template < class TDS > diff --git a/Triangulation_3/include/CGAL/Triangulation_ds_vertex_base_3.h b/Triangulation_3/include/CGAL/Triangulation_ds_vertex_base_3.h index 27f62fb3c0d..483422aed4e 100644 --- a/Triangulation_3/include/CGAL/Triangulation_ds_vertex_base_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_ds_vertex_base_3.h @@ -64,6 +64,8 @@ public: private: Cell_handle _c; +public: + std::size_t ts; }; template < class TDS >