merge experimental-packages/Triangulation_2-unrecursive

.gitattributes

@@ -4297,6 +4297,9 @@ Testsuite/test/makefile2 -text
 Testsuite/test/run_testsuite_with_cmake -text
 Triangulation_2/benchmark/Triangulation_2/CDT_with_intersection_2.cpp -text
 Triangulation_2/benchmark/Triangulation_2/Delaunay_remove.cpp -text
+Triangulation_2/benchmark/Triangulation_2/create_test_directories -text
+Triangulation_2/benchmark/Triangulation_2/launch_tests -text
+Triangulation_2/benchmark/Triangulation_2/sum -text
 Triangulation_2/demo/Triangulation_2/qt3/help/cindex.html svneol=native#text/html
 Triangulation_2/demo/Triangulation_2/qt3/help/cinput_point_layer.gif -text svneol=unset#image/gif
 Triangulation_2/demo/Triangulation_2/qt3/help/conflict_zone.gif -text svneol=unset#image/gif

Triangulation_2/benchmark/Triangulation_2/Triangulation_benchmark_2.cpp

@@ -1,4 +1,3 @@
-
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 #define BENCH_CLASS_LOCAL
@@ -35,8 +34,22 @@ int main(int argc, char **argv)
   Delaunay delaunay;
   delaunay.insert(points.begin(), points.end());
   t.stop();
-
+  std::cerr << t.time() << " seconds (";
-  std::cerr << t.time() << " seconds" << std::endl;
+#ifdef CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP 
+  std::cerr << "Old implementation)\n";
+#elif not defined(BENCH_CLASS_LOCAL)
+  std::cerr << "Thread-local vector)\n";
+#elif defined (BENCH_STACK)
+  std::cerr << "Local stack)\n";
+#elif defined (BENCH_VECTOR)
+  std::cerr << "Local vector)\n";
+#elif defined (BENCH_SMALL_VECTOR)
+  std::cerr << "Local \"home-made\" small vector)\n";
+#elif defined (BENCH_ARRAY)
+  std::cerr << "Local fixed size array)\n";
+#else
+#  error
+#endif
             
  return 0;
 }

Triangulation_2/benchmark/Triangulation_2/create_test_directories

@@ -0,0 +1,2 @@
+#!/bin/sh
+for o in STACK VECTOR SMALL_VECTOR ARRAY; do mkdir build-$o; pushd build-$o; cmake ..; cmake -D CMAKE_CXX_FLAGS:STRING="-DBENCH_CLASS_LOCAL -DBENCH_$o -frounding-math -Wall" -D CMAKE_BUILD_TYPE:STRING=Release .; make VERBOSE=1 ; popd; done

Triangulation_2/benchmark/Triangulation_2/launch_tests

