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/examples/Mesh_3/mesh_determinism_example.cpp

@@ -1,5 +1,7 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
+std::size_t TS;
+
 #include <CGAL/Mesh_triangulation_3.h>
 #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 #include <CGAL/Mesh_criteria_3.h>
@@ -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;

Mesh_3/examples/Mesh_3/mesh_polyhedral_domain.cpp

@@ -1,5 +1,10 @@
+
+
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
+std::size_t TS;
+
 #include <CGAL/Mesh_triangulation_3.h>
 #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 #include <CGAL/Mesh_criteria_3.h>
@@ -26,8 +31,10 @@ typedef CGAL::Mesh_criteria_3<Tr> 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<C3t3>(domain, criteria, no_perturb(), no_exude());
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(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();
+  }
+}

Mesh_3/examples/Mesh_3/mesh_polyhedral_domain_with_features.cpp

@@ -1,5 +1,6 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
+std::size_t TS;
 #include <CGAL/Mesh_triangulation_3.h>
 #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 #include <CGAL/Mesh_criteria_3.h>
@@ -22,8 +23,10 @@ typedef CGAL::Mesh_criteria_3<Tr> 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<C3t3>(domain, criteria);
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(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();
+  }
 }

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
 

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*>      Facet_handle_set;
-  typedef std::set<Halfedge*>   He_handle_set;
+  typedef std::set<Facet_handle>      Facet_handle_set;
+  typedef std::set<Halfedge_handle>   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 <typename P_>
 void
 Detect_features_in_polyhedra<P_>::
-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 ( &current_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(&current_he) == 0 )
+        if ( he_to_explore.erase(current_he) == 0 )
         {
-          he_to_explore.insert(&*(current_he.opposite()));
+          he_to_explore.insert(current_he->opposite());
         }
       }
     }

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 <typename Handle>
   bool operator()(const Handle& va, const Handle& vb) const {
-    return &*va < &*vb;
+    return va->ts < vb->ts;
   }
 };
 

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)

Mesh_3/include/CGAL/Mesh_polyhedron_3.h

@@ -44,6 +44,10 @@ private:
   typedef CGAL::HalfedgeDS_vertex_base<Refs, Tag, Point> 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<Refs,Tprev,Tvertex,Tface>
 {
 private:
   bool feature_edge;
+
+public:
+  std::size_t ts;
+
 public:
   
   Mesh_polyhedron_halfedge() 
@@ -95,6 +103,10 @@ public CGAL::HalfedgeDS_face_base<Refs,T_,Pln_>
 {
 private:
   Patch_id_ patch_id_;
+
+public:
+  std::size_t ts;
+
 public:
 
   typedef Patch_id_ Patch_id;

Mesh_3/include/CGAL/Mesh_vertex_base_3.h

@@ -158,6 +158,10 @@ public:
       Get_io_signature<int>()() + "+" +
       Get_io_signature<Index>()();
   }
+
+public:
+  std::size_t ts;
+
 private:
 
   int number_of_incident_facets_;

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<Point_3> random_point(1.);
 
+  Random_points_on_sphere_3<Point_3> random_point(1.);
   int i = n;
 #ifdef CGAL_MESH_3_VERBOSE
   std::cerr << "construct initial points:" << std::endl;

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<GT_,P_,TA_,Tag_,E_tag_>::
+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<GT_,P_,TA_,Tag_,E_tag_>::
 detect_features(FT angle_in_degree, Polyhedron& p)
 {
+  initialize_ts(p);
   // Get sharp features
   Mesh_3::detect_features(p,angle_in_degree);
   

Mesh_3/include/CGAL/perturb_mesh_3.h

@@ -77,9 +77,9 @@ default_perturbation_vector(const C3T3&,
   
   std::vector<Perturbation*> 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;
 }

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);
     }
 

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 <class Vb, class Cb >
+bool 
+  operator< (typename Triangulation_data_structure_3<Vb,Cb>::Cell_handle a, typename Triangulation_data_structure_3<Vb,Cb>::Cell_handle b)
+  {
+    std::cerr << "operator<" << std::endl;
+    return a->ts < b->ts;
+}
+
+
 } //namespace CGAL
 
 #endif // CGAL_TRIANGULATION_DATA_STRUCTURE_3_H

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 >

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 >