diff --git a/.gitattributes b/.gitattributes index 6fb24e2472e..8b1689770b9 100644 --- a/.gitattributes +++ b/.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 diff --git a/Triangulation_3/benchmark/Triangulation_3/simple_2.cpp b/Triangulation_3/benchmark/Triangulation_3/simple_2.cpp new file mode 100644 index 00000000000..5fdb34bd6f8 --- /dev/null +++ b/Triangulation_3/benchmark/Triangulation_3/simple_2.cpp @@ -0,0 +1,50 @@ +#include +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Delaunay_triangulation_3 DT; +typedef K::Point_3 Point_3; +typedef CGAL::Timer Timer; +typedef CGAL::Creator_uniform_3 Creator; + +int main(int,char** argv) +{ + int n=atoi(argv[1]); + std::vector points; + points.reserve( n ); + + CGAL::Random rng(0); + CGAL::Random_points_in_sphere_3 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; +} diff --git a/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h b/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h index d0ca2a9bdf4..3a47663531b 100644 --- a/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_data_structure_3.h @@ -133,7 +133,7 @@ public: typedef Triple Edge; typedef Triangulation_simplex_3 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 typename Triangulation_data_structure_3::Cell_handle Triangulation_data_structure_3:: -create_star_3(Vertex_handle v, Cell_handle c, int li, - int prev_ind2) +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() ); @@ -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 typename Triangulation_data_structure_3::Cell_handle Triangulation_data_structure_3:: -create_star_3(Vertex_handle v, Cell_handle c, int li, - int prev_ind2) +recursive_create_star_3(Vertex_handle v, Cell_handle c, int li, + 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 +typename Triangulation_data_structure_3::Cell_handle +Triangulation_data_structure_3:: +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 - static boost::thread_specific_ptr< std::vector > stack_safe_ptr; - if (stack_safe_ptr.get() == NULL) { - stack_safe_ptr.reset(new std::vector()); - stack_safe_ptr.get()->reserve(64); - } - std::vector& adjacency_info_stack=* stack_safe_ptr.get(); -#else - static std::vector adjacency_info_stack; - adjacency_info_stack.reserve(64); -#endif + std::stack adjacency_info_stack; 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.pop_back(); + adjacency_info_stack.top().update_variables(zzz,cnew,ii,c,li,prev_ind2); + adjacency_info_stack.pop(); set_adjacency(nnn, zzz, cnew, ii); } }