@@ -0,0 +1,6 @@
+#!/bin/zsh
+for n in 1 2 3 4 5 6; do
+  for f in build*/*(.x); do
+     $f < ~/trunk/Generator/benchmark/Generator/10000000.pts 
+  done
+done          

Triangulation_2/benchmark/Triangulation_2/sum

@@ -0,0 +1,15 @@
+#!/usr/bin/env perl
+
+use strict;
+
+my %sum;
+my $c;
+
+while(<>) {
+    /([0-9.]+) seconds \((.*)\)/ && { $sum{$2} += $1 }
+}
+
+foreach $c (sort { $sum{$a} cmp $sum{$b} } (keys(%sum)))
+{
+    print "$sum{$c} : $c\n"
+} 

Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h

@@ -268,6 +268,7 @@ public:
 		      Face_handle fh, 
 		      int i,
 		      std::pair<OutputItFaces,OutputItBoundaryEdges>  pit)  const {
+#ifdef CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
    Face_handle fn = fh->neighbor(i);
    
    if ( fh->is_constrained(i) || ! test_conflict(p,fn)) {
@@ -278,6 +279,25 @@ public:
      pit = propagate_conflicts(p,fn,ccw(j),pit);
      pit = propagate_conflicts(p,fn,cw(j), pit);
    }
+#else // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
+    std::stack<std::pair<Face_handle, int> > stack;
+    stack.push(std::make_pair(fh, i));
+    while(!stack.empty()) {
+      const Face_handle fh = stack.top().first;
+      const int i = stack.top().second;
+      stack.pop();
+      const Face_handle fn = fh->neighbor(i);
+      if ( fh->is_constrained(i) || ! test_conflict(p,fn)) {
+        *(pit.second)++ = Edge(fn, fn->index(fh));
+      } else {
+        *(pit.first)++ = fn;
+        int j = fn->index(fh);
+        stack.push(std::make_pair(fn, cw(j)));
+        stack.push(std::make_pair(fn,ccw(j)));
+      }
+    }
+#endif // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
+
    return pit;
  }
 
@@ -384,7 +404,7 @@ public:
 template < class Gt, class Tds, class Itag >
 bool 
 Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
-is_flipable(Face_handle f, int i, bool perturb) const
+is_flipable(Face_handle f, int i, bool perturb /* = true */) const
   // determines if edge (f,i) can be flipped 
 {
   Face_handle ni = f->neighbor(i); 
@@ -469,12 +489,28 @@ Constrained_Delaunay_triangulation_2<Gt,Tds,Itag>::
 propagating_flip(Face_handle& f,int i)
 // similar to the corresponding function in Delaunay_triangulation_2.h 
 { 
+#ifdef CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
   if (!is_flipable(f,i)) return;
   Face_handle ni = f->neighbor(i); 
   flip(f, i); // flip for constrained triangulations
   propagating_flip(f,i); 
   i = ni->index(f->vertex(i)); 
   propagating_flip(ni,i); 
+#else // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
+  const Vertex_handle v = f->vertex(i);
+  std::stack<Face_handle> stack;
+  stack.push(f);
+  while(!stack.empty()) {
+    Face_handle f = stack.top();
+    stack.pop();
+    const int i = f->index(v);
+    if (!is_flipable(f,i,true)) continue;
+    const Face_handle n = f->neighbor(i);
+    flip(f, i); // flip for constrained triangulations
+    stack.push(n);
+    stack.push(f);
+  }
+#endif // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
 } 
 
  template < class Gt, class Tds, class Itag > 

Triangulation_2/include/CGAL/Delaunay_triangulation_2.h

@@ -27,6 +27,8 @@
 
 #include <CGAL/Triangulation_2.h>
 #include <CGAL/iterator.h>
+#include <CGAL/Small_stack.h>
+#include <vector>
 
 #ifndef CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
 #include <CGAL/Spatial_sort_traits_adapter_2.h>
@@ -36,6 +38,10 @@
 #include <boost/mpl/and.hpp>
 #endif //CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
 
+#ifdef CGAL_HAS_THREADS
+#  include <boost/thread/tss.hpp>
+#endif
+
 namespace CGAL {
 
 template < class Gt, 
@@ -175,7 +181,23 @@ public:
 #endif
 
 private:
-  void propagating_flip(Face_handle& f,int i);
+
+#ifdef BENCH_CLASS_LOCAL
+#if defined( BENCH_STACK )
+ mutable std::stack<Face_handle> stack;
+#elif defined( BENCH_VECTOR )
+mutable  std::vector<Face_handle> stack;
+#elif defined( BENCH_SMALL_VECTOR )
+ mutable Small_stack<Face_handle> stack;
+#elif defined( BENCH_ARRAY )
+ mutable Face_handle stack[10];
+ mutable int TOP;
+#else
+#error "Choose a stack type"
+#endif
+#endif // #ifdef BENCH_CLASS_LOCAL
+
+  void propagating_flip(const Face_handle f,const int i);
   void remove_2D(Vertex_handle v );
 
 // auxilliary functions for remove
@@ -438,10 +460,11 @@ private:
   template <class OutputItFaces, class OutputItBoundaryEdges> 
   std::pair<OutputItFaces,OutputItBoundaryEdges>
   propagate_conflicts (const Point  &p,
-		       Face_handle fh, 
+		       const Face_handle fh, 
-		       int i,
+		       const int i,
 		       std::pair<OutputItFaces,OutputItBoundaryEdges>
 		       pit)  const {
+#ifdef CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
     Face_handle fn = fh->neighbor(i);
     if (! test_conflict(p,fn)) {
       *(pit.second)++ = Edge(fn, fn->index(fh));
@@ -451,6 +474,41 @@ private:
       pit = propagate_conflicts(p,fn,ccw(j),pit);
       pit = propagate_conflicts(p,fn,cw(j), pit);
     }
+#else // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
+
+#ifdef CGAL_HAS_THREADS  
+    static boost::thread_specific_ptr< std::vector<std::pair<Face_handle,int> > > stack_safe_ptr;
+  if (stack_safe_ptr.get() == NULL) {
+    stack_safe_ptr.reset(new std::vector<std::pair<Face_handle,int> >());
+  }
+  std::vector<std::pair<Face_handle,int> >& stack=* stack_safe_ptr.get();
+#else
+  static std::stack<std::pair<Face_handle, int> > stack;
+#endif
+
+    stack.push_back(std::make_pair(fh, i));
+#ifdef CGAL_PROFILE
+    std::size_t S = 0;
+#endif
+    while(!stack.empty()) {
+#ifdef CGAL_PROFILE
+      S = (std::max)(S, stack.size());
+#endif
+      const Face_handle fh = stack.back().first;
+      const int i = stack.back().second;
+      stack.pop_back();
+      const Face_handle fn = fh->neighbor(i);
+      if (! test_conflict(p,fn)) {
+        *(pit.second)++ = Edge(fn, fn->index(fh));
+      } else {
+        *(pit.first)++ = fn;
+        int j = fn->index(fh);
+        stack.push_back(std::make_pair(fn, cw(j)));
+        stack.push_back(std::make_pair(fn,ccw(j)));
+      }
+    }
+  CGAL_HISTOGRAM_PROFILER("propagate_conflicts stack size ", S);
+#endif // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATE_CONFLICTS
     return pit;
   }
 
@@ -876,18 +934,100 @@ restore_Delaunay(Vertex_handle v)
 template < class Gt, class Tds >
 void
 Delaunay_triangulation_2<Gt,Tds>::
-propagating_flip(Face_handle& f,int i)
+propagating_flip(const Face_handle f, const int i)
 {
+
+#ifdef CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
   Face_handle n = f->neighbor(i);
-      
   if ( ON_POSITIVE_SIDE != 
        side_of_oriented_circle(n,  f->vertex(i)->point(), true) ) {          
     return;
   }
   this->flip(f, i);
   propagating_flip(f,i);
-  i = n->index(f->vertex(i));
+  propagating_flip(n,n->index(f->vertex(i)));
-  propagating_flip(n,i);
+#else // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
+  
+#ifndef BENCH_CLASS_LOCAL  
+#ifdef CGAL_HAS_THREADS  
+  static boost::thread_specific_ptr< std::vector<Face_handle> > stack_safe_ptr;
+  if (stack_safe_ptr.get() == NULL) {
+    stack_safe_ptr.reset(new std::vector<Face_handle>());
+  }
+  std::vector<Face_handle>& stack=* stack_safe_ptr.get();
+#else
+  static std::vector<Face_handle> stack;
+#endif
+#endif
+
+  
+#if  defined ( BENCH_CLASS_LOCAL ) && defined (BENCH_ARRAY )
+  const Vertex_handle v = f->vertex(i);
+
+  const Point& p = v->point();
+  stack[0] = f;
+  TOP=0;
+  while(TOP>-1) {
+    const Face_handle& f = stack[TOP];
+
+    const int i = f->index(v);
+    const Face_handle& n = f->neighbor(i);
+      
+    if ( ON_POSITIVE_SIDE != 
+         side_of_oriented_circle(n, p, true) ) {
+      --TOP;      
+    } else {
+      this->flip(f, i);
+      ++TOP;
+      stack[TOP] = n;
+    }
+  }
+#endif
+  
+#if ( ! defined ( BENCH_CLASS_LOCAL ) ) ||   defined( BENCH_VECTOR ) || defined( BENCH_SMALL_VECTOR )
+  const Vertex_handle& v = f->vertex(i);
+  const Point& p = v->point();
+
+  stack.push_back(f);
+  while(!stack.empty()) {
+    const Face_handle& f = stack.back();
+
+    const int i = f->index(v);
+    const Face_handle& n = f->neighbor(i);
+      
+    if ( ON_POSITIVE_SIDE != 
+         side_of_oriented_circle(n, p, true) ) {
+      stack.pop_back();      
+    } else {
+      this->flip(f, i);
+      stack.push_back(n);
+    }
+  }
+#endif
+  
+
+#if  defined ( BENCH_CLASS_LOCAL ) &&  defined( BENCH_STACK )
+  const Vertex_handle& v = f->vertex(i);
+  const Point& p = v->point();
+
+  stack.push(f);
+  while(!stack.empty()) {
+    const Face_handle& f = stack.top();
+
+    const int i = f->index(v);
+    const Face_handle& n = f->neighbor(i);
+      
+    if ( ON_POSITIVE_SIDE != 
+         side_of_oriented_circle(n, p, true) ) {
+      stack.pop();      
+    } else {
+      this->flip(f, i);
+      stack.push(n);
+    }
+  }
+#endif
+  
+#endif // NO CGAL_TRIANGULATION_2_USE_OLD_PROPAGATING_FLIP
 }
 
 

Triangulation_2/include/CGAL/Small_stack.h

@@ -0,0 +1,56 @@
+
+#include <vector>
+
+namespace CGAL {
+  template <typename T, int N=10>
+class Small_stack {
+
+private:
+  std::vector<T> V;
+
+  int TOP;
+  
+  T ts[N];
+
+public:
+  
+  Small_stack()
+    : TOP(-1)
+  {}
+
+  void push_back(const T& t)
+  {
+    if(TOP < (N-1)){
+      ++TOP;
+      ts[TOP] = t;
+    } else {
+      V.push_back(t);
+    }
+  }
+
+  const T& back() const
+  {
+    if(TOP < N){
+      return ts[TOP];
+    } else {
+      return V.back();
+    }    
+  }
+
+  void pop_back()
+  {
+    if(TOP < N){
+      --TOP;
+    } else {
+      V.pop_back();
+    }    
+  }
+
+  bool empty() const
+  {
+    return TOP == -1;
+  }
+  
+};
+
+} // namespace