make create_star_3 non recursive only after 100 recursive calls.

The speed observed on random data set produced by simple_2.cpp is equivalent
to the original recursive version and faster than the version better this commit

.gitattributes

@@ -4448,6 +4448,7 @@ Triangulation_2/examples/Triangulation_2/info_insert_with_zip_iterator_2.cpp -te
 Triangulation_2/examples/Triangulation_2/polygon_triangulation.cpp -text
 Triangulation_2/examples/Triangulation_2/print_cropped_voronoi.cpp -text
 Triangulation_2/test/Triangulation_2/test_delaunay_triangulation_proj.cpp -text
+Triangulation_3/benchmark/Triangulation_3/simple_2.cpp -text
 Triangulation_3/demo/Triangulation_3/CMakeLists.txt -text
 Triangulation_3/demo/Triangulation_3/MainWindow.cpp -text
 Triangulation_3/demo/Triangulation_3/MainWindow.h -text

Triangulation_3/benchmark/Triangulation_3/simple_2.cpp

@@ -0,0 +1,50 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/point_generators_3.h>
+
+#include <CGAL/Timer.h>
+#include <iostream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
+typedef CGAL::Delaunay_triangulation_3<K> DT;
+typedef K::Point_3 Point_3;
+typedef CGAL::Timer Timer;
+typedef CGAL::Creator_uniform_3<double,Point_3>  Creator;
+
+int main(int,char** argv)
+{
+  int n=atoi(argv[1]);
+  std::vector<Point_3> points;
+  points.reserve( n );
+
+  CGAL::Random rng(0);
+  CGAL::Random_points_in_sphere_3<Point_3,Creator> g( 1,rng);
+  CGAL::cpp11::copy_n( g, n, std::back_inserter(points));
+
+  Timer timer;
+  timer.start();
+  DT dt;
+  dt.insert(points.begin(), points.end());
+  timer.stop();
+  std::size_t N = dt.number_of_vertices();
+  
+  std::cerr << N << std::endl;
+  std::cout << timer.time() << " seconds (last"
+#ifdef CGAL_TDS_USE_RECURSIVE_CREATE_STAR_3
+  "CGAL_TDS_USE_RECURSIVE_CREATE_STAR_3"
+#endif
+#ifdef ONLY_NON_RECURSIVE
+  "ONLY_NON_RECURSIVE"
+#endif
+#ifdef RECURSIVE_FOLLOWED_BY_UNRECURSIVE
+  "RECURSIVE_FOLLOWED_BY_UNRECURSIVE"
+#endif
+#ifdef UNRECURSIVE_WITH_LOCAL_STACK
+  "UNRECURSIVE_WITH_LOCAL_STACK"
+#endif
+#ifdef RECURSIVE_FOLLOWED_BY_UNRECURSIVE_WITH_LOCAL_STACK
+  "RECURSIVE_FOLLOWED_BY_UNRECURSIVE_WITH_LOCAL_STACK"
+#endif
+  << ")\n";
+  return 0;
+}

Triangulation_3/include/CGAL/Triangulation_data_structure_3.h

@@ -133,7 +133,7 @@ public:
   typedef Triple<Cell_handle, int, int>            Edge;
 
   typedef Triangulation_simplex_3<Tds>             Simplex;
-#ifndef CGAL_TDS_USE_OLD_CREATE_STAR_3
+//#ifndef CGAL_TDS_USE_RECURSIVE_CREATE_STAR_3
   //internally used for create_star_3 (faster than a tuple)
   struct iAdjacency_info{
     int v1;
@@ -155,7 +155,7 @@ public:
       a6=v6;
     }
   };
-#endif  
+//#endif  
  
 
 public:
@@ -392,8 +392,19 @@ private:
 		   int i2, int j2, int next2,
 		   Vertex_handle v3);
 
+  #ifdef CGAL_TDS_USE_RECURSIVE_CREATE_STAR_3
   Cell_handle create_star_3(Vertex_handle v, Cell_handle c,
 	                    int li, int prev_ind2 = -1);
+  #else
+  Cell_handle recursive_create_star_3(Vertex_handle v, Cell_handle c, int li, int prev_ind2,int depth);
+  Cell_handle non_recursive_create_star_3(Vertex_handle v, Cell_handle c, int li, int prev_ind2);
+
+  Cell_handle create_star_3(Vertex_handle v, Cell_handle c,
+	                    int li, int prev_ind2 = -1)
+  {
+    return recursive_create_star_3(v,c,li,prev_ind2,0);
+  }
+  #endif
 
   Cell_handle create_star_2(Vertex_handle v,
 	                    Cell_handle c, int li);
@@ -1062,12 +1073,11 @@ private:
   // counts AND checks the validity
 };
 
-#ifdef CGAL_TDS_USE_OLD_CREATE_STAR_3
+#ifdef CGAL_TDS_USE_RECURSIVE_CREATE_STAR_3
 template <class Vb, class Cb >
 typename Triangulation_data_structure_3<Vb,Cb>::Cell_handle
 Triangulation_data_structure_3<Vb,Cb>::
-create_star_3(Vertex_handle v, Cell_handle c, int li,
+create_star_3(Vertex_handle v, Cell_handle c, int li, int prev_ind2)
-	      int prev_ind2)
 {
     CGAL_triangulation_precondition( dimension() == 3);
     CGAL_triangulation_precondition( c->tds_data().is_in_conflict() );
@@ -1121,20 +1131,71 @@ create_star_3(Vertex_handle v, Cell_handle c, int li,
     return cnew;
 }
 #else
-//This function used to be recursive.
-//This design resulted in call stack overflow in some cases.
-//When we rewrote it, to emulate the class stack, we use a c++ stack.
-//The stack used is a enriched std::vector that call reserve(N) when constructed
-//(faster that testing capacity>N at each call of create_star_3).
-//We use the class iAdjacency_info instead of a tuple because
-//at the moment we made the change it was faster like this.
-//The stack is a static variable of the function because it was also
-//faster when we tested it.
 template <class Vb, class Cb >
 typename Triangulation_data_structure_3<Vb,Cb>::Cell_handle
 Triangulation_data_structure_3<Vb,Cb>::
-create_star_3(Vertex_handle v, Cell_handle c, int li,
+recursive_create_star_3(Vertex_handle v, Cell_handle c, int li,
-	      int prev_ind2)
+                        int prev_ind2, int depth)
+{
+    if ( depth == 100 ) return non_recursive_create_star_3(v,c,li,prev_ind2);
+    CGAL_triangulation_precondition( dimension() == 3);
+    CGAL_triangulation_precondition( c->tds_data().is_in_conflict() );
+    CGAL_triangulation_precondition( ! c->neighbor(li)->tds_data().is_in_conflict() );
+
+    Cell_handle cnew = create_cell(c->vertex(0),
+				   c->vertex(1),
+				   c->vertex(2),
+				   c->vertex(3));
+    cnew->set_vertex(li, v);
+    Cell_handle c_li = c->neighbor(li);
+    set_adjacency(cnew, li, c_li, c_li->index(c));
+
+    // Look for the other neighbors of cnew.
+    for (int ii=0; ii<4; ++ii) {
+      if (ii == prev_ind2 || cnew->neighbor(ii) != Cell_handle())
+	  continue;
+      cnew->vertex(ii)->set_cell(cnew);
+
+      // Indices of the vertices of cnew such that ii,vj1,vj2,li positive.
+      Vertex_handle vj1 = c->vertex(next_around_edge(ii, li));
+      Vertex_handle vj2 = c->vertex(next_around_edge(li, ii));
+      Cell_handle cur = c;
+      int zz = ii;
+      Cell_handle n = cur->neighbor(zz);
+      // turn around the oriented edge vj1 vj2
+      while ( n->tds_data().is_in_conflict() ) {
+	CGAL_triangulation_assertion( n != c );
+	cur = n;
+	zz = next_around_edge(n->index(vj1), n->index(vj2));
+	n = cur->neighbor(zz);
+      }
+      // Now n is outside region, cur is inside.
+
+      n->tds_data().clear(); // Reset the flag for boundary cells.
+
+      int jj1 = n->index(vj1);
+      int jj2 = n->index(vj2);
+      Vertex_handle vvv = n->vertex(next_around_edge(jj1, jj2));
+      Cell_handle nnn = n->neighbor(next_around_edge(jj2, jj1));
+      int zzz = nnn->index(vvv);
+      if (nnn == cur) {
+	// Neighbor relation is reciprocal, ie
+	// the cell we are looking for is not yet created.
+	nnn = recursive_create_star_3(v, nnn, zz, zzz,depth+1);
+      }
+
+      set_adjacency(nnn, zzz, cnew, ii);
+    }
+
+    return cnew;
+}
+
+//We use the class iAdjacency_info instead of a tuple because
+//at the moment we made the change it was faster like this.
+template <class Vb, class Cb >
+typename Triangulation_data_structure_3<Vb,Cb>::Cell_handle
+Triangulation_data_structure_3<Vb,Cb>::
+non_recursive_create_star_3(Vertex_handle v, Cell_handle c, int li, int prev_ind2)
 {
     CGAL_triangulation_precondition( dimension() == 3);
     CGAL_triangulation_precondition( c->tds_data().is_in_conflict() );
@@ -1148,17 +1209,7 @@ create_star_3(Vertex_handle v, Cell_handle c, int li,
     Cell_handle c_li = c->neighbor(li);
     set_adjacency(cnew, li, c_li, c_li->index(c));
     
-#ifdef CGAL_HAS_THREADS  
+    std::stack<iAdjacency_info> adjacency_info_stack;
-  static boost::thread_specific_ptr< std::vector<iAdjacency_info> > stack_safe_ptr;
-    if (stack_safe_ptr.get() == NULL) {
-      stack_safe_ptr.reset(new std::vector<iAdjacency_info>());
-      stack_safe_ptr.get()->reserve(64);
-    }
-  std::vector<iAdjacency_info>& adjacency_info_stack=* stack_safe_ptr.get();  
-#else
-    static std::vector<iAdjacency_info> adjacency_info_stack;
-    adjacency_info_stack.reserve(64);
-#endif  
   
     int ii=0;
     do
@@ -1193,7 +1244,7 @@ create_star_3(Vertex_handle v, Cell_handle c, int li,
           // Neighbor relation is reciprocal, ie
           // the cell we are looking for is not yet created.
           //re-run the loop
-          adjacency_info_stack.push_back( iAdjacency_info(zzz,cnew,ii,c,li,prev_ind2) );
+          adjacency_info_stack.push( iAdjacency_info(zzz,cnew,ii,c,li,prev_ind2) );
           c=nnn;  li=zz; prev_ind2=zzz; ii=0;
           //copy-pasted from beginning to avoid if ii==0
           CGAL_triangulation_precondition( c->tds_data().is_in_conflict() );
@@ -1210,8 +1261,8 @@ create_star_3(Vertex_handle v, Cell_handle c, int li,
         if ( adjacency_info_stack.empty() ) return cnew;
         Cell_handle nnn=cnew; 
         int zzz;
-        adjacency_info_stack.back().update_variables(zzz,cnew,ii,c,li,prev_ind2);
+        adjacency_info_stack.top().update_variables(zzz,cnew,ii,c,li,prev_ind2);
-        adjacency_info_stack.pop_back();
+        adjacency_info_stack.pop();
         set_adjacency(nnn, zzz, cnew, ii);
       }
     }