From 5c9e18a2f637f9b93e4b32b5cf4bb6193a2002ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 1 Nov 2017 15:09:18 +0100 Subject: [PATCH] Improved readability of Regular/Delaunay_triangulation_3 The quest continues --- .../include/CGAL/Delaunay_triangulation_3.h | 1454 +++--- .../include/CGAL/Regular_triangulation_3.h | 4154 +++++++++-------- .../include/CGAL/Triangulation_3.h | 10 +- 3 files changed, 2864 insertions(+), 2754 deletions(-) diff --git a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h index a574313cf2b..9e71c0ec71e 100644 --- a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h @@ -26,7 +26,6 @@ #include - #include #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING @@ -88,17 +87,16 @@ template < class Gt, class Tds_, class Delaunay_triangulation_3 : public Triangulation_3 { - typedef Delaunay_triangulation_3 Self; + typedef Delaunay_triangulation_3 Self; + public: - typedef Triangulation_3 Tr_Base; + typedef Triangulation_3 Tr_Base; - typedef typename Tr_Base::Triangulation_data_structure - Triangulation_data_structure; - typedef Gt Geom_traits; - typedef Compact_location Location_policy; + typedef typename Tr_Base::Triangulation_data_structure Triangulation_data_structure; + typedef Gt Geom_traits; + typedef Compact_location Location_policy; - typedef typename Tr_Base::Lock_data_structure Lock_data_structure; + typedef typename Tr_Base::Lock_data_structure Lock_data_structure; typedef typename Gt::Point_3 Point; typedef typename Gt::Segment_3 Segment; @@ -133,14 +131,14 @@ public: typedef typename Tr_Base::All_cells_iterator All_cells_iterator; - typedef typename Tr_Base::size_type size_type; - typedef typename Tr_Base::Locate_type Locate_type; + typedef typename Tr_Base::size_type size_type; + typedef typename Tr_Base::Locate_type Locate_type; //Tag to distinguish Delaunay from regular triangulations - typedef Tag_false Weighted_tag; + typedef Tag_false Weighted_tag; // Tag to distinguish periodic triangulations from others - typedef Tag_false Periodic_tag; + typedef Tag_false Periodic_tag; #ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 using Tr_Base::cw; @@ -178,60 +176,50 @@ public: protected: - Oriented_side - side_of_oriented_sphere(const Point &p0, const Point &p1, const Point &p2, - const Point &p3, const Point &t, bool perturb = false) const; + Oriented_side side_of_oriented_sphere(const Point& p0, const Point& p1, + const Point& p2, const Point& p3, + const Point& t, bool perturb = false) const; - Bounded_side - coplanar_side_of_bounded_circle(const Point &p, const Point &q, - const Point &r, const Point &s, bool perturb = false) const; + Bounded_side coplanar_side_of_bounded_circle(const Point& p, const Point& q, + const Point& r, const Point& s, bool perturb = false) const; // for dual: - Point - construct_circumcenter(const Point &p, const Point &q, const Point &r) const + Point construct_circumcenter(const Point& p, const Point& q, const Point& r) const { - return geom_traits().construct_circumcenter_3_object()(p, q, r); + return geom_traits().construct_circumcenter_3_object()(p, q, r); } - Line - construct_equidistant_line(const Point &p1, const Point &p2, - const Point &p3) const + Line construct_equidistant_line(const Point& p1, const Point& p2, const Point& p3) const { - return geom_traits().construct_equidistant_line_3_object()(p1, p2, p3); + return geom_traits().construct_equidistant_line_3_object()(p1, p2, p3); } - Ray - construct_ray(const Point &p, const Line &l) const + Ray construct_ray(const Point& p, const Line& l) const { - return geom_traits().construct_ray_3_object()(p, l); + return geom_traits().construct_ray_3_object()(p, l); } - Object - construct_object(const Point &p) const + Object construct_object(const Point& p) const { - return geom_traits().construct_object_3_object()(p); + return geom_traits().construct_object_3_object()(p); } - Object - construct_object(const Segment &s) const + Object construct_object(const Segment& s) const { - return geom_traits().construct_object_3_object()(s); + return geom_traits().construct_object_3_object()(s); } - Object - construct_object(const Ray &r) const + Object construct_object(const Ray& r) const { - return geom_traits().construct_object_3_object()(r); + return geom_traits().construct_object_3_object()(r); } - bool - less_distance(const Point &p, const Point &q, const Point &r) const + bool less_distance(const Point& p, const Point& q, const Point& r) const { - return geom_traits().compare_distance_3_object()(p, q, r) == SMALLER; + return geom_traits().compare_distance_3_object()(p, q, r) == SMALLER; } public: - Delaunay_triangulation_3(const Gt& gt = Gt(), Lock_data_structure *lock_ds = NULL) : Tr_Base(gt, lock_ds) {} @@ -242,19 +230,18 @@ public: // Create a 3D triangulation from 4 points which must be well-oriented // AND non-coplanar - Delaunay_triangulation_3(const Point &p0, const Point &p1, - const Point &p2, const Point &p3, + Delaunay_triangulation_3(const Point& p0, const Point& p1, + const Point& p2, const Point& p3, const Gt& gt = Gt(), Lock_data_structure *lock_ds = NULL) : Tr_Base(p0, p1, p2, p3, gt, lock_ds) {} - // copy constructor duplicates vertices and cells - Delaunay_triangulation_3(const Delaunay_triangulation_3 & tr) + Delaunay_triangulation_3(const Delaunay_triangulation_3& tr) : Tr_Base(tr) { - CGAL_triangulation_postcondition( is_valid() ); + CGAL_triangulation_postcondition(is_valid()); } template < typename InputIterator > @@ -262,7 +249,7 @@ public: const Gt& gt = Gt(), Lock_data_structure *lock_ds = NULL) : Tr_Base(gt, lock_ds) { - insert(first, last); + insert(first, last); } template < typename InputIterator > @@ -271,68 +258,69 @@ public: const Gt& gt = Gt()) : Tr_Base(gt, lock_ds) { - insert(first, last); + insert(first, last); } private: #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - std::vector + std::vector add_temporary_points_on_far_sphere(const size_t num_points) { - std::vector far_sphere_vertices; + std::vector far_sphere_vertices; - const size_t MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS = 1000000; - if (num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) + const size_t MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS = 1000000; + if(num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) + { + // Add temporary vertices on a "far sphere" to reduce contention on + // the infinite vertex + + // Get bbox + const Bbox_3& bbox = *this->get_bbox(); + + // Compute radius for far sphere + const double& xdelta = bbox.xmax() - bbox.xmin(); + const double& ydelta = bbox.ymax() - bbox.ymin(); + const double& zdelta = bbox.zmax() - bbox.zmin(); + const double radius = 1.3 * 0.5 * std::sqrt(xdelta*xdelta + + ydelta*ydelta + + zdelta*zdelta); + + // WARNING - TODO: this code has to be fixed because Vector_3 is not + // required by the traits concept + const typename Gt::Vector_3 center(bbox.xmin() + 0.5*xdelta, + bbox.ymin() + 0.5*ydelta, + bbox.zmin() + 0.5*zdelta); + Random_points_on_sphere_3 random_point(radius); + const int NUM_PSEUDO_INFINITE_VERTICES = static_cast( + tbb::task_scheduler_init::default_num_threads() * 3.5); + std::vector points_on_far_sphere; + + points_on_far_sphere.reserve(NUM_PSEUDO_INFINITE_VERTICES); + far_sphere_vertices.reserve(NUM_PSEUDO_INFINITE_VERTICES); + + for(int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point) + points_on_far_sphere.push_back(*random_point + center); + + spatial_sort(points_on_far_sphere.begin(), + points_on_far_sphere.end(), + geom_traits()); + + typename std::vector::const_iterator it_p = points_on_far_sphere.begin(); + typename std::vector::const_iterator it_p_end = points_on_far_sphere.end(); + + Vertex_handle hint; + for(; it_p != it_p_end ; ++it_p) { - // Add temporary vertices on a "far sphere" to reduce contention on - // the infinite vertex - - // Get bbox - const Bbox_3 &bbox = *this->get_bbox(); - // Compute radius for far sphere - const double& xdelta = bbox.xmax() - bbox.xmin(); - const double& ydelta = bbox.ymax() - bbox.ymin(); - const double& zdelta = bbox.zmax() - bbox.zmin(); - const double radius = 1.3 * 0.5 * std::sqrt(xdelta*xdelta + - ydelta*ydelta + - zdelta*zdelta); - // WARNING - TODO: this code has to be fixed because Vector_3 is not - // required by the traits concept - const typename Gt::Vector_3 center( - bbox.xmin() + 0.5*xdelta, - bbox.ymin() + 0.5*ydelta, - bbox.zmin() + 0.5*zdelta); - Random_points_on_sphere_3 random_point(radius); - const int NUM_PSEUDO_INFINITE_VERTICES = static_cast( - tbb::task_scheduler_init::default_num_threads() * 3.5); - std::vector points_on_far_sphere; - - points_on_far_sphere.reserve(NUM_PSEUDO_INFINITE_VERTICES); - far_sphere_vertices.reserve(NUM_PSEUDO_INFINITE_VERTICES); - - for (int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point) - points_on_far_sphere.push_back(*random_point + center); - - spatial_sort(points_on_far_sphere.begin(), - points_on_far_sphere.end(), - geom_traits()); - - typename std::vector::const_iterator it_p = points_on_far_sphere.begin(); - typename std::vector::const_iterator it_p_end = points_on_far_sphere.end(); - - Vertex_handle hint; - for ( ; it_p != it_p_end ; ++it_p) - { - hint = insert(*it_p, hint); - far_sphere_vertices.push_back(hint); - } + hint = insert(*it_p, hint); + far_sphere_vertices.push_back(hint); } + } - return far_sphere_vertices; + return far_sphere_vertices; } void remove_temporary_points_on_far_sphere( - const std::vector & far_sphere_vertices) + const std::vector& far_sphere_vertices) { if(!far_sphere_vertices.empty()) { @@ -346,19 +334,16 @@ public: #ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, InputIterator last, - typename boost::enable_if< - boost::is_convertible< - typename std::iterator_traits::value_type, - Point - > - >::type* = NULL - ) + std::ptrdiff_t insert(InputIterator first, InputIterator last, + typename boost::enable_if< + boost::is_convertible< + typename std::iterator_traits::value_type, + Point + > + >::type* = NULL) #else template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, InputIterator last) + std::ptrdiff_t insert(InputIterator first, InputIterator last) #endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO { #ifdef CGAL_TRIANGULATION_3_PROFILING @@ -366,38 +351,36 @@ public: #endif size_type n = number_of_vertices(); - std::vector points (first, last); - spatial_sort (points.begin(), points.end(), geom_traits()); + std::vector points(first, last); + spatial_sort(points.begin(), points.end(), geom_traits()); // Parallel #ifdef CGAL_LINKED_WITH_TBB - if (this->is_parallel()) + if(this->is_parallel()) { size_t num_points = points.size(); Vertex_handle hint; - + #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE std::vector far_sphere_vertices = add_temporary_points_on_far_sphere(num_points); #endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - + size_t i = 0; // Insert "num_points_seq" points sequentially // (or more if dim < 3 after that) - size_t num_points_seq = (std::min)(num_points, (size_t)100); - while (dimension() < 3 || i < num_points_seq) + size_t num_points_seq = (std::min)(num_points, (size_t) 100); + while(dimension() < 3 || i < num_points_seq) { hint = insert(points[i], hint); ++i; } tbb::enumerable_thread_specific tls_hint(hint); - tbb::parallel_for( - tbb::blocked_range( i, num_points ), - Insert_point(*this, points, tls_hint) - ); - + tbb::parallel_for(tbb::blocked_range(i, num_points), + Insert_point(*this, points, tls_hint)); + #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE remove_temporary_points_on_far_sphere(far_sphere_vertices); #endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE @@ -407,9 +390,11 @@ public: #endif // CGAL_LINKED_WITH_TBB { Vertex_handle hint; - for (typename std::vector::const_iterator p = points.begin(), end = points.end(); - p != end; ++p) + for(typename std::vector::const_iterator p = points.begin(), + end = points.end(); p != end; ++p) + { hint = insert(*p, hint); + } } #ifdef CGAL_TRIANGULATION_3_PROFILING @@ -419,142 +404,137 @@ public: return number_of_vertices() - n; } - #ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO private: //top stands for tuple-or-pair template const Point& top_get_first(const std::pair& pair) const { return pair.first; } + template const Info& top_get_second(const std::pair& pair) const { return pair.second; } + template const Point& top_get_first(const boost::tuple& tuple) const { return boost::get<0>(tuple); } + template const Info& top_get_second(const boost::tuple& tuple) const { return boost::get<1>(tuple); } template - std::ptrdiff_t insert_with_info(InputIterator first,InputIterator last) + std::ptrdiff_t insert_with_info(InputIterator first, InputIterator last) { size_type n = number_of_vertices(); std::vector indices; std::vector points; std::vector infos; std::size_t index=0; - for (InputIterator it=first;it!=last;++it){ + for(InputIterator it=first;it!=last;++it) + { Tuple_or_pair value=*it; - points.push_back( top_get_first(value) ); - infos.push_back ( top_get_second(value) ); + points.push_back(top_get_first(value)); + infos.push_back(top_get_second(value)); indices.push_back(index++); } typedef typename Pointer_property_map::type Pmap; typedef Spatial_sort_traits_adapter_3 Search_traits; - spatial_sort(indices.begin(), - indices.end(), + spatial_sort(indices.begin(), indices.end(), Search_traits(make_property_map(points),geom_traits())); #ifdef CGAL_LINKED_WITH_TBB - if(this->is_parallel()){ - - size_t num_points = points.size(); - Vertex_handle hint; + if(this->is_parallel()) + { + size_t num_points = points.size(); + Vertex_handle hint; #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - std::vector far_sphere_vertices = + std::vector far_sphere_vertices = add_temporary_points_on_far_sphere(num_points); #endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - size_t i = 0; - // Insert "num_points_seq" points sequentially - // (or more if dim < 3 after that) - size_t num_points_seq = (std::min)(num_points, (size_t)100); - while (dimension() < 3 || i < num_points_seq) - { - hint = insert(points[indices[i]], hint); - if (hint != Vertex_handle()) hint->info() = infos[indices[i]]; - ++i; - } + size_t i = 0; + // Insert "num_points_seq" points sequentially + // (or more if dim < 3 after that) + size_t num_points_seq = (std::min)(num_points, (size_t)100); + while(dimension() < 3 || i < num_points_seq) + { + hint = insert(points[indices[i]], hint); + if(hint != Vertex_handle()) hint->info() = infos[indices[i]]; + ++i; + } - tbb::enumerable_thread_specific tls_hint(hint); - tbb::parallel_for( - tbb::blocked_range( i, num_points ), - Insert_point_with_info(*this, points, infos, indices, tls_hint) - ); + tbb::enumerable_thread_specific tls_hint(hint); + tbb::parallel_for(tbb::blocked_range(i, num_points), + Insert_point_with_info(*this, points, infos, indices, tls_hint)); #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE remove_temporary_points_on_far_sphere(far_sphere_vertices); #endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - - } - // Sequential - else + } + // Sequential + else #endif + { + Vertex_handle hint; + for(typename std::vector::const_iterator + it = indices.begin(), end = indices.end(); it != end; ++it) { - Vertex_handle hint; - for (typename std::vector::const_iterator - it = indices.begin(), end = indices.end(); - it != end; ++it) { - hint = insert(points[*it], hint); - if (hint != Vertex_handle()) hint->info() = infos[*it]; - } + hint = insert(points[*it], hint); + if(hint != Vertex_handle()) + hint->info() = infos[*it]; } + } return number_of_vertices() - n; } public: - template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, - InputIterator last, - typename boost::enable_if< - boost::is_convertible< - typename std::iterator_traits::value_type, - std::pair::type> - > >::type* =NULL - ) + std::ptrdiff_t insert(InputIterator first, InputIterator last, + typename boost::enable_if< + boost::is_convertible< + typename std::iterator_traits::value_type, + std::pair::type> + > >::type* =NULL) { return insert_with_info< std::pair::type> >(first,last); } template std::ptrdiff_t - insert( boost::zip_iterator< boost::tuple > first, + insert(boost::zip_iterator< boost::tuple > first, boost::zip_iterator< boost::tuple > last, typename boost::enable_if< boost::mpl::and_< boost::is_convertible< typename std::iterator_traits::value_type, Point >, boost::is_convertible< typename std::iterator_traits::value_type, typename internal::Info_check::type > > - >::type* =NULL - ) + >::type* =NULL) { - return insert_with_info< boost::tuple::type> >(first,last); + return insert_with_info< boost::tuple::type> >(first,last); } #endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - Vertex_handle insert(const Point & p, Vertex_handle hint, - bool *could_lock_zone = NULL) + Vertex_handle insert(const Point& p, Vertex_handle hint, bool *could_lock_zone = NULL) { return insert(p, hint == Vertex_handle() ? this->infinite_cell() : hint->cell(), could_lock_zone); } - Vertex_handle insert(const Point & p, Cell_handle start = Cell_handle(), + Vertex_handle insert(const Point& p, Cell_handle start = Cell_handle(), bool *could_lock_zone = NULL); - Vertex_handle insert(const Point & p, Locate_type lt, + Vertex_handle insert(const Point& p, Locate_type lt, Cell_handle c, int li, int, - bool *could_lock_zone = NULL); + bool *could_lock_zone = NULL); public: // internal methods - template - Vertex_handle insert_and_give_new_cells(const Point &p, + Vertex_handle insert_and_give_new_cells(const Point& p, OutputItCells fit, - Cell_handle start = Cell_handle() ); + Cell_handle start = Cell_handle()); template Vertex_handle insert_and_give_new_cells(const Point& p, @@ -568,9 +548,8 @@ public: // internal methods OutputItCells fit); public: - #ifndef CGAL_NO_DEPRECATED_CODE - CGAL_DEPRECATED Vertex_handle move_point(Vertex_handle v, const Point & p); + CGAL_DEPRECATED Vertex_handle move_point(Vertex_handle v, const Point& p); #endif template - find_conflicts(const Point &p, Cell_handle c, + find_conflicts(const Point& p, Cell_handle c, OutputIteratorBoundaryFacets bfit, OutputIteratorCells cit, OutputIteratorInternalFacets ifit, bool *could_lock_zone = NULL) const { - CGAL_triangulation_precondition(dimension() >= 2); + CGAL_triangulation_precondition(dimension() >= 2); - std::vector cells; - cells.reserve(32); - std::vector facets; - facets.reserve(64); + std::vector cells; + cells.reserve(32); + std::vector facets; + facets.reserve(64); - if (dimension() == 2) { - Conflict_tester_2 tester(p, this); - ifit = Tr_Base::find_conflicts - (c, tester, - make_triple(std::back_inserter(facets), - std::back_inserter(cells), - ifit), could_lock_zone).third; - } - else { - Conflict_tester_3 tester(p, this); - ifit = Tr_Base::find_conflicts - (c, tester, - make_triple(std::back_inserter(facets), - std::back_inserter(cells), - ifit), could_lock_zone).third; - } + if(dimension() == 2) + { + Conflict_tester_2 tester(p, this); + ifit = Tr_Base::find_conflicts(c, tester, + make_triple(std::back_inserter(facets), + std::back_inserter(cells), + ifit), could_lock_zone).third; + } + else + { + Conflict_tester_3 tester(p, this); + ifit = Tr_Base::find_conflicts(c, tester, + make_triple(std::back_inserter(facets), + std::back_inserter(cells), + ifit), could_lock_zone).third; + } - // Reset the conflict flag on the boundary. - for(typename std::vector::iterator fit=facets.begin(); - fit != facets.end(); ++fit) { - fit->first->neighbor(fit->second)->tds_data().clear(); - *bfit++ = *fit; - } + // Reset the conflict flag on the boundary. + for(typename std::vector::iterator fit=facets.begin(); + fit != facets.end(); ++fit) + { + fit->first->neighbor(fit->second)->tds_data().clear(); + *bfit++ = *fit; + } - // Reset the conflict flag in the conflict cells. - for(typename std::vector::iterator ccit=cells.begin(); - ccit != cells.end(); ++ccit) { - (*ccit)->tds_data().clear(); - *cit++ = *ccit; - } - return make_triple(bfit, cit, ifit); + // Reset the conflict flag in the conflict cells. + for(typename std::vector::iterator ccit=cells.begin(); + ccit != cells.end(); ++ccit) + { + (*ccit)->tds_data().clear(); + *cit++ = *ccit; + } + return make_triple(bfit, cit, ifit); } template std::pair - find_conflicts(const Point &p, Cell_handle c, + find_conflicts(const Point& p, Cell_handle c, OutputIteratorBoundaryFacets bfit, OutputIteratorCells cit, bool *could_lock_zone = NULL) const @@ -656,31 +637,36 @@ public: vertices_on_conflict_zone_boundary(const Point&p, Cell_handle c, OutputIterator res) const { - CGAL_triangulation_precondition(dimension() >= 2); + CGAL_triangulation_precondition(dimension() >= 2); - // Get the facets on the boundary of the hole. - std::vector facets; - find_conflicts(p, c, std::back_inserter(facets), - Emptyset_iterator(), Emptyset_iterator()); + // Get the facets on the boundary of the hole. + std::vector facets; + find_conflicts(p, c, std::back_inserter(facets), + Emptyset_iterator(), Emptyset_iterator()); - // Then extract uniquely the vertices. - std::set vertices; - if (dimension() == 3) { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.insert(i->first->vertex((i->second+1)&3)); - vertices.insert(i->first->vertex((i->second+2)&3)); - vertices.insert(i->first->vertex((i->second+3)&3)); - } - } else { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.insert(i->first->vertex(cw(i->second))); - vertices.insert(i->first->vertex(ccw(i->second))); - } + // Then extract uniquely the vertices. + std::set vertices; + if(dimension() == 3) + { + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.insert(i->first->vertex((i->second+1)&3)); + vertices.insert(i->first->vertex((i->second+2)&3)); + vertices.insert(i->first->vertex((i->second+3)&3)); } + } + else + { + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.insert(i->first->vertex(cw(i->second))); + vertices.insert(i->first->vertex(ccw(i->second))); + } + } - return std::copy(vertices.begin(), vertices.end(), res); + return std::copy(vertices.begin(), vertices.end(), res); } // REMOVE @@ -706,19 +692,17 @@ public: // Parallel #ifdef CGAL_LINKED_WITH_TBB - if (this->is_parallel()) + if(this->is_parallel()) { // TODO: avoid that by asking for random-access iterators? std::vector vertices(first, beyond); tbb::concurrent_vector vertices_to_remove_sequentially; - tbb::parallel_for( - tbb::blocked_range( 0, vertices.size()), - Remove_point(*this, vertices, vertices_to_remove_sequentially) - ); + tbb::parallel_for(tbb::blocked_range(0, vertices.size()), + Remove_point(*this, vertices, vertices_to_remove_sequentially)); // Do the rest sequentially - for ( typename tbb::concurrent_vector::const_iterator + for(typename tbb::concurrent_vector::const_iterator it = vertices_to_remove_sequentially.begin(), it_end = vertices_to_remove_sequentially.end() ; it != it_end @@ -731,7 +715,8 @@ public: else #endif // CGAL_LINKED_WITH_TBB { - while (first != beyond) { + while(first != beyond) + { remove(*first); ++first; } @@ -748,123 +733,111 @@ public: size_type remove_cluster(InputIterator first, InputIterator beyond) { Self tmp; - Vertex_remover remover (tmp); + Vertex_remover remover(tmp); return Tr_Base::remove(first, beyond, remover); } // MOVE - Vertex_handle move_if_no_collision(Vertex_handle v, const Point &p); + Vertex_handle move_if_no_collision(Vertex_handle v, const Point& p); + Vertex_handle move(Vertex_handle v, const Point& p); - Vertex_handle move(Vertex_handle v, const Point &p); - - // return new cells (internal) + // return new cells(internal) template Vertex_handle move_if_no_collision_and_give_new_cells(Vertex_handle v, - const Point &p, + const Point& p, OutputItCells fit); private: - Bounded_side side_of_sphere(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, - const Point &p, bool perturb) const; + const Point& p, bool perturb) const; public: - // Queries - Bounded_side - side_of_sphere(Cell_handle c, const Point & p, - bool perturb = false) const + Bounded_side side_of_sphere(Cell_handle c, const Point& p, bool perturb = false) const { - return side_of_sphere(c->vertex(0), c->vertex(1), - c->vertex(2), c->vertex(3), p, perturb); + return side_of_sphere(c->vertex(0), c->vertex(1), + c->vertex(2), c->vertex(3), p, perturb); } - Bounded_side - side_of_circle( const Facet & f, const Point & p, bool perturb = false) const + Bounded_side side_of_circle(const Facet& f, const Point& p, bool perturb = false) const { - return side_of_circle(f.first, f.second, p, perturb); + return side_of_circle(f.first, f.second, p, perturb); } - Bounded_side - side_of_circle( Cell_handle c, int i, const Point & p, - bool perturb = false) const; + Bounded_side side_of_circle(Cell_handle c, int i, const Point& p, bool perturb = false) const; - Vertex_handle - nearest_vertex_in_cell(const Point& p, Cell_handle c) const; - - Vertex_handle - nearest_vertex(const Point& p, Cell_handle c = Cell_handle()) const; + Vertex_handle nearest_vertex_in_cell(const Point& p, Cell_handle c) const; + Vertex_handle nearest_vertex(const Point& p, Cell_handle c = Cell_handle()) const; bool is_Gabriel(Cell_handle c, int i) const; bool is_Gabriel(Cell_handle c, int i, int j) const; bool is_Gabriel(const Facet& f)const ; bool is_Gabriel(const Edge& e) const; - bool is_delaunay_after_displacement(Vertex_handle v, - const Point &p) const; + bool is_delaunay_after_displacement(Vertex_handle v, const Point& p) const; // Dual functions Point dual(Cell_handle c) const; - - Object dual(const Facet & f) const - { return dual( f.first, f.second ); } - + Object dual(const Facet& f) const { return dual(f.first, f.second); } Object dual(Cell_handle c, int i) const; Line dual_support(Cell_handle c, int i) const; bool is_valid(bool verbose = false, int level = 0) const; - bool is_valid(Cell_handle c, bool verbose = false, int level = 0) const; template < class Stream> - Stream& draw_dual(Stream & os) - { - for (Finite_facets_iterator fit = finite_facets_begin(), - end = finite_facets_end(); - fit != end; ++fit) { - Object o = dual(*fit); - if (const Segment *s = object_cast(&o)) os << *s; - else if (const Ray *r = object_cast(&o)) os << *r; - else if (const Point *p = object_cast(&o)) os << *p; - } - return os; + Stream& draw_dual(Stream& os) + { + for(Finite_facets_iterator fit = finite_facets_begin(), + end = finite_facets_end(); + fit != end; ++fit) { + Object o = dual(*fit); + if(const Segment* s = object_cast(&o)) + os << *s; + else if(const Ray* r = object_cast(&o)) + os << *r; + else if(const Point* p = object_cast(&o)) + os << *p; } + return os; + } protected: Vertex_handle - nearest_vertex(const Point &p, Vertex_handle v, Vertex_handle w) const + nearest_vertex(const Point& p, Vertex_handle v, Vertex_handle w) const { - // In case of equality, v is returned. - CGAL_triangulation_precondition(v != w); + // In case of equality, v is returned. + CGAL_triangulation_precondition(v != w); - if (is_infinite(v)) - return w; - if (is_infinite(w)) - return v; - return less_distance(p, w->point(), v->point()) ? w : v; + if(is_infinite(v)) + return w; + if(is_infinite(w)) + return v; + return less_distance(p, w->point(), v->point()) ? w : v; } class Conflict_tester_3 { - const Point &p; + const Point& p; const Self *t; public: - - Conflict_tester_3(const Point &pt, const Self *tr) + Conflict_tester_3(const Point& pt, const Self *tr) : p(pt), t(tr) {} bool operator()(const Cell_handle c) const { return t->side_of_sphere(c, p, true) == ON_BOUNDED_SIDE; } - Oriented_side compare_weight(const Point &, const Point &) const + + Oriented_side compare_weight(const Point& , const Point& ) const { return ZERO; } + bool test_initial_cell(Cell_handle) const { return true; @@ -873,22 +846,23 @@ protected: class Conflict_tester_2 { - const Point &p; + const Point& p; const Self *t; public: - - Conflict_tester_2(const Point &pt, const Self *tr) + Conflict_tester_2(const Point& pt, const Self *tr) : p(pt), t(tr) {} bool operator()(const Cell_handle c) const { return t->side_of_circle(c, 3, p, true) == ON_BOUNDED_SIDE; } - Oriented_side compare_weight(const Point &, const Point &) const + + Oriented_side compare_weight(const Point& , const Point& ) const { return ZERO; } + bool test_initial_cell(Cell_handle) const { return true; @@ -902,12 +876,12 @@ protected: template void process_cells_in_conflict(InputIterator, InputIterator) const {} - void reinsert_vertices(Vertex_handle ) {} - Vertex_handle replace_vertex(Cell_handle c, int index, - const Point &) { + void reinsert_vertices(Vertex_handle) {} + Vertex_handle replace_vertex(Cell_handle c, int index, const Point& ) + { return c->vertex(index); } - void hide_point(Cell_handle, const Point &) {} + void hide_point(Cell_handle, const Point& ) {} }; #ifdef CGAL_LINKED_WITH_TBB @@ -918,46 +892,45 @@ protected: typedef typename DT::Point Point; typedef typename DT::Vertex_handle Vertex_handle; - DT & m_dt; - const std::vector & m_points; - tbb::enumerable_thread_specific & m_tls_hint; + DT& m_dt; + const std::vector& m_points; + tbb::enumerable_thread_specific& m_tls_hint; public: // Constructor - Insert_point(DT & dt, - const std::vector & points, - tbb::enumerable_thread_specific & tls_hint) + Insert_point(DT& dt, + const std::vector& points, + tbb::enumerable_thread_specific& tls_hint) : m_dt(dt), m_points(points), m_tls_hint(tls_hint) {} // Constructor - Insert_point(const Insert_point &ip) + Insert_point(const Insert_point& ip) : m_dt(ip.m_dt), m_points(ip.m_points), m_tls_hint(ip.m_tls_hint) {} // operator() - void operator()( const tbb::blocked_range& r ) const + void operator()(const tbb::blocked_range& r) const { #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING static Profile_branch_counter_3 bcounter( "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); #endif - Vertex_handle &hint = m_tls_hint.local(); - for( std::size_t i_point = r.begin() ; i_point != r.end() ; ++i_point) + Vertex_handle& hint = m_tls_hint.local(); + for(std::size_t i_point = r.begin() ; i_point != r.end() ; ++i_point) { bool success = false; while(!success) { - if (m_dt.try_lock_vertex(hint) && m_dt.try_lock_point(m_points[i_point])) + if(m_dt.try_lock_vertex(hint) && m_dt.try_lock_point(m_points[i_point])) { bool could_lock_zone; - Vertex_handle new_hint = m_dt.insert( - m_points[i_point], hint, &could_lock_zone); + Vertex_handle new_hint = m_dt.insert(m_points[i_point], hint, &could_lock_zone); m_dt.unlock_all_elements(); - if (could_lock_zone) + if(could_lock_zone) { hint = new_hint; success = true; @@ -994,57 +967,56 @@ protected: typedef typename DT::Vertex_handle Vertex_handle; typedef typename DT::Triangulation_data_structure::Vertex::Info Info; - DT & m_dt; - const std::vector & m_points; - const std::vector & m_infos; - const std::vector & m_indices; - tbb::enumerable_thread_specific & m_tls_hint; + DT& m_dt; + const std::vector& m_points; + const std::vector& m_infos; + const std::vector& m_indices; + tbb::enumerable_thread_specific& m_tls_hint; public: // Constructor - Insert_point_with_info(DT & dt, - const std::vector & points, - const std::vector & infos, - const std::vector & indices, - tbb::enumerable_thread_specific & tls_hint) - : m_dt(dt), m_points(points), m_infos(infos), m_indices(indices), + Insert_point_with_info(DT& dt, + const std::vector& points, + const std::vector& infos, + const std::vector& indices, + tbb::enumerable_thread_specific& tls_hint) + : m_dt(dt), m_points(points), m_infos(infos), m_indices(indices), m_tls_hint(tls_hint) {} // Constructor - Insert_point_with_info(const Insert_point_with_info &ip) + Insert_point_with_info(const Insert_point_with_info& ip) : m_dt(ip.m_dt), m_points(ip.m_points), m_infos(ip.m_infos), m_indices(ip.m_indices), m_tls_hint(ip.m_tls_hint) {} // operator() - void operator()( const tbb::blocked_range& r ) const + void operator()(const tbb::blocked_range& r) const { #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING static Profile_branch_counter_3 bcounter( "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert_with_info]"); #endif - Vertex_handle &hint = m_tls_hint.local(); - for( std::size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) + Vertex_handle& hint = m_tls_hint.local(); + for(std::size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) { bool success = false; std::ptrdiff_t i_point = m_indices[i_idx]; - const Point &p = m_points[i_point]; + const Point& p = m_points[i_point]; while(!success) { - if (m_dt.try_lock_vertex(hint) && m_dt.try_lock_point(p)) + if(m_dt.try_lock_vertex(hint) && m_dt.try_lock_point(p)) { bool could_lock_zone; - Vertex_handle new_hint = m_dt.insert( - p, hint, &could_lock_zone); + Vertex_handle new_hint = m_dt.insert(p, hint, &could_lock_zone); m_dt.unlock_all_elements(); - if (could_lock_zone) + if(could_lock_zone) { hint = new_hint; - if (hint!=Vertex_handle()) hint->info()=m_infos[i_point]; + if(hint!=Vertex_handle()) hint->info()=m_infos[i_point]; success = true; #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING ++bcounter; @@ -1078,30 +1050,29 @@ protected: typedef typename DT::Point Point; typedef typename DT::Vertex_handle Vertex_handle; - DT & m_dt; - const std::vector & m_vertices; - tbb::concurrent_vector & m_vertices_to_remove_sequentially; + DT& m_dt; + const std::vector& m_vertices; + tbb::concurrent_vector& m_vertices_to_remove_sequentially; public: // Constructor - Remove_point(DT & dt, - const std::vector & vertices, - tbb::concurrent_vector & - vertices_to_remove_sequentially) - : m_dt(dt), m_vertices(vertices), + Remove_point(DT& dt, + const std::vector& vertices, + tbb::concurrent_vector& vertices_to_remove_sequentially) + : m_dt(dt), m_vertices(vertices), m_vertices_to_remove_sequentially(vertices_to_remove_sequentially) {} // Constructor - Remove_point(const Remove_point &rp) + Remove_point(const Remove_point& rp) : m_dt(rp.m_dt), m_vertices(rp.m_vertices), m_vertices_to_remove_sequentially(rp.m_vertices_to_remove_sequentially) {} // operator() - void operator()( const tbb::blocked_range& r ) const + void operator()(const tbb::blocked_range& r) const { - for( size_t i_vertex = r.begin() ; i_vertex != r.end() ; ++i_vertex) + for(size_t i_vertex = r.begin() ; i_vertex != r.end() ; ++i_vertex) { Vertex_handle v = m_vertices[i_vertex]; bool could_lock_zone, needs_to_be_done_sequentially; @@ -1110,9 +1081,10 @@ protected: needs_to_be_done_sequentially = !m_dt.remove(v, &could_lock_zone); m_dt.unlock_all_elements(); - } while (!could_lock_zone); + } + while(!could_lock_zone); - if (needs_to_be_done_sequentially) + if(needs_to_be_done_sequentially) m_vertices_to_remove_sequentially.push_back(v); } } @@ -1134,16 +1106,16 @@ protected: template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -insert(const Point & p, Cell_handle start, bool *could_lock_zone) +insert(const Point& p, Cell_handle start, bool *could_lock_zone) { Locate_type lt; int li, lj; // Parallel - if (could_lock_zone) + if(could_lock_zone) { Cell_handle c = locate(p, lt, li, lj, start, could_lock_zone); - if (*could_lock_zone) + if(*could_lock_zone) return insert(p, lt, c, li, lj, could_lock_zone); else return Vertex_handle(); @@ -1154,33 +1126,32 @@ insert(const Point & p, Cell_handle start, bool *could_lock_zone) Cell_handle c = locate(p, lt, li, lj, start); return insert(p, lt, c, li, lj); } - } template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -insert(const Point & p, Locate_type lt, Cell_handle c, int li, int lj, - bool *could_lock_zone) +insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj, bool *could_lock_zone) { - switch (dimension()) { - case 3: + switch(dimension()) + { + case 3: { Conflict_tester_3 tester(p, this); Vertex_handle v = insert_in_conflict(p, lt, c, li, lj, tester, hidden_point_visitor, could_lock_zone); return v; }// dim 3 - case 2: + case 2: { Conflict_tester_2 tester(p, this); return insert_in_conflict(p, lt, c, li, lj, tester, hidden_point_visitor, could_lock_zone); }//dim 2 - default : - // dimension <= 1 - // Do not use the generic insert. - return Tr_Base::insert(p, c); + default : + // dimension <= 1 + // Do not use the generic insert. + return Tr_Base::insert(p, c); } } @@ -1188,21 +1159,24 @@ template < class Gt, class Tds, class Lds > template typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -insert_and_give_new_cells(const Point &p, - OutputItCells fit, - Cell_handle start) +insert_and_give_new_cells(const Point &p, OutputItCells fit, Cell_handle start) { Vertex_handle v = insert(p, start); int dimension = this->dimension(); - if(dimension == 3) this->incident_cells(v, fit); + if(dimension == 3) + { + this->incident_cells(v, fit); + } else if(dimension == 2) { Cell_handle c = v->cell(), end = c; - do { + do + { *fit++ = c; int i = c->index(v); c = c->neighbor((i+1)%3); - } while(c != end); + } + while(c != end); } else if(dimension == 1) { @@ -1210,7 +1184,10 @@ insert_and_give_new_cells(const Point &p, *fit++ = c; *fit++ = c->neighbor((~(c->index(v)))&1); } - else *fit++ = v->cell(); // dimension = 0 + else // dimension = 0 + { + *fit++ = v->cell(); + } return v; } @@ -1218,21 +1195,24 @@ template < class Gt, class Tds, class Lds > template typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -insert_and_give_new_cells(const Point& p, - OutputItCells fit, - Vertex_handle hint) +insert_and_give_new_cells(const Point& p, OutputItCells fit, Vertex_handle hint) { Vertex_handle v = insert(p, hint); int dimension = this->dimension(); - if(dimension == 3) this->incident_cells(v, fit); + if(dimension == 3) + { + this->incident_cells(v, fit); + } else if(dimension == 2) { Cell_handle c = v->cell(), end = c; - do { + do + { *fit++ = c; int i = c->index(v); c = c->neighbor((i+1)%3); - } while(c != end); + } + while(c != end); } else if(dimension == 1) { @@ -1240,7 +1220,10 @@ insert_and_give_new_cells(const Point& p, *fit++ = c; *fit++ = c->neighbor((~(c->index(v)))&1); } - else *fit++ = v->cell(); // dimension = 0 + else // dimension = 0 + { + *fit++ = v->cell(); + } return v; } @@ -1255,15 +1238,20 @@ insert_and_give_new_cells(const Point& p, { Vertex_handle v = insert(p, lt, c, li, lj); int dimension = this->dimension(); - if(dimension == 3) this->incident_cells(v, fit); + if(dimension == 3) + { + this->incident_cells(v, fit); + } else if(dimension == 2) { Cell_handle c = v->cell(), end = c; - do { + do + { *fit++ = c; int i = c->index(v); c = c->neighbor((i+1)%3); - } while(c != end); + } + while(c != end); } else if(dimension == 1) { @@ -1271,7 +1259,10 @@ insert_and_give_new_cells(const Point& p, *fit++ = c; *fit++ = c->neighbor((~(c->index(v)))&1); } - else *fit++ = v->cell(); // dimension = 0 + else // dimension = 0 + { + *fit++ = v->cell(); + } return v; } @@ -1279,31 +1270,34 @@ insert_and_give_new_cells(const Point& p, template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -move_point(Vertex_handle v, const Point & p) +move_point(Vertex_handle v, const Point& p) { - CGAL_triangulation_precondition(! is_infinite(v)); - CGAL_triangulation_expensive_precondition(is_vertex(v)); + CGAL_triangulation_precondition(! is_infinite(v)); + CGAL_triangulation_expensive_precondition(is_vertex(v)); - // Dummy implementation for a start. + // Dummy implementation for a start. - // Remember an incident vertex to restart - // the point location after the removal. - Cell_handle c = v->cell(); - Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0); - CGAL_triangulation_assertion(old_neighbor != v); + // Remember an incident vertex to restart + // the point location after the removal. + Cell_handle c = v->cell(); + Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0); + CGAL_triangulation_assertion(old_neighbor != v); - remove(v); + remove(v); - if (dimension() <= 0) - return insert(p); - return insert(p, old_neighbor->cell()); + if(dimension() <= 0) + return insert(p); + + return insert(p, old_neighbor->cell()); } #endif template template -class Delaunay_triangulation_3::Vertex_remover { +class Delaunay_triangulation_3::Vertex_remover +{ typedef DelaunayTriangulation_3 Delaunay; + public: typedef Nullptr_t Hidden_points_iterator; @@ -1315,16 +1309,20 @@ public: Hidden_points_iterator hidden_points_begin() { return NULL; } Hidden_points_iterator hidden_points_end() { return NULL; } - Bounded_side side_of_bounded_circle(const Point &p, const Point &q, - const Point &r, const Point &s, bool perturb = false) const { + Bounded_side side_of_bounded_circle(const Point& p, const Point& q, + const Point& r, const Point& s, + bool perturb = false) const + { return tmp.coplanar_side_of_bounded_circle(p,q,r,s,perturb); } }; template template -class Delaunay_triangulation_3::Vertex_inserter { +class Delaunay_triangulation_3::Vertex_inserter +{ typedef DelaunayTriangulation_3 Delaunay; + public: typedef Nullptr_t Hidden_points_iterator; @@ -1336,16 +1334,18 @@ public: Hidden_points_iterator hidden_points_begin() { return NULL; } Hidden_points_iterator hidden_points_end() { return NULL; } - Vertex_handle insert(const Point& p, - Locate_type lt, Cell_handle c, int li, int lj) { + Vertex_handle insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj) + { return tmp.insert(p, lt, c, li, lj); } - Vertex_handle insert(const Point& p, Cell_handle c) { + Vertex_handle insert(const Point& p, Cell_handle c) + { return tmp.insert(p, c); } - Vertex_handle insert(const Point& p) { + Vertex_handle insert(const Point& p) + { return tmp.insert(p); } }; @@ -1356,7 +1356,7 @@ Delaunay_triangulation_3:: remove(Vertex_handle v) { Self tmp; - Vertex_remover remover (tmp); + Vertex_remover remover(tmp); Tr_Base::remove(v, remover); CGAL_triangulation_expensive_postcondition(is_valid()); @@ -1368,7 +1368,7 @@ Delaunay_triangulation_3:: remove(Vertex_handle v, bool *could_lock_zone) { Self tmp; - Vertex_remover remover (tmp); + Vertex_remover remover(tmp); bool ret = Tr_Base::remove(v, remover, could_lock_zone); CGAL_triangulation_expensive_postcondition(is_valid()); @@ -1378,11 +1378,11 @@ remove(Vertex_handle v, bool *could_lock_zone) template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -move_if_no_collision(Vertex_handle v, const Point &p) +move_if_no_collision(Vertex_handle v, const Point& p) { Self tmp; - Vertex_remover remover (tmp); - Vertex_inserter inserter (*this); + Vertex_remover remover(tmp); + Vertex_inserter inserter(*this); Vertex_handle res = Tr_Base::move_if_no_collision(v,p,remover,inserter); CGAL_triangulation_expensive_postcondition(is_valid()); @@ -1392,12 +1392,15 @@ move_if_no_collision(Vertex_handle v, const Point &p) template typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -move(Vertex_handle v, const Point &p) { +move(Vertex_handle v, const Point& p) +{ CGAL_triangulation_precondition(!is_infinite(v)); - if(v->point() == p) return v; + if(v->point() == p) + return v; + Self tmp; - Vertex_remover remover (tmp); - Vertex_inserter inserter (*this); + Vertex_remover remover(tmp); + Vertex_inserter inserter(*this); return Tr_Base::move(v,p,remover,inserter); } @@ -1408,7 +1411,7 @@ Delaunay_triangulation_3:: remove_and_give_new_cells(Vertex_handle v, OutputItCells fit) { Self tmp; - Vertex_remover remover (tmp); + Vertex_remover remover(tmp); Tr_Base::remove_and_give_new_cells(v,remover,fit); CGAL_triangulation_expensive_postcondition(is_valid()); @@ -1418,12 +1421,12 @@ template < class Gt, class Tds, class Lds > template typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: -move_if_no_collision_and_give_new_cells(Vertex_handle v, const Point &p, - OutputItCells fit) +move_if_no_collision_and_give_new_cells(Vertex_handle v, const Point& p, + OutputItCells fit) { Self tmp; - Vertex_remover remover (tmp); - Vertex_inserter inserter (*this); + Vertex_remover remover(tmp); + Vertex_inserter inserter(*this); Vertex_handle res = Tr_Base::move_if_no_collision_and_give_new_cells(v,p, remover,inserter,fit); @@ -1435,96 +1438,95 @@ move_if_no_collision_and_give_new_cells(Vertex_handle v, const Point &p, template < class Gt, class Tds, class Lds > Oriented_side Delaunay_triangulation_3:: -side_of_oriented_sphere(const Point &p0, const Point &p1, const Point &p2, - const Point &p3, const Point &p, bool perturb) const +side_of_oriented_sphere(const Point& p0, const Point& p1, const Point& p2, + const Point& p3, const Point& p, bool perturb) const { - CGAL_triangulation_precondition( orientation(p0, p1, p2, p3) == POSITIVE ); + CGAL_triangulation_precondition(orientation(p0, p1, p2, p3) == POSITIVE); - Oriented_side os = - geom_traits().side_of_oriented_sphere_3_object()(p0, p1, p2, p3, p); + Oriented_side os = + geom_traits().side_of_oriented_sphere_3_object()(p0, p1, p2, p3, p); - if (os != ON_ORIENTED_BOUNDARY || !perturb) - return os; + if(os != ON_ORIENTED_BOUNDARY || !perturb) + return os; - // We are now in a degenerate case => we do a symbolic perturbation. + // We are now in a degenerate case => we do a symbolic perturbation. - // We sort the points lexicographically. - const Point * points[5] = {&p0, &p1, &p2, &p3, &p}; - std::sort(points, points + 5, typename Tr_Base::Perturbation_order(this)); + // We sort the points lexicographically. + const Point * points[5] = {&p0, &p1, &p2, &p3, &p}; + std::sort(points, points + 5, typename Tr_Base::Perturbation_order(this)); - // We successively look whether the leading monomial, then 2nd monomial - // of the determinant has non null coefficient. - // 2 iterations are enough (cf paper) - for (int i=4; i>2; --i) { - if (points[i] == &p) - return ON_NEGATIVE_SIDE; // since p0 p1 p2 p3 are non coplanar - // and positively oriented - Orientation o; - if (points[i] == &p3 && (o = orientation(p0,p1,p2,p)) != COPLANAR ) - return o; - if (points[i] == &p2 && (o = orientation(p0,p1,p,p3)) != COPLANAR ) - return o; - if (points[i] == &p1 && (o = orientation(p0,p,p2,p3)) != COPLANAR ) - return o; - if (points[i] == &p0 && (o = orientation(p,p1,p2,p3)) != COPLANAR ) - return o; - } + // We successively look whether the leading monomial, then 2nd monomial + // of the determinant has non null coefficient. + // 2 iterations are enough (cf paper) + for(int i=4; i>2; --i) + { + if(points[i] == &p) + return ON_NEGATIVE_SIDE; // since p0 p1 p2 p3 are non coplanar + // and positively oriented + Orientation o; + if(points[i] == &p3 && (o = orientation(p0,p1,p2,p)) != COPLANAR) + return o; + if(points[i] == &p2 && (o = orientation(p0,p1,p,p3)) != COPLANAR) + return o; + if(points[i] == &p1 && (o = orientation(p0,p,p2,p3)) != COPLANAR) + return o; + if(points[i] == &p0 && (o = orientation(p,p1,p2,p3)) != COPLANAR) + return o; + } - CGAL_triangulation_assertion(false); - return ON_NEGATIVE_SIDE; + CGAL_triangulation_assertion(false); + return ON_NEGATIVE_SIDE; } template < class Gt, class Tds, class Lds > Bounded_side Delaunay_triangulation_3:: -coplanar_side_of_bounded_circle(const Point &p0, const Point &p1, - const Point &p2, const Point &p, bool perturb) const +coplanar_side_of_bounded_circle(const Point& p0, const Point& p1, + const Point& p2, const Point& p, bool perturb) const { - // In dim==2, we should even be able to assert orient == POSITIVE. - CGAL_triangulation_precondition( coplanar_orientation(p0, p1, p2) - != COLLINEAR ); + // In dim==2, we should even be able to assert orient == POSITIVE. + CGAL_triangulation_precondition(coplanar_orientation(p0, p1, p2) + != COLLINEAR); - Bounded_side bs = + Bounded_side bs = geom_traits().coplanar_side_of_bounded_circle_3_object()(p0, p1, p2, p); - if (bs != ON_BOUNDARY || !perturb) - return bs; + if(bs != ON_BOUNDARY || !perturb) + return bs; - // We are now in a degenerate case => we do a symbolic perturbation. + // We are now in a degenerate case => we do a symbolic perturbation. - // We sort the points lexicographically. - const Point * points[4] = {&p0, &p1, &p2, &p}; - std::sort(points, points + 4, typename Tr_Base::Perturbation_order(this)); + // We sort the points lexicographically. + const Point * points[4] = {&p0, &p1, &p2, &p}; + std::sort(points, points + 4, typename Tr_Base::Perturbation_order(this)); - Orientation local = coplanar_orientation(p0, p1, p2); + Orientation local = coplanar_orientation(p0, p1, p2); - // we successively look whether the leading monomial, then 2nd monimial, - // then 3rd monomial, of the determinant which has non null coefficient - // [syl] : TODO : Probably it can be stopped earlier like the 3D version - for (int i=3; i>0; --i) { - if (points[i] == &p) - return Bounded_side(NEGATIVE); // since p0 p1 p2 are non collinear - // but not necessarily positively oriented - Orientation o; - if (points[i] == &p2 - && (o = coplanar_orientation(p0,p1,p)) != COLLINEAR ) - // [syl] : TODO : I'm not sure of the signs here (nor the rest :) - return Bounded_side(o*local); - if (points[i] == &p1 - && (o = coplanar_orientation(p0,p,p2)) != COLLINEAR ) - return Bounded_side(o*local); - if (points[i] == &p0 - && (o = coplanar_orientation(p,p1,p2)) != COLLINEAR ) - return Bounded_side(o*local); - } + // we successively look whether the leading monomial, then 2nd monimial, + // then 3rd monomial, of the determinant which has non null coefficient + // [syl] : TODO : Probably it can be stopped earlier like the 3D version + for(int i=3; i>0; --i) + { + if(points[i] == &p) + return Bounded_side(NEGATIVE); // since p0 p1 p2 are non collinear + // but not necessarily positively oriented + Orientation o; + if(points[i] == &p2 && (o = coplanar_orientation(p0,p1,p)) != COLLINEAR) + // [syl] : TODO : I'm not sure of the signs here (nor the rest :) + return Bounded_side(o*local); + if(points[i] == &p1 && (o = coplanar_orientation(p0,p,p2)) != COLLINEAR) + return Bounded_side(o*local); + if(points[i] == &p0 && (o = coplanar_orientation(p,p1,p2)) != COLLINEAR) + return Bounded_side(o*local); + } - // case when the first non null coefficient is the coefficient of - // the 4th monomial - // moreover, the tests (points[] == &p) were false up to here, so the - // monomial corresponding to p is the only monomial with non-zero - // coefficient, it is equal to coplanar_orient(p0,p1,p2) == positive - // so, no further test is required - return Bounded_side(-local); //ON_UNBOUNDED_SIDE; + // case when the first non null coefficient is the coefficient of + // the 4th monomial + // moreover, the tests (points[] == &p) were false up to here, so the + // monomial corresponding to p is the only monomial with non-zero + // coefficient, it is equal to coplanar_orient(p0,p1,p2) == positive + // so, no further test is required + return Bounded_side(-local); //ON_UNBOUNDED_SIDE; } template < class Gt, class Tds, class Lds > @@ -1532,46 +1534,54 @@ Bounded_side Delaunay_triangulation_3:: side_of_sphere(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, - const Point &p, bool perturb) const + const Point& p, bool perturb) const { - CGAL_triangulation_precondition( dimension() == 3 ); + CGAL_triangulation_precondition(dimension() == 3); - if (is_infinite(v0)) { - Orientation o = orientation(v2->point(), v1->point(), v3->point(), p); - if (o != COPLANAR) - return Bounded_side(o); - return coplanar_side_of_bounded_circle(v2->point(), v1->point(), v3->point(), p, perturb); - } + if(is_infinite(v0)) + { + Orientation o = orientation(v2->point(), v1->point(), v3->point(), p); + if(o != COPLANAR) + return Bounded_side(o); - if (is_infinite(v1)) { - Orientation o = orientation(v2->point(), v3->point(), v0->point(), p); - if (o != COPLANAR) - return Bounded_side(o); - return coplanar_side_of_bounded_circle(v2->point(), v3->point(), v0->point(), p, perturb); - } + return coplanar_side_of_bounded_circle(v2->point(), v1->point(), v3->point(), p, perturb); + } - if (is_infinite(v2)) { - Orientation o = orientation(v1->point(), v0->point(), v3->point(), p); - if (o != COPLANAR) - return Bounded_side(o); - return coplanar_side_of_bounded_circle(v1->point(), v0->point(), v3->point(), p, perturb); - } + if(is_infinite(v1)) + { + Orientation o = orientation(v2->point(), v3->point(), v0->point(), p); + if(o != COPLANAR) + return Bounded_side(o); - if (is_infinite(v3)) { - Orientation o = orientation(v0->point(), v1->point(), v2->point(), p); - if (o != COPLANAR) - return Bounded_side(o); - return coplanar_side_of_bounded_circle(v0->point(), v1->point(), v2->point(), p, perturb); - } + return coplanar_side_of_bounded_circle(v2->point(), v3->point(), v0->point(), p, perturb); + } - return (Bounded_side) side_of_oriented_sphere(v0->point(), v1->point(), v2->point(), v3->point(), p, perturb); + if(is_infinite(v2)) + { + Orientation o = orientation(v1->point(), v0->point(), v3->point(), p); + if(o != COPLANAR) + return Bounded_side(o); + + return coplanar_side_of_bounded_circle(v1->point(), v0->point(), v3->point(), p, perturb); + } + + if(is_infinite(v3)) + { + Orientation o = orientation(v0->point(), v1->point(), v2->point(), p); + if(o != COPLANAR) + return Bounded_side(o); + + return coplanar_side_of_bounded_circle(v0->point(), v1->point(), v2->point(), p, perturb); + } + + return (Bounded_side) side_of_oriented_sphere(v0->point(), v1->point(), v2->point(), v3->point(), p, perturb); } template < class Gt, class Tds, class Lds > Bounded_side Delaunay_triangulation_3:: -side_of_circle(Cell_handle c, int i, - const Point & p, bool perturb) const +side_of_circle(Cell_handle c, int i, const Point& p, bool perturb) const +{ // precondition : dimension >=2 // in dimension 3, - for a finite facet // returns ON_BOUNDARY if the point lies on the circle, @@ -1585,53 +1595,53 @@ side_of_circle(Cell_handle c, int i, // finite edge (endpoints included) // ON_BOUNDED_SIDE for a point in the open half-plane // ON_UNBOUNDED_SIDE elsewhere -{ - CGAL_triangulation_precondition( dimension() >= 2 ); + + CGAL_triangulation_precondition(dimension() >= 2); int i3 = 5; - if ( dimension() == 2 ) { - CGAL_triangulation_precondition( i == 3 ); + if(dimension() == 2) + { + CGAL_triangulation_precondition(i == 3); // the triangulation is supposed to be valid, ie the facet // with vertices 0 1 2 in this order is positively oriented - if ( ! c->has_vertex( infinite_vertex(), i3 ) ) - return coplanar_side_of_bounded_circle( c->vertex(0)->point(), + if(! c->has_vertex(infinite_vertex(), i3)) + return coplanar_side_of_bounded_circle(c->vertex(0)->point(), c->vertex(1)->point(), c->vertex(2)->point(), p, perturb); // else infinite facet // v1, v2 finite vertices of the facet such that v1,v2,infinite // is positively oriented - Vertex_handle v1 = c->vertex( ccw(i3) ), - v2 = c->vertex( cw(i3) ); + Vertex_handle v1 = c->vertex(ccw(i3)), + v2 = c->vertex(cw(i3)); CGAL_triangulation_assertion(coplanar_orientation(v1->point(), v2->point(), mirror_vertex(c, i3)->point()) == NEGATIVE); Orientation o = coplanar_orientation(v1->point(), v2->point(), p); - if ( o != COLLINEAR ) - return Bounded_side( o ); + if(o != COLLINEAR) + return Bounded_side(o); // because p is in f iff // it does not lie on the same side of v1v2 as vn int i_e; Locate_type lt; // case when p collinear with v1v2 - return side_of_segment( p, - v1->point(), v2->point(), - lt, i_e ); + return side_of_segment(p, v1->point(), v2->point(), lt, i_e); } // else dimension == 3 - CGAL_triangulation_precondition( i >= 0 && i < 4 ); - if ( ( ! c->has_vertex(infinite_vertex(),i3) ) || ( i3 != i ) ) { + CGAL_triangulation_precondition(i >= 0 && i < 4); + if((! c->has_vertex(infinite_vertex(),i3)) || (i3 != i)) + { // finite facet // initialization of i0 i1 i2, vertices of the facet positively // oriented (if the triangulation is valid) int i0 = (i>0) ? 0 : 1; int i1 = (i>1) ? 1 : 2; int i2 = (i>2) ? 2 : 3; - CGAL_triangulation_precondition( coplanar( c->vertex(i0)->point(), + CGAL_triangulation_precondition(coplanar(c->vertex(i0)->point(), c->vertex(i1)->point(), c->vertex(i2)->point(), - p ) ); - return coplanar_side_of_bounded_circle( c->vertex(i0)->point(), + p)); + return coplanar_side_of_bounded_circle(c->vertex(i0)->point(), c->vertex(i1)->point(), c->vertex(i2)->point(), p, perturb); @@ -1640,40 +1650,39 @@ side_of_circle(Cell_handle c, int i, //else infinite facet // v1, v2 finite vertices of the facet such that v1,v2,infinite // is positively oriented - Vertex_handle v1 = c->vertex( next_around_edge(i3,i) ), - v2 = c->vertex( next_around_edge(i,i3) ); + Vertex_handle v1 = c->vertex(next_around_edge(i3,i)), + v2 = c->vertex(next_around_edge(i,i3)); Orientation o = (Orientation) - (coplanar_orientation( v1->point(), v2->point(), + (coplanar_orientation(v1->point(), v2->point(), c->vertex(i)->point()) * - coplanar_orientation( v1->point(), v2->point(), p )); + coplanar_orientation(v1->point(), v2->point(), p)); // then the code is duplicated from 2d case - if ( o != COLLINEAR ) - return Bounded_side( -o ); + if(o != COLLINEAR) + return Bounded_side(-o); // because p is in f iff // it is not on the same side of v1v2 as c->vertex(i) int i_e; Locate_type lt; // case when p collinear with v1v2 - return side_of_segment( p, - v1->point(), v2->point(), - lt, i_e ); + return side_of_segment(p, v1->point(), v2->point(), lt, i_e); } template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: nearest_vertex_in_cell(const Point& p, Cell_handle c) const -// Returns the finite vertex of the cell c which is the closest to p. { - CGAL_triangulation_precondition(dimension() >= 0); +// Returns the finite vertex of the cell c which is the closest to p. + CGAL_triangulation_precondition(dimension() >= 0); - Vertex_handle nearest = nearest_vertex(p, c->vertex(0), c->vertex(1)); - if (dimension() >= 2) { - nearest = nearest_vertex(p, nearest, c->vertex(2)); - if (dimension() == 3) - nearest = nearest_vertex(p, nearest, c->vertex(3)); - } - return nearest; + Vertex_handle nearest = nearest_vertex(p, c->vertex(0), c->vertex(1)); + if(dimension() >= 2) + { + nearest = nearest_vertex(p, nearest, c->vertex(2)); + if(dimension() == 3) + nearest = nearest_vertex(p, nearest, c->vertex(3)); + } + return nearest; } template < class Gt, class Tds, class Lds > @@ -1681,44 +1690,50 @@ typename Delaunay_triangulation_3::Vertex_handle Delaunay_triangulation_3:: nearest_vertex(const Point& p, Cell_handle start) const { - if (number_of_vertices() == 0) - return Vertex_handle(); + if(number_of_vertices() == 0) + return Vertex_handle(); - // Use a brute-force algorithm if dimension < 3. - if (dimension() < 3) { - Finite_vertices_iterator vit = finite_vertices_begin(); - Vertex_handle res = vit; - ++vit; - for (Finite_vertices_iterator end = finite_vertices_end(); vit != end; ++vit) - res = nearest_vertex(p, res, vit); - return res; + // Use a brute-force algorithm if dimension < 3. + if(dimension() < 3) + { + Finite_vertices_iterator vit = finite_vertices_begin(); + Vertex_handle res = vit; + ++vit; + for(Finite_vertices_iterator end = finite_vertices_end(); vit != end; ++vit) + res = nearest_vertex(p, res, vit); + return res; + } + + Locate_type lt; + int li, lj; + Cell_handle c = locate(p, lt, li, lj, start); + if(lt == Tr_Base::VERTEX) + return c->vertex(li); + + // - start with the closest vertex from the located cell. + // - repeatedly take the nearest of its incident vertices if any + // - if not, we're done. + Vertex_handle nearest = nearest_vertex_in_cell(p, c); + std::vector vs; + vs.reserve(32); + while(true) + { + Vertex_handle tmp = nearest; + adjacent_vertices(nearest, std::back_inserter(vs)); + for(typename std::vector::const_iterator + vsit = vs.begin(); vsit != vs.end(); ++vsit) + { + tmp = nearest_vertex(p, tmp, *vsit); } - Locate_type lt; - int li, lj; - Cell_handle c = locate(p, lt, li, lj, start); - if (lt == Tr_Base::VERTEX) - return c->vertex(li); + if(tmp == nearest) + break; - // - start with the closest vertex from the located cell. - // - repeatedly take the nearest of its incident vertices if any - // - if not, we're done. - Vertex_handle nearest = nearest_vertex_in_cell(p, c); - std::vector vs; - vs.reserve(32); - while (true) { - Vertex_handle tmp = nearest; - adjacent_vertices(nearest, std::back_inserter(vs)); - for (typename std::vector::const_iterator - vsit = vs.begin(); vsit != vs.end(); ++vsit) - tmp = nearest_vertex(p, tmp, *vsit); - if (tmp == nearest) - break; - vs.clear(); - nearest = tmp; - } + vs.clear(); + nearest = tmp; + } - return nearest; + return nearest; } // This is not a fast version. @@ -1731,12 +1746,14 @@ nearest_vertex(const Point& p, Cell_handle start) const template bool Delaunay_triangulation_3:: -is_delaunay_after_displacement(Vertex_handle v, const Point &p) const +is_delaunay_after_displacement(Vertex_handle v, const Point& p) const { CGAL_triangulation_precondition(!this->is_infinite(v)); CGAL_triangulation_precondition(this->dimension() == 2); CGAL_triangulation_precondition(!this->test_dim_down(v)); - if(v->point() == p) return true; + if(v->point() == p) + return true; + Point ant = v->point(); v->set_point(p); @@ -1752,8 +1769,7 @@ is_delaunay_after_displacement(Vertex_handle v, const Point &p) const Cell_handle c = cells[i]; if(this->is_infinite(c)) continue; if(this->orientation(c->vertex(0)->point(), c->vertex(1)->point(), - c->vertex(2)->point(), c->vertex(3)->point()) - != POSITIVE) + c->vertex(2)->point(), c->vertex(3)->point()) != POSITIVE) { v->set_point(ant); return false; @@ -1767,20 +1783,24 @@ is_delaunay_after_displacement(Vertex_handle v, const Point &p) const size = facets.size(); for(std::size_t i=0; ineighbor(j); int mj = this->mirror_index(c, j); Vertex_handle h1 = c->vertex(j); - if(this->is_infinite(h1)) { - if(this->side_of_sphere(c, cj->vertex(mj)->point(), true) - != ON_UNBOUNDED_SIDE) { + if(this->is_infinite(h1)) + { + if(this->side_of_sphere(c, cj->vertex(mj)->point(), true) != ON_UNBOUNDED_SIDE) + { v->set_point(ant); return false; } - } else { - if(this->side_of_sphere(cj, h1->point(), true) != ON_UNBOUNDED_SIDE) { + } + else + { + if(this->side_of_sphere(cj, h1->point(), true) != ON_UNBOUNDED_SIDE) + { v->set_point(ant); return false; } @@ -1805,25 +1825,25 @@ Delaunay_triangulation_3:: is_Gabriel(Cell_handle c, int i) const { CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i)); - typename Geom_traits::Side_of_bounded_sphere_3 - side_of_bounded_sphere = + typename Geom_traits::Side_of_bounded_sphere_3 side_of_bounded_sphere = geom_traits().side_of_bounded_sphere_3_object(); - if ((!is_infinite(c->vertex(i))) && - side_of_bounded_sphere ( - c->vertex(vertex_triple_index(i,0))->point(), - c->vertex(vertex_triple_index(i,1))->point(), - c->vertex(vertex_triple_index(i,2))->point(), - c->vertex(i)->point()) == ON_BOUNDED_SIDE ) return false; - Cell_handle neighbor = c->neighbor(i); + if((!is_infinite(c->vertex(i))) && + side_of_bounded_sphere (c->vertex(vertex_triple_index(i,0))->point(), + c->vertex(vertex_triple_index(i,1))->point(), + c->vertex(vertex_triple_index(i,2))->point(), + c->vertex(i)->point()) == ON_BOUNDED_SIDE) + return false; + + Cell_handle neighbor = c->neighbor(i); int in = neighbor->index(c); - if ((!is_infinite(neighbor->vertex(in))) && - side_of_bounded_sphere( - c->vertex(vertex_triple_index(i,0))->point(), - c->vertex(vertex_triple_index(i,1))->point(), - c->vertex(vertex_triple_index(i,2))->point(), - neighbor->vertex(in)->point()) == ON_BOUNDED_SIDE ) return false; + if((!is_infinite(neighbor->vertex(in))) && + side_of_bounded_sphere(c->vertex(vertex_triple_index(i,0))->point(), + c->vertex(vertex_triple_index(i,1))->point(), + c->vertex(vertex_triple_index(i,2))->point(), + neighbor->vertex(in)->point()) == ON_BOUNDED_SIDE) + return false; return true; } @@ -1842,25 +1862,24 @@ Delaunay_triangulation_3:: is_Gabriel(Cell_handle c, int i, int j) const { CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i,j)); - typename Geom_traits::Side_of_bounded_sphere_3 - side_of_bounded_sphere = + typename Geom_traits::Side_of_bounded_sphere_3 side_of_bounded_sphere = geom_traits().side_of_bounded_sphere_3_object(); - Facet_circulator fcirc = incident_facets(c,i,j), - fdone(fcirc); + Facet_circulator fcirc = incident_facets(c,i,j), fdone(fcirc); Vertex_handle v1 = c->vertex(i); Vertex_handle v2 = c->vertex(j); - do { + do + { // test whether the vertex of cc opposite to *fcirc // is inside the sphere defined by the edge e = (s, i,j) Cell_handle cc = (*fcirc).first; int ii = (*fcirc).second; - if (!is_infinite(cc->vertex(ii)) && - side_of_bounded_sphere( v1->point(), - v2->point(), - cc->vertex(ii)->point()) - == ON_BOUNDED_SIDE ) return false; - } while(++fcirc != fdone); + if(!is_infinite(cc->vertex(ii)) && + side_of_bounded_sphere(v1->point(), v2->point(), cc->vertex(ii)->point()) + == ON_BOUNDED_SIDE) return false; + } + while(++fcirc != fdone); + return true; } @@ -1870,43 +1889,48 @@ Delaunay_triangulation_3:: dual(Cell_handle c) const { CGAL_triangulation_precondition(dimension()==3); - CGAL_triangulation_precondition( ! is_infinite(c) ); + CGAL_triangulation_precondition(! is_infinite(c)); return c->circumcenter(geom_traits()); } - template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Object Delaunay_triangulation_3:: dual(Cell_handle c, int i) const { CGAL_triangulation_precondition(dimension()>=2); - CGAL_triangulation_precondition( ! is_infinite(c,i) ); + CGAL_triangulation_precondition(! is_infinite(c,i)); - if ( dimension() == 2 ) { - CGAL_triangulation_precondition( i == 3 ); - return construct_object( construct_circumcenter(c->vertex(0)->point(), - c->vertex(1)->point(), - c->vertex(2)->point()) ); + if(dimension() == 2) + { + CGAL_triangulation_precondition(i == 3); + return construct_object(construct_circumcenter(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point())); } // dimension() == 3 Cell_handle n = c->neighbor(i); - if ( ! is_infinite(c) && ! is_infinite(n) ) - return construct_object(construct_segment( dual(c), dual(n) )); + if(! is_infinite(c) && ! is_infinite(n)) + return construct_object(construct_segment(dual(c), dual(n))); // either n or c is infinite int in; - if ( is_infinite(c) ) + if(is_infinite(c)) + { in = n->index(c); - else { + } + else + { n = c; in = i; } + // n now denotes a finite cell, either c or c->neighbor(i) int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; - if ( (in&1) == 1 ) + if((in&1) == 1) std::swap(ind[0], ind[1]); + // in=0: 1 2 3 // in=1: 3 2 0 // in=2: 3 0 1 @@ -1915,66 +1939,71 @@ dual(Cell_handle c, int i) const const Point& q = n->vertex(ind[1])->point(); const Point& r = n->vertex(ind[2])->point(); - Line l = construct_equidistant_line( p, q, r ); - return construct_object(construct_ray( dual(n), l)); + Line l = construct_equidistant_line(p, q, r); + return construct_object(construct_ray(dual(n), l)); } - - template < class Gt, class Tds, class Lds > typename Delaunay_triangulation_3::Line Delaunay_triangulation_3:: dual_support(Cell_handle c, int i) const { CGAL_triangulation_precondition(dimension()>=2); - CGAL_triangulation_precondition( ! is_infinite(c,i) ); + CGAL_triangulation_precondition(! is_infinite(c,i)); - if ( dimension() == 2 ) { - CGAL_triangulation_precondition( i == 3 ); - return construct_equidistant_line( c->vertex(0)->point(), + if(dimension() == 2) + { + CGAL_triangulation_precondition(i == 3); + return construct_equidistant_line(c->vertex(0)->point(), c->vertex(1)->point(), - c->vertex(2)->point() ); + c->vertex(2)->point()); } - return construct_equidistant_line( c->vertex((i+1)&3)->point(), + return construct_equidistant_line(c->vertex((i+1)&3)->point(), c->vertex((i+2)&3)->point(), - c->vertex((i+3)&3)->point() ); + c->vertex((i+3)&3)->point()); } - template < class Gt, class Tds, class Lds > bool Delaunay_triangulation_3:: is_valid(bool verbose, int level) const { - if ( ! tds().is_valid(verbose,level) ) { - if (verbose) - std::cerr << "invalid data structure" << std::endl; + if(! tds().is_valid(verbose,level)) + { + if(verbose) + std::cerr << "invalid data structure" << std::endl; + CGAL_triangulation_assertion(false); return false; } - if ( infinite_vertex() == Vertex_handle() ) { - if (verbose) - std::cerr << "no infinite vertex" << std::endl; + if(infinite_vertex() == Vertex_handle()) + { + if(verbose) + std::cerr << "no infinite vertex" << std::endl; + CGAL_triangulation_assertion(false); return false; } - switch ( dimension() ) { - case 3: + switch(dimension()) + { + case 3: { - for(Finite_cells_iterator it = finite_cells_begin(), end = finite_cells_end(); it != end; ++it) { + for(Finite_cells_iterator it = finite_cells_begin(), end = finite_cells_end(); it != end; ++it) + { is_valid_finite(it); - for(int i=0; i<4; i++ ) { - if ( !is_infinite - (it->neighbor(i)->vertex(it->neighbor(i)->index(it))) ) { - if ( side_of_sphere - (it, - it->neighbor(i)->vertex(it->neighbor(i)->index(it))->point()) - == ON_BOUNDED_SIDE ) { - if (verbose) + for(int i=0; i<4; i++) + { + if(!is_infinite(it->neighbor(i)->vertex(it->neighbor(i)->index(it)))) + { + if(side_of_sphere(it, it->neighbor(i)->vertex( + it->neighbor(i)->index(it))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) std::cerr << "non-empty sphere " << std::endl; + CGAL_triangulation_assertion(false); return false; } @@ -1983,21 +2012,24 @@ is_valid(bool verbose, int level) const } break; } - case 2: + case 2: { - for(Finite_facets_iterator it = finite_facets_begin(), end = finite_facets_end(); it != end; ++it) { + for(Finite_facets_iterator it = finite_facets_begin(), end = finite_facets_end(); it != end; ++it) + { is_valid_finite((*it).first); - for(int i=0; i<3; i++ ) { - if( !is_infinite - ((*it).first->neighbor(i)->vertex( (((*it).first)->neighbor(i)) - ->index((*it).first))) ) { - if ( side_of_circle ( (*it).first, 3, - (*it).first->neighbor(i)-> - vertex( (((*it).first)->neighbor(i)) - ->index((*it).first) )->point() ) - == ON_BOUNDED_SIDE ) { - if (verbose) + for(int i=0; i<3; i++) + { + if(!is_infinite((*it).first->neighbor(i)->vertex( + (((*it).first)->neighbor(i))->index( + (*it).first)))) + { + if(side_of_circle((*it).first, 3, (*it).first->neighbor(i)-> + vertex((((*it).first)->neighbor(i))->index( + (*it).first))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) std::cerr << "non-empty circle " << std::endl; + CGAL_triangulation_assertion(false); return false; } @@ -2006,15 +2038,18 @@ is_valid(bool verbose, int level) const } break; } - case 1: + case 1: { for(Finite_edges_iterator it = finite_edges_begin(), end = finite_edges_end(); it != end; ++it) + { is_valid_finite((*it).first); + } break; } } - if (verbose) + if(verbose) std::cerr << "Delaunay valid triangulation" << std::endl; + return true; } @@ -2023,26 +2058,31 @@ bool Delaunay_triangulation_3:: is_valid(Cell_handle c, bool verbose, int level) const { - if ( ! Tr_Base::is_valid(c,verbose,level) ) { - if (verbose) { + if(! Tr_Base::is_valid(c, verbose, level)) + { + if(verbose) + { std::cerr << "combinatorically invalid cell" ; - for (int i=0; i <= dimension(); i++ ) + for(int i=0; i <= dimension(); i++) std::cerr << c->vertex(i)->point() << ", " ; std::cerr << std::endl; } CGAL_triangulation_assertion(false); return false; } - switch ( dimension() ) { - case 3: + switch(dimension()) + { + case 3: { - if ( ! is_infinite(c) ) { - is_valid_finite(c,verbose,level); - for (int i=0; i<4; i++ ) { - if (side_of_sphere(c, c->vertex((c->neighbor(i))->index(c))->point()) - == ON_BOUNDED_SIDE ) { - if (verbose) - std::cerr << "non-empty sphere " << std::endl; + if(! is_infinite(c)) + { + is_valid_finite(c, verbose, level); + for(int i=0; i<4; i++) { + if(side_of_sphere(c, c->vertex((c->neighbor(i))->index(c))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) + std::cerr << "non-empty sphere " << std::endl; + CGAL_triangulation_assertion(false); return false; } @@ -2050,14 +2090,17 @@ is_valid(Cell_handle c, bool verbose, int level) const } break; } - case 2: + case 2: { - if ( ! is_infinite(c,3) ) { - for (int i=0; i<2; i++ ) { - if (side_of_circle(c, 3, c->vertex(c->neighbor(i)->index(c))->point()) - == ON_BOUNDED_SIDE ) { - if (verbose) - std::cerr << "non-empty circle " << std::endl; + if(! is_infinite(c,3)) + { + for(int i=0; i<2; i++) + { + if(side_of_circle(c, 3, c->vertex(c->neighbor(i)->index(c))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) + std::cerr << "non-empty circle " << std::endl; + CGAL_triangulation_assertion(false); return false; } @@ -2066,8 +2109,9 @@ is_valid(Cell_handle c, bool verbose, int level) const break; } } - if (verbose) - std::cerr << "Delaunay valid cell" << std::endl; + if(verbose) + std::cerr << "Delaunay valid cell" << std::endl; + return true; } diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index 28bfcb5108d..c6fb7ecbab3 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -74,1237 +74,1264 @@ namespace CGAL { - /************************************************ +/************************************************ * * Regular_triangulation_3 class * ************************************************/ - template < class Gt, class Tds_ = Default, class Lock_data_structure_ = Default > - class Regular_triangulation_3 +template < class Gt, class Tds_ = Default, class Lock_data_structure_ = Default > +class Regular_triangulation_3 : public Triangulation_3< - Gt, - typename Default::Get, - Regular_triangulation_cell_base_3 > >::type, - Lock_data_structure_> - { - private: - typedef typename Default::Get, - Regular_triangulation_cell_base_3 > >::type Tds; + Gt, + typename Default::Get, + Regular_triangulation_cell_base_3 > >::type, + Lock_data_structure_> +{ +private: + typedef typename Default::Get, + Regular_triangulation_cell_base_3 > + >::type Tds; - typedef Regular_triangulation_3 Self; + typedef Regular_triangulation_3 Self; - public: - typedef Triangulation_3 Tr_Base; +public: + typedef Triangulation_3 Tr_Base; - typedef Gt Geom_traits; - typedef Tds Triangulation_data_structure; + typedef Gt Geom_traits; + typedef Tds Triangulation_data_structure; - typedef Geom_traits Traits; - typedef typename Tr_Base::Concurrency_tag Concurrency_tag; - typedef typename Tr_Base::Lock_data_structure Lock_data_structure; + typedef Geom_traits Traits; + typedef typename Tr_Base::Concurrency_tag Concurrency_tag; + typedef typename Tr_Base::Lock_data_structure Lock_data_structure; - typedef typename Tr_Base::Vertex_handle Vertex_handle; - typedef typename Tr_Base::Cell_handle Cell_handle; - typedef typename Tr_Base::Vertex Vertex; - typedef typename Tr_Base::Cell Cell; - typedef typename Tr_Base::Facet Facet; - typedef typename Tr_Base::Edge Edge; + typedef typename Tr_Base::Vertex_handle Vertex_handle; + typedef typename Tr_Base::Cell_handle Cell_handle; + typedef typename Tr_Base::Vertex Vertex; + typedef typename Tr_Base::Cell Cell; + typedef typename Tr_Base::Facet Facet; + typedef typename Tr_Base::Edge Edge; - typedef typename Tr_Base::size_type size_type; - typedef typename Tr_Base::Locate_type Locate_type; - typedef typename Tr_Base::Cell_iterator Cell_iterator; - typedef typename Tr_Base::Facet_iterator Facet_iterator; - typedef typename Tr_Base::Edge_iterator Edge_iterator; - typedef typename Tr_Base::Facet_circulator Facet_circulator; + typedef typename Tr_Base::size_type size_type; + typedef typename Tr_Base::Locate_type Locate_type; + typedef typename Tr_Base::Cell_iterator Cell_iterator; + typedef typename Tr_Base::Facet_iterator Facet_iterator; + typedef typename Tr_Base::Edge_iterator Edge_iterator; + typedef typename Tr_Base::Facet_circulator Facet_circulator; - typedef typename Tr_Base::Finite_vertices_iterator Finite_vertices_iterator; - typedef typename Tr_Base::Finite_cells_iterator Finite_cells_iterator; - typedef typename Tr_Base::Finite_facets_iterator Finite_facets_iterator; - typedef typename Tr_Base::Finite_edges_iterator Finite_edges_iterator; - typedef typename Tr_Base::All_cells_iterator All_cells_iterator; + typedef typename Tr_Base::Finite_vertices_iterator Finite_vertices_iterator; + typedef typename Tr_Base::Finite_cells_iterator Finite_cells_iterator; + typedef typename Tr_Base::Finite_facets_iterator Finite_facets_iterator; + typedef typename Tr_Base::Finite_edges_iterator Finite_edges_iterator; + typedef typename Tr_Base::All_cells_iterator All_cells_iterator; - // Traits are not supposed to define Bare_point, but leaving below - // for backward compatibility - typedef typename boost::mpl::eval_if_c< - internal::Has_nested_type_Bare_point::value, - typename internal::Bare_point_type, - boost::mpl::identity - >::type Bare_point; - typedef typename Gt::Weighted_point_3 Weighted_point; + // Traits are not supposed to define Bare_point, but leaving below + // for backward compatibility + typedef typename boost::mpl::eval_if_c< + internal::Has_nested_type_Bare_point::value, + typename internal::Bare_point_type, + boost::mpl::identity + >::type Bare_point; + typedef typename Gt::Weighted_point_3 Weighted_point; - typedef typename Gt::Segment_3 Segment; - typedef typename Gt::Triangle_3 Triangle; - typedef typename Gt::Tetrahedron_3 Tetrahedron; + typedef typename Gt::Segment_3 Segment; + typedef typename Gt::Triangle_3 Triangle; + typedef typename Gt::Tetrahedron_3 Tetrahedron; - // types for dual: - typedef typename Gt::Line_3 Line; - typedef typename Gt::Ray_3 Ray; - typedef typename Gt::Plane_3 Plane; - typedef typename Gt::Object_3 Object; + // types for dual: + typedef typename Gt::Line_3 Line; + typedef typename Gt::Ray_3 Ray; + typedef typename Gt::Plane_3 Plane; + typedef typename Gt::Object_3 Object; - //Tag to distinguish Delaunay from regular triangulations - typedef Tag_true Weighted_tag; + //Tag to distinguish Delaunay from regular triangulations + typedef Tag_true Weighted_tag; - // Tag to distinguish periodic triangulations from others - typedef Tag_false Periodic_tag; + // Tag to distinguish periodic triangulations from others + typedef Tag_false Periodic_tag; #ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 - using Tr_Base::geom_traits; + using Tr_Base::geom_traits; #endif - using Tr_Base::adjacent_vertices; - using Tr_Base::cw; - using Tr_Base::ccw; - using Tr_Base::construct_point; - using Tr_Base::coplanar_orientation; - using Tr_Base::dimension; - using Tr_Base::find_conflicts; - using Tr_Base::finite_facets_begin; - using Tr_Base::finite_facets_end; - using Tr_Base::finite_vertices_begin; - using Tr_Base::finite_vertices_end; - using Tr_Base::finite_cells_begin; - using Tr_Base::finite_cells_end; - using Tr_Base::finite_edges_begin; - using Tr_Base::finite_edges_end; - using Tr_Base::incident_facets; - using Tr_Base::insert_in_conflict; - using Tr_Base::infinite_vertex; - using Tr_Base::is_infinite; - using Tr_Base::is_valid; - using Tr_Base::is_valid_finite; - using Tr_Base::locate; - using Tr_Base::mirror_vertex; - using Tr_Base::mirror_index; - using Tr_Base::next_around_edge; - using Tr_Base::number_of_vertices; - using Tr_Base::orientation; - using Tr_Base::point; - using Tr_Base::side_of_segment; - using Tr_Base::side_of_edge; - using Tr_Base::tds; - using Tr_Base::vertex_triple_index; + using Tr_Base::adjacent_vertices; + using Tr_Base::cw; + using Tr_Base::ccw; + using Tr_Base::construct_point; + using Tr_Base::coplanar_orientation; + using Tr_Base::dimension; + using Tr_Base::find_conflicts; + using Tr_Base::finite_facets_begin; + using Tr_Base::finite_facets_end; + using Tr_Base::finite_vertices_begin; + using Tr_Base::finite_vertices_end; + using Tr_Base::finite_cells_begin; + using Tr_Base::finite_cells_end; + using Tr_Base::finite_edges_begin; + using Tr_Base::finite_edges_end; + using Tr_Base::incident_facets; + using Tr_Base::insert_in_conflict; + using Tr_Base::infinite_vertex; + using Tr_Base::is_infinite; + using Tr_Base::is_valid; + using Tr_Base::is_valid_finite; + using Tr_Base::locate; + using Tr_Base::mirror_vertex; + using Tr_Base::mirror_index; + using Tr_Base::next_around_edge; + using Tr_Base::number_of_vertices; + using Tr_Base::orientation; + using Tr_Base::point; + using Tr_Base::side_of_segment; + using Tr_Base::side_of_edge; + using Tr_Base::tds; + using Tr_Base::vertex_triple_index; - Regular_triangulation_3(const Gt & gt = Gt(), Lock_data_structure *lock_ds = NULL) - : Tr_Base(gt, lock_ds), hidden_point_visitor(this) - { } + Regular_triangulation_3(const Gt& gt = Gt(), Lock_data_structure *lock_ds = NULL) + : Tr_Base(gt, lock_ds), hidden_point_visitor(this) + { } - Regular_triangulation_3(Lock_data_structure *lock_ds, const Gt & gt = Gt()) - : Tr_Base(lock_ds, gt), hidden_point_visitor(this) - { } + Regular_triangulation_3(Lock_data_structure *lock_ds, const Gt& gt = Gt()) + : Tr_Base(lock_ds, gt), hidden_point_visitor(this) + { } - Regular_triangulation_3(const Regular_triangulation_3 & rt) - : Tr_Base(rt), hidden_point_visitor(this) - { - CGAL_triangulation_postcondition( is_valid() ); - } + Regular_triangulation_3(const Regular_triangulation_3& rt) + : Tr_Base(rt), hidden_point_visitor(this) + { + CGAL_triangulation_postcondition(is_valid()); + } - //insertion - template < typename InputIterator > - Regular_triangulation_3(InputIterator first, InputIterator last, - const Gt & gt = Gt(), Lock_data_structure *lock_ds = NULL) - : Tr_Base(gt, lock_ds), hidden_point_visitor(this) - { - insert(first, last); - } + //insertion + template < typename InputIterator > + Regular_triangulation_3(InputIterator first, InputIterator last, + const Gt& gt = Gt(), Lock_data_structure *lock_ds = NULL) + : Tr_Base(gt, lock_ds), hidden_point_visitor(this) + { + insert(first, last); + } - template < typename InputIterator > - Regular_triangulation_3(InputIterator first, InputIterator last, - Lock_data_structure *lock_ds, const Gt & gt = Gt()) - : Tr_Base(gt, lock_ds), hidden_point_visitor(this) - { - insert(first, last); - } - - private: + template < typename InputIterator > + Regular_triangulation_3(InputIterator first, InputIterator last, + Lock_data_structure *lock_ds, const Gt& gt = Gt()) + : Tr_Base(gt, lock_ds), hidden_point_visitor(this) + { + insert(first, last); + } +private: #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - std::vector - add_temporary_points_on_far_sphere(const size_t num_points) + std::vector + add_temporary_points_on_far_sphere(const size_t num_points) + { + std::vector far_sphere_vertices; + + const size_t MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS = 1000000; + if(num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) { - std::vector far_sphere_vertices; + // Add temporary vertices on a "far sphere" to reduce contention on + // the infinite vertex - const size_t MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS = 1000000; - if (num_points >= MIN_NUM_POINTS_FOR_FAR_SPHERE_POINTS) - { - // Add temporary vertices on a "far sphere" to reduce contention on - // the infinite vertex + // Get bbox + const Bbox_3& bbox = *this->get_bbox(); + // Compute radius for far sphere + const double& xdelta = bbox.xmax() - bbox.xmin(); + const double& ydelta = bbox.ymax() - bbox.ymin(); + const double& zdelta = bbox.zmax() - bbox.zmin(); + const double radius = 1.3 * 0.5 * std::sqrt(xdelta*xdelta + + ydelta*ydelta + + zdelta*zdelta); - // Get bbox - const Bbox_3 &bbox = *this->get_bbox(); - // Compute radius for far sphere - const double& xdelta = bbox.xmax() - bbox.xmin(); - const double& ydelta = bbox.ymax() - bbox.ymin(); - const double& zdelta = bbox.zmax() - bbox.zmin(); - const double radius = 1.3 * 0.5 * std::sqrt(xdelta*xdelta + - ydelta*ydelta + - zdelta*zdelta); - - // WARNING - TODO: this code has to be fixed because Vector_3 is not - // required by the traits concept - const typename Gt::Vector_3 center( - bbox.xmin() + 0.5*xdelta, - bbox.ymin() + 0.5*ydelta, - bbox.zmin() + 0.5*zdelta); - Random_points_on_sphere_3 random_point(radius); - const int NUM_PSEUDO_INFINITE_VERTICES = static_cast( - tbb::task_scheduler_init::default_num_threads() * 3.5); - typename Gt::Construct_weighted_point_3 cwp = + // WARNING - TODO @fixme this code has to be fixed because Vector_3 is not + // required by the traits concept + const typename Gt::Vector_3 center(bbox.xmin() + 0.5*xdelta, + bbox.ymin() + 0.5*ydelta, + bbox.zmin() + 0.5*zdelta); + Random_points_on_sphere_3 random_point(radius); + const int NUM_PSEUDO_INFINITE_VERTICES = static_cast( + tbb::task_scheduler_init::default_num_threads() * 3.5); + typename Gt::Construct_weighted_point_3 cwp = geom_traits().construct_weighted_point_3_object(); - std::vector points_on_far_sphere; - for (int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point) - points_on_far_sphere.push_back(cwp(*random_point + center)); - - // Spatial sorting can only be applied to bare points, so we need an adaptor - typedef typename Geom_traits::Construct_point_3 Construct_point_3; - typedef typename boost::result_of::type Ret; - typedef CGAL::internal::boost_::function_property_map fpmap; - typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; - - spatial_sort(points_on_far_sphere.begin(), points_on_far_sphere.end(), - Search_traits_3( - CGAL::internal::boost_::make_function_property_map( - geom_traits().construct_point_3_object()), geom_traits())); - - typename std::vector::const_iterator it_p = - points_on_far_sphere.begin(); - typename std::vector::const_iterator it_p_end = - points_on_far_sphere.end(); - - for (; it_p != it_p_end ; ++it_p) - { - Locate_type lt; - Cell_handle c, hint; - int li, lj; - - c = locate(*it_p, lt, li, lj, hint); - Vertex_handle v = insert(*it_p, lt, c, li, lj); - hint = (v == Vertex_handle() ? c : v->cell()); - - far_sphere_vertices.push_back(v); - } - } - - return far_sphere_vertices; - } - - void remove_temporary_points_on_far_sphere( - const std::vector & far_sphere_vertices) - { - if (!far_sphere_vertices.empty()) - { - // Remove the temporary vertices on far sphere - remove(far_sphere_vertices.begin(), far_sphere_vertices.end()); - } - } -#endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - - public: - -#ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, InputIterator last, - typename boost::enable_if< - boost::is_convertible< - typename std::iterator_traits::value_type, - Weighted_point - > - >::type* = NULL - ) -#else - template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, InputIterator last) -#endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - { -#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING - static Profile_branch_counter_3 bcounter( - "early withdrawals / late withdrawals / successes [Regular_tri_3::insert]"); -#endif - -#ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; -#endif - - size_type n = number_of_vertices(); - std::vector points(first, last); + std::vector points_on_far_sphere; + for(int i = 0 ; i < NUM_PSEUDO_INFINITE_VERTICES ; ++i, ++random_point) + points_on_far_sphere.push_back(cwp(*random_point + center)); // Spatial sorting can only be applied to bare points, so we need an adaptor - // @todo Unary_function_to_property_map makes a copy (get() returns a value_type) but - // we could hope to get a const & to the bare point. Unfortunately, the lazy - // kernel creates temporaries and prevent it. typedef typename Geom_traits::Construct_point_3 Construct_point_3; typedef typename boost::result_of::type Ret; typedef CGAL::internal::boost_::function_property_map fpmap; typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; - spatial_sort(points.begin(), points.end(), + spatial_sort(points_on_far_sphere.begin(), points_on_far_sphere.end(), Search_traits_3( CGAL::internal::boost_::make_function_property_map( - geom_traits().construct_point_3_object()), geom_traits())); + geom_traits().construct_point_3_object()), geom_traits())); + + typename std::vector::const_iterator it_p = + points_on_far_sphere.begin(); + typename std::vector::const_iterator it_p_end = + points_on_far_sphere.end(); + + for(; it_p != it_p_end ; ++it_p) + { + Locate_type lt; + Cell_handle c, hint; + int li, lj; + + c = locate(*it_p, lt, li, lj, hint); + Vertex_handle v = insert(*it_p, lt, c, li, lj); + hint = (v == Vertex_handle() ? c : v->cell()); + + far_sphere_vertices.push_back(v); + } + } + + return far_sphere_vertices; + } + + void remove_temporary_points_on_far_sphere( + const std::vector& far_sphere_vertices) + { + if(!far_sphere_vertices.empty()) + { + // Remove the temporary vertices on far sphere + remove(far_sphere_vertices.begin(), far_sphere_vertices.end()); + } + } +#endif // CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE + +public: +#ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO + template < class InputIterator > + std::ptrdiff_t insert(InputIterator first, InputIterator last, + typename boost::enable_if< + boost::is_convertible< + typename std::iterator_traits::value_type, + Weighted_point> >::type* = NULL) +#else + template < class InputIterator > + std::ptrdiff_t insert(InputIterator first, InputIterator last) +#endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO + { +#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING + static Profile_branch_counter_3 bcounter( + "early withdrawals / late withdrawals / successes [Regular_tri_3::insert]"); +#endif + +#ifdef CGAL_TRIANGULATION_3_PROFILING + WallClockTimer t; +#endif + + size_type n = number_of_vertices(); + std::vector points(first, last); + + // Spatial sorting can only be applied to bare points, so we need an adaptor + // @todo Unary_function_to_property_map makes a copy (get() returns a value_type) but + // we could hope to get a const& to the bare point. Unfortunately, the lazy + // kernel creates temporaries and prevent it. + typedef typename Geom_traits::Construct_point_3 Construct_point_3; + typedef typename boost::result_of::type Ret; + typedef CGAL::internal::boost_::function_property_map fpmap; + typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; + + spatial_sort(points.begin(), points.end(), + Search_traits_3( + CGAL::internal::boost_::make_function_property_map( + geom_traits().construct_point_3_object()), geom_traits())); // Parallel #ifdef CGAL_LINKED_WITH_TBB - if (this->is_parallel()) - { - size_t num_points = points.size(); - Cell_handle hint; - + if(this->is_parallel()) + { + size_t num_points = points.size(); + Cell_handle hint; + #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - std::vector far_sphere_vertices = + std::vector far_sphere_vertices = add_temporary_points_on_far_sphere(num_points); #endif - - size_t i = 0; - // Insert "num_points_seq" points sequentially - // (or more if dim < 3 after that) - size_t num_points_seq = (std::min)(num_points, (size_t)100); - while (dimension() < 3 || i < num_points_seq) - { - Locate_type lt; - Cell_handle c; - int li, lj; - c = locate (points[i], lt, li, lj, hint); - Vertex_handle v = insert (points[i], lt, c, li, lj); + size_t i = 0; + // Insert "num_points_seq" points sequentially + // (or more if dim < 3 after that) + size_t num_points_seq = (std::min)(num_points, (size_t)100); + while(dimension() < 3 || i < num_points_seq) + { + Locate_type lt; + Cell_handle c; + int li, lj; - hint = (v == Vertex_handle() ? c : v->cell()); - ++i; - } + c = locate(points[i], lt, li, lj, hint); + Vertex_handle v = insert (points[i], lt, c, li, lj); - tbb::enumerable_thread_specific tls_hint(hint->vertex(0)); - tbb::parallel_for( - tbb::blocked_range( i, num_points ), - Insert_point(*this, points, tls_hint) - ); + hint = (v == Vertex_handle() ? c : v->cell()); + ++i; + } + + tbb::enumerable_thread_specific tls_hint(hint->vertex(0)); + tbb::parallel_for(tbb::blocked_range(i, num_points), + Insert_point(*this, points, tls_hint)); #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - remove_temporary_points_on_far_sphere(far_sphere_vertices); + remove_temporary_points_on_far_sphere(far_sphere_vertices); #endif - } - // Sequential - else -#endif // CGAL_LINKED_WITH_TBB - { - Cell_handle hint; - for (typename std::vector::const_iterator p = points.begin(), - end = points.end(); p != end; ++p) - { - Locate_type lt; - Cell_handle c; - int li, lj; - c = locate (*p, lt, li, lj, hint); - - Vertex_handle v = insert (*p, lt, c, li, lj); - - hint = v == Vertex_handle() ? c : v->cell(); - } - } -#ifdef CGAL_TRIANGULATION_3_PROFILING - std::cerr << "Points inserted in " << t.elapsed() << " seconds." << std::endl; -#endif - return number_of_vertices() - n; } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB + { + Cell_handle hint; + for(typename std::vector::const_iterator p = points.begin(), + end = points.end(); p != end; ++p) + { + Locate_type lt; + Cell_handle c; + int li, lj; + c = locate(*p, lt, li, lj, hint); + Vertex_handle v = insert(*p, lt, c, li, lj); + + hint = v == Vertex_handle() ? c : v->cell(); + } + } +#ifdef CGAL_TRIANGULATION_3_PROFILING + std::cerr << "Points inserted in " << t.elapsed() << " seconds." << std::endl; +#endif + return number_of_vertices() - n; + } #ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - private: - //top stands for tuple-or-pair - template - const Weighted_point& top_get_first(const std::pair& pair) const { return pair.first; } - template - const Info& top_get_second(const std::pair& pair) const { return pair.second; } - template - const Weighted_point& top_get_first(const boost::tuple& tuple) const { return boost::get<0>(tuple); } - template - const Info& top_get_second(const boost::tuple& tuple) const { return boost::get<1>(tuple); } +private: + //top stands for tuple-or-pair + template + const Weighted_point& top_get_first(const std::pair& pair) const { return pair.first; } - // Functor to go from an index of a container of Weighted_point to - // the corresponding Bare_point - template - struct Index_to_Bare_point + template + const Info& top_get_second(const std::pair& pair) const { return pair.second; } + + template + const Weighted_point& top_get_first(const boost::tuple& tuple) const { return boost::get<0>(tuple); } + + template + const Info& top_get_second(const boost::tuple& tuple) const { return boost::get<1>(tuple); } + + // Functor to go from an index of a container of Weighted_point to + // the corresponding Bare_point + template + struct Index_to_Bare_point + { + const Bare_point& operator()(const std::size_t& i) const { - const Bare_point& operator()(const std::size_t& i) const - { - return cp(c[i]); - } + return cp(c[i]); + } - Index_to_Bare_point(const Container& c, const Construct_bare_point& cp) - : c(c), cp(cp) { } + Index_to_Bare_point(const Container& c, const Construct_bare_point& cp) + : c(c), cp(cp) { } - const Container& c; - const Construct_bare_point cp; - }; + const Container& c; + const Construct_bare_point cp; + }; - template - std::ptrdiff_t insert_with_info(InputIterator first,InputIterator last) + template + std::ptrdiff_t insert_with_info(InputIterator first,InputIterator last) + { + size_type n = number_of_vertices(); + std::vector indices; + std::vector points; + std::vector infos; + std::size_t index=0; + for(InputIterator it=first;it!=last;++it) { - size_type n = number_of_vertices(); - std::vector indices; - std::vector points; - std::vector infos; - std::size_t index=0; - for (InputIterator it=first;it!=last;++it){ - Tuple_or_pair pair = *it; - points.push_back( top_get_first(pair) ); - infos.push_back ( top_get_second(pair) ); - indices.push_back(index++); - } + Tuple_or_pair pair = *it; + points.push_back(top_get_first(pair)); + infos.push_back(top_get_second(pair)); + indices.push_back(index++); + } - // We need to sort the points and their info at the same time through - // the `indices` vector AND spatial sort can only handle Gt::Point_3. - typedef typename Geom_traits::Construct_point_3 Construct_point_3; - typedef Index_to_Bare_point > Access_bare_point; - typedef typename boost::result_of::type Ret; - typedef CGAL::internal::boost_::function_property_map fpmap; - typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; + // We need to sort the points and their info at the same time through + // the `indices` vector AND spatial sort can only handle Gt::Point_3. + typedef typename Geom_traits::Construct_point_3 Construct_point_3; + typedef Index_to_Bare_point > Access_bare_point; + typedef typename boost::result_of::type Ret; + typedef CGAL::internal::boost_::function_property_map fpmap; + typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; - Access_bare_point accessor(points, geom_traits().construct_point_3_object()); - spatial_sort(indices.begin(), indices.end(), - Search_traits_3( - CGAL::internal::boost_::make_function_property_map< - std::size_t, Ret, Access_bare_point>(accessor), - geom_traits())); + Access_bare_point accessor(points, geom_traits().construct_point_3_object()); + spatial_sort(indices.begin(), indices.end(), + Search_traits_3( + CGAL::internal::boost_::make_function_property_map< + std::size_t, Ret, Access_bare_point>(accessor), + geom_traits())); #ifdef CGAL_LINKED_WITH_TBB - if (this->is_parallel()) - { - size_t num_points = points.size(); - Cell_handle hint; + if(this->is_parallel()) + { + size_t num_points = points.size(); + Cell_handle hint; #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - std::vector far_sphere_vertices = + std::vector far_sphere_vertices = add_temporary_points_on_far_sphere(num_points); #endif - size_t i = 0; - // Insert "num_points_seq" points sequentially - // (or more if dim < 3 after that) - size_t num_points_seq = (std::min)(num_points, (size_t)100); - while (dimension() < 3 || i < num_points_seq) + size_t i = 0; + // Insert "num_points_seq" points sequentially + // (or more if dim < 3 after that) + size_t num_points_seq = (std::min)(num_points, (size_t)100); + while(dimension() < 3 || i < num_points_seq) + { + Locate_type lt; + Cell_handle c; + int li, lj; + c = locate(points[indices[i]], lt, li, lj, hint); + + Vertex_handle v = insert(points[indices[i]], lt, c, li, lj); + if(v != Vertex_handle()) { - Locate_type lt; - Cell_handle c; - int li, lj; - c = locate(points[indices[i]], lt, li, lj, hint); - - Vertex_handle v = insert(points[indices[i]], lt, c, li, lj); - if (v != Vertex_handle()){ - v->info() = infos[indices[i]]; - hint = v->cell(); - } - else - hint = c; - - ++i; + v->info() = infos[indices[i]]; + hint = v->cell(); } + else + hint = c; - tbb::enumerable_thread_specific tls_hint(hint->vertex(0)); - tbb::parallel_for( - tbb::blocked_range(i, num_points), - Insert_point_with_info(*this, points, infos, indices, tls_hint) - ); + ++i; + } + + tbb::enumerable_thread_specific tls_hint(hint->vertex(0)); + tbb::parallel_for(tbb::blocked_range(i, num_points), + Insert_point_with_info(*this, points, infos, indices, tls_hint)); #ifdef CGAL_CONCURRENT_TRIANGULATION_3_ADD_TEMPORARY_POINTS_ON_FAR_SPHERE - remove_temporary_points_on_far_sphere(far_sphere_vertices); + remove_temporary_points_on_far_sphere(far_sphere_vertices); #endif - } - // Sequential - else + } + // Sequential + else #endif // CGAL_LINKED_WITH_TBB + { + Cell_handle hint; + for(typename std::vector::const_iterator + it = indices.begin(), end = indices.end(); + it != end; ++it) { - Cell_handle hint; - for (typename std::vector::const_iterator - it = indices.begin(), end = indices.end(); - it != end; ++it) - { - Locate_type lt; - Cell_handle c; - int li, lj; - c = locate (points[*it], lt, li, lj, hint); + Locate_type lt; + Cell_handle c; + int li, lj; + c = locate(points[*it], lt, li, lj, hint); - Vertex_handle v = insert (points[*it], lt, c, li, lj); - if (v!=Vertex_handle()){ - v->info()=infos[*it]; - hint=v->cell(); - } - else - hint=c; + Vertex_handle v = insert(points[*it], lt, c, li, lj); + if(v!=Vertex_handle()) + { + v->info()=infos[*it]; + hint=v->cell(); + } + else + { + hint = c; } } - - return number_of_vertices() - n; } - public: + return number_of_vertices() - n; + } - template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, - InputIterator last, - typename boost::enable_if< - boost::is_convertible< - typename std::iterator_traits::value_type, - std::pair::type> - > - >::type* = NULL - ) - {return insert_with_info< std::pair::type> >(first,last);} +public: - template - std::ptrdiff_t - insert( boost::zip_iterator< boost::tuple > first, - boost::zip_iterator< boost::tuple > last, - typename boost::enable_if< - boost::mpl::and_< - typename boost::is_convertible< typename std::iterator_traits::value_type, Weighted_point >, - typename boost::is_convertible< typename std::iterator_traits::value_type, typename internal::Info_check::type > - > - >::type* =NULL - ) - {return insert_with_info< boost::tuple::type> >(first,last);} + template < class InputIterator > + std::ptrdiff_t insert(InputIterator first, + InputIterator last, + typename boost::enable_if< + boost::is_convertible< + typename std::iterator_traits::value_type, + std::pair::type> + > + >::type* = NULL) + { + return insert_with_info< + std::pair::type> + >(first,last); + } + + template + std::ptrdiff_t + insert(boost::zip_iterator< boost::tuple > first, + boost::zip_iterator< boost::tuple > last, + typename boost::enable_if< + boost::mpl::and_< + typename boost::is_convertible< typename std::iterator_traits::value_type, Weighted_point >, + typename boost::is_convertible< typename std::iterator_traits::value_type, typename internal::Info_check::type > + > >::type* =NULL) + { + return insert_with_info< + boost::tuple::type> + >(first,last); + } #endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - Vertex_handle insert(const Weighted_point & p, Vertex_handle hint, - bool *could_lock_zone = NULL) + Vertex_handle insert(const Weighted_point& p, Vertex_handle hint, + bool *could_lock_zone = NULL) + { + return insert(p, + hint == Vertex_handle() ? this->infinite_cell() : hint->cell(), + could_lock_zone); + } + + Vertex_handle insert(const Weighted_point& p, + Cell_handle start = Cell_handle(), bool *could_lock_zone = NULL); + + Vertex_handle insert(const Weighted_point& p, Locate_type lt, + Cell_handle c, int li, int, bool *could_lock_zone = NULL); + + template + Vertex_handle insert_in_hole(const Weighted_point& p, + CellIt cell_begin, CellIt cell_end, + Cell_handle begin, int i); + + template + Vertex_handle insert_in_hole(const Weighted_point& p, + CellIt cell_begin, CellIt cell_end, + Cell_handle begin, int i, Vertex_handle newv); + + template + Triple + find_conflicts(const Weighted_point& p, Cell_handle c, + OutputIteratorBoundaryFacets bfit, + OutputIteratorCells cit, + OutputIteratorInternalFacets ifit, + bool *could_lock_zone = NULL, + const Facet *this_facet_must_be_in_the_cz = NULL, + bool *the_facet_is_in_its_cz = NULL) const + { + CGAL_triangulation_precondition(dimension() >= 2); + + std::vector cells; + cells.reserve(32); + std::vector facets; + facets.reserve(64); + + if(dimension() == 2) { - return insert(p, - hint == Vertex_handle() ? this->infinite_cell() : hint->cell(), - could_lock_zone); + Conflict_tester_2 tester(p, this); + if(! tester(c)) + return make_triple(bfit, cit, ifit); + + ifit = Tr_Base::find_conflicts(c, tester, + make_triple(std::back_inserter(facets), + std::back_inserter(cells), + ifit), + could_lock_zone, + this_facet_must_be_in_the_cz, + the_facet_is_in_its_cz).third; + } + else + { + Conflict_tester_3 tester(p, this); + if(! tester(c)) + return make_triple(bfit, cit, ifit); + + ifit = Tr_Base::find_conflicts(c, tester, + make_triple(std::back_inserter(facets), + std::back_inserter(cells), + ifit), + could_lock_zone, + this_facet_must_be_in_the_cz, + the_facet_is_in_its_cz).third; } - Vertex_handle insert(const Weighted_point & p, - Cell_handle start = Cell_handle(), bool *could_lock_zone = NULL); - - Vertex_handle insert(const Weighted_point & p, Locate_type lt, - Cell_handle c, int li, int, bool *could_lock_zone = NULL); - - template - Vertex_handle - insert_in_hole(const Weighted_point & p, CellIt cell_begin, CellIt cell_end, - Cell_handle begin, int i); - - template - Vertex_handle - insert_in_hole(const Weighted_point & p, CellIt cell_begin, CellIt cell_end, - Cell_handle begin, int i, Vertex_handle newv); - - template - Triple - find_conflicts(const Weighted_point &p, Cell_handle c, - OutputIteratorBoundaryFacets bfit, - OutputIteratorCells cit, - OutputIteratorInternalFacets ifit - , bool *could_lock_zone = NULL - , const Facet *this_facet_must_be_in_the_cz = NULL - , bool *the_facet_is_in_its_cz = NULL - ) const + // Reset the conflict flag on the boundary. + for(typename std::vector::iterator fit = facets.begin(); + fit != facets.end(); ++fit) { - CGAL_triangulation_precondition(dimension() >= 2); - - std::vector cells; - cells.reserve(32); - std::vector facets; - facets.reserve(64); - - if (dimension() == 2) { - Conflict_tester_2 tester(p, this); - if (! tester (c)) return make_triple (bfit, cit, ifit); - ifit = Tr_Base::find_conflicts - (c, tester, - make_triple(std::back_inserter(facets), - std::back_inserter(cells), - ifit) - , could_lock_zone - , this_facet_must_be_in_the_cz - , the_facet_is_in_its_cz - ).third; - } - else { - Conflict_tester_3 tester(p, this); - if (! tester (c)) return make_triple (bfit, cit, ifit); - ifit = Tr_Base::find_conflicts - (c, tester, - make_triple(std::back_inserter(facets), - std::back_inserter(cells), - ifit) - , could_lock_zone - , this_facet_must_be_in_the_cz - , the_facet_is_in_its_cz - ).third; - } - - // Reset the conflict flag on the boundary. - for(typename std::vector::iterator fit=facets.begin(); - fit != facets.end(); ++fit) { - fit->first->neighbor(fit->second)->tds_data().clear(); - *bfit++ = *fit; - } - - // Reset the conflict flag in the conflict cells. - for(typename std::vector::iterator ccit=cells.begin(); - ccit != cells.end(); ++ccit) { - (*ccit)->tds_data().clear(); - *cit++ = *ccit; - } - return make_triple(bfit, cit, ifit); + fit->first->neighbor(fit->second)->tds_data().clear(); + *bfit++ = *fit; } - template - std::pair - find_conflicts(const Weighted_point &p, Cell_handle c, - OutputIteratorBoundaryFacets bfit, - OutputIteratorCells cit - , bool *could_lock_zone = NULL - ) const + // Reset the conflict flag in the conflict cells. + for(typename std::vector::iterator ccit = cells.begin(); + ccit != cells.end(); ++ccit) { - Triple t = find_conflicts(p, c, bfit, cit, - Emptyset_iterator(), - could_lock_zone); - return std::make_pair(t.first, t.second); + (*ccit)->tds_data().clear(); + *cit++ = *ccit; } + return make_triple(bfit, cit, ifit); + } - // Returns the vertices on the interior of the conflict hole. - template - OutputIterator - vertices_inside_conflict_zone(const Weighted_point&p, Cell_handle c, - OutputIterator res) const - { - CGAL_triangulation_precondition(dimension() >= 2); + template + std::pair + find_conflicts(const Weighted_point& p, Cell_handle c, + OutputIteratorBoundaryFacets bfit, + OutputIteratorCells cit, + bool *could_lock_zone = NULL) const + { + Triple t = find_conflicts(p, c, bfit, cit, + Emptyset_iterator(), + could_lock_zone); + return std::make_pair(t.first, t.second); + } - // Get the facets on the boundary of the hole, and the cells of the hole - std::vector cells; - std::vector facets; - find_conflicts(p, c, std::back_inserter(facets), - std::back_inserter(cells), Emptyset_iterator()); + // Returns the vertices on the interior of the conflict hole. + template + OutputIterator vertices_inside_conflict_zone(const Weighted_point&p, Cell_handle c, + OutputIterator res) const + { + CGAL_triangulation_precondition(dimension() >= 2); - // Put all vertices on the hole in 'vertices' - const int d = dimension(); - std::set vertices; - for (typename std::vector::const_iterator - it = cells.begin(), + // Get the facets on the boundary of the hole, and the cells of the hole + std::vector cells; + std::vector facets; + find_conflicts(p, c, std::back_inserter(facets), + std::back_inserter(cells), Emptyset_iterator()); + + // Put all vertices on the hole in 'vertices' + const int d = dimension(); + std::set vertices; + for(typename std::vector::const_iterator it = cells.begin(), end = cells.end(); it != end; ++it) - { - for(int i = 0; i <= d; ++i) { - vertices.insert((*it)->vertex(i)); - } - } - // Then extract the vertices of the boundary and remove them from - // 'vertices' - if (dimension() == 3) { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.erase(i->first->vertex((i->second+1)&3)); - vertices.erase(i->first->vertex((i->second+2)&3)); - vertices.erase(i->first->vertex((i->second+3)&3)); - } - } else { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.erase(i->first->vertex(cw(i->second))); - vertices.erase(i->first->vertex(ccw(i->second))); - } - } - - return std::copy(vertices.begin(), vertices.end(), res); + { + for(int i = 0; i <= d; ++i) + vertices.insert((*it)->vertex(i)); } + // Then extract the vertices of the boundary and remove them from + // 'vertices' + if(dimension() == 3) + { + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.erase(i->first->vertex((i->second+1)&3)); + vertices.erase(i->first->vertex((i->second+2)&3)); + vertices.erase(i->first->vertex((i->second+3)&3)); + } + } + else + { + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.erase(i->first->vertex(cw(i->second))); + vertices.erase(i->first->vertex(ccw(i->second))); + } + } + + return std::copy(vertices.begin(), vertices.end(), res); + } #ifndef CGAL_NO_DEPRECATED_CODE - // Returns the vertices on the boundary of the conflict hole. - template - OutputIterator - vertices_in_conflict(const Weighted_point&p, Cell_handle c, OutputIterator res) const - { - return vertices_on_conflict_zone_boundary(p, c, res); - } + // Returns the vertices on the boundary of the conflict hole. + template + OutputIterator vertices_in_conflict(const Weighted_point&p, Cell_handle c, + OutputIterator res) const + { + return vertices_on_conflict_zone_boundary(p, c, res); + } #endif // CGAL_NO_DEPRECATED_CODE - // Returns the vertices on the boundary of the conflict hole. - template - OutputIterator - vertices_on_conflict_zone_boundary(const Weighted_point&p, Cell_handle c, - OutputIterator res) const + // Returns the vertices on the boundary of the conflict hole. + template + OutputIterator vertices_on_conflict_zone_boundary(const Weighted_point&p, + Cell_handle c, + OutputIterator res) const + { + CGAL_triangulation_precondition(dimension() >= 2); + + // Get the facets on the boundary of the hole. + std::vector facets; + find_conflicts(p, c, std::back_inserter(facets), + Emptyset_iterator(), Emptyset_iterator()); + + // Then extract uniquely the vertices. + std::set vertices; + if(dimension() == 3) { - CGAL_triangulation_precondition(dimension() >= 2); - - // Get the facets on the boundary of the hole. - std::vector facets; - find_conflicts(p, c, std::back_inserter(facets), - Emptyset_iterator(), Emptyset_iterator()); - - // Then extract uniquely the vertices. - std::set vertices; - if (dimension() == 3) { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.insert(i->first->vertex((i->second+1)&3)); - vertices.insert(i->first->vertex((i->second+2)&3)); - vertices.insert(i->first->vertex((i->second+3)&3)); - } - } else { - for (typename std::vector::const_iterator i = facets.begin(); - i != facets.end(); ++i) { - vertices.insert(i->first->vertex(cw(i->second))); - vertices.insert(i->first->vertex(ccw(i->second))); - } + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.insert(i->first->vertex((i->second+1)&3)); + vertices.insert(i->first->vertex((i->second+2)&3)); + vertices.insert(i->first->vertex((i->second+3)&3)); + } + } else { + for(typename std::vector::const_iterator i = facets.begin(); + i != facets.end(); ++i) + { + vertices.insert(i->first->vertex(cw(i->second))); + vertices.insert(i->first->vertex(ccw(i->second))); } - - return std::copy(vertices.begin(), vertices.end(), res); } - void remove (Vertex_handle v); - // Concurrency-safe - // See Triangulation_3::remove for more information - bool remove (Vertex_handle v, bool *could_lock_zone); + return std::copy(vertices.begin(), vertices.end(), res); + } - template < typename InputIterator > - size_type remove(InputIterator first, InputIterator beyond) - { - CGAL_triangulation_precondition(!this->does_repeat_in_range(first, beyond)); - size_type n = number_of_vertices(); + void remove(Vertex_handle v); + // Concurrency-safe + // See Triangulation_3::remove for more information + bool remove(Vertex_handle v, bool *could_lock_zone); + + template < typename InputIterator > + size_type remove(InputIterator first, InputIterator beyond) + { + CGAL_triangulation_precondition(!this->does_repeat_in_range(first, beyond)); + size_type n = number_of_vertices(); #ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; + WallClockTimer t; #endif - // Parallel + // Parallel #ifdef CGAL_LINKED_WITH_TBB - if (this->is_parallel()) - { - // TODO: avoid that by asking for ramdom-access iterators? - std::vector vertices(first, beyond); - tbb::concurrent_vector vertices_to_remove_sequentially; + if(this->is_parallel()) + { + // TODO: avoid that by asking for ramdom-access iterators? + std::vector vertices(first, beyond); + tbb::concurrent_vector vertices_to_remove_sequentially; - tbb::parallel_for( - tbb::blocked_range( 0, vertices.size()), - Remove_point(*this, vertices, vertices_to_remove_sequentially) - ); + tbb::parallel_for(tbb::blocked_range(0, vertices.size()), + Remove_point(*this, vertices, vertices_to_remove_sequentially)); - // Do the rest sequentially - for ( typename tbb::concurrent_vector::const_iterator - it = vertices_to_remove_sequentially.begin(), - it_end = vertices_to_remove_sequentially.end() + // Do the rest sequentially + for(typename tbb::concurrent_vector::const_iterator + it = vertices_to_remove_sequentially.begin(), + it_end = vertices_to_remove_sequentially.end() ; it != it_end ; ++it) - { - remove(*it); - } - } - // Sequential - else -#endif // CGAL_LINKED_WITH_TBB { - while (first != beyond) { - remove(*first); - ++first; - } + remove(*it); } + } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB + { + while(first != beyond) + { + remove(*first); + ++first; + } + } #ifdef CGAL_TRIANGULATION_3_PROFILING - std::cerr << "Points removed in " << t.elapsed() << " seconds." << std::endl; + std::cerr << "Points removed in " << t.elapsed() << " seconds." << std::endl; #endif - return n - number_of_vertices(); - } + return n - number_of_vertices(); + } - - template - void remove_and_give_new_cells(Vertex_handle v, OutputItCells cit) + template + void remove_and_give_new_cells(Vertex_handle v, OutputItCells cit) + { + Self tmp; + Vertex_remover remover(tmp); + Tr_Base::remove_and_give_new_cells(v, remover, cit); + + CGAL_triangulation_expensive_postcondition(is_valid()); + } + + // DISPLACEMENT + Vertex_handle move_point(Vertex_handle v, const Weighted_point& p); + + // Displacement works only for regular triangulation + // without hidden points at any time + Vertex_handle move_if_no_collision(Vertex_handle v, const Weighted_point& p); + Vertex_handle move(Vertex_handle v, const Weighted_point& p); + + // REMOVE CLUSTER - works only when Regular has no hidden point at all + // "regular as Delaunay" + template < typename InputIterator > + size_type remove_cluster(InputIterator first, InputIterator beyond) + { + Self tmp; + Vertex_remover remover(tmp); + return Tr_Base::remove(first, beyond, remover); + } + +protected: + + Oriented_side side_of_oriented_power_sphere(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p3, + const Weighted_point& p, + bool perturb = false) const; + + Oriented_side side_of_oriented_power_circle(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p, + bool perturb = false) const; + + Bounded_side side_of_bounded_power_circle(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p, + bool perturb = false) const; + + Bounded_side side_of_bounded_power_segment(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p, + bool perturb = false) const; + +public: + // Queries + Bounded_side side_of_power_sphere(Cell_handle c, const Weighted_point& p, + bool perturb = false) const; + + Bounded_side side_of_power_circle(const Facet& f, const Weighted_point& p, + bool /* perturb */ = false) const + { + return side_of_power_circle(f.first, f.second, p); + } + + Bounded_side side_of_power_circle(Cell_handle c, int i, const Weighted_point& p, + bool perturb = false) const; + + Bounded_side side_of_power_segment(Cell_handle c, const Weighted_point& p, + bool perturb = false) const; + + // Undocumented, needed for Mesh_3 (because of Periodic_3_mesh_3) + bool greater_or_equal_power_distance(const Bare_point& p, + const Weighted_point& q, + const Weighted_point& r) const; + + double compute_power_distance_to_power_sphere(const Cell_handle c, + const int i) const; + double compute_power_distance_to_power_sphere(const Cell_handle c, + const Vertex_handle v) const; + + Vertex_handle nearest_power_vertex_in_cell(const Bare_point& p, + Cell_handle c) const; + + Vertex_handle nearest_power_vertex(const Bare_point& p, Cell_handle c = Cell_handle()) const; + + bool is_Gabriel(Cell_handle c, int i) const; + bool is_Gabriel(Cell_handle c, int i, int j) const; + bool is_Gabriel(const Facet& f)const ; + bool is_Gabriel(const Edge& e) const; + bool is_Gabriel(Vertex_handle v) const; + + // Dual functions + Bare_point dual(Cell_handle c) const; + Object dual(Cell_handle c, int i) const; + Object dual(const Facet& facet) const; + + void dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const; + void dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const; + void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const; + void dual_ray(Cell_handle c, int i, Ray& ray) const; + void dual_ray(const Facet& facet, Ray& ray) const; + void dual_ray_exact(const Facet& facet, Ray& ray) const; + + template < class Stream> + Stream& draw_dual(Stream& os) const; + + bool is_valid(bool verbose = false, int level = 0) const; + +protected: + bool less_power_distance(const Bare_point& p, + const Weighted_point& q, + const Weighted_point& r) const + { + return geom_traits().compare_power_distance_3_object()(p, q, r) == SMALLER; + } + + Bare_point construct_weighted_circumcenter(const Weighted_point& p, + const Weighted_point& q, + const Weighted_point& r, + const Weighted_point& s) const + { + return geom_traits().construct_weighted_circumcenter_3_object()(p,q,r,s); + } + + Bare_point construct_weighted_circumcenter(const Weighted_point& p, + const Weighted_point& q, + const Weighted_point& r) const + { + return geom_traits().construct_weighted_circumcenter_3_object()(p,q,r); + } + + Segment construct_segment(const Bare_point& p, const Bare_point& q) const + { + return geom_traits().construct_segment_3_object()(p, q); + } + + Line construct_perpendicular_line(const Plane& pl, const Bare_point& p) const + { + return geom_traits().construct_perpendicular_line_3_object()(pl, p); + } + + Plane construct_plane(const Bare_point& p, const Bare_point& q, const Bare_point& r) const + { + return geom_traits().construct_plane_3_object()(p, q, r); + } + + Ray construct_ray(const Bare_point& p, const Line& l) const + { + return geom_traits().construct_ray_3_object()(p, l); + } + + Object construct_object(const Bare_point& p) const + { + return geom_traits().construct_object_3_object()(p); + } + + Object construct_object(const Segment& s) const + { + return geom_traits().construct_object_3_object()(s); + } + + Object construct_object(const Ray& r) const + { + return geom_traits().construct_object_3_object()(r); + } + + Vertex_handle nearest_power_vertex(const Bare_point& p, + Vertex_handle v, + Vertex_handle w) const + { + // In case of equality, v is returned. + CGAL_triangulation_precondition(v != w); + if(is_infinite(v)) + return w; + + if(is_infinite(w)) + return v; + + return less_power_distance(p, w->point(), v->point()) ? w : v; + } + + Oriented_side power_test(const Weighted_point& p, const Weighted_point& q) const + { + CGAL_triangulation_precondition(this->equal(p, q)); + return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q); + } + + Oriented_side power_test(const Weighted_point& p, const Weighted_point& q, + const Weighted_point& r) const + { + CGAL_triangulation_precondition(this->collinear(p, q, r)); + return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r); + } + + Oriented_side power_test(const Weighted_point& p, const Weighted_point& q, + const Weighted_point& r, const Weighted_point& s) const + { + CGAL_triangulation_precondition(this->coplanar(p, q, r, s)); + return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r, s); + } + + Oriented_side power_test(const Weighted_point& p, const Weighted_point& q, + const Weighted_point& r, const Weighted_point& s, + const Weighted_point& t) const + { + return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r, s, t); + } + + bool in_conflict_3(const Weighted_point& p, const Cell_handle c) const + { + return side_of_power_sphere(c, p, true) == ON_BOUNDED_SIDE; + } + + bool in_conflict_2(const Weighted_point& p, const Cell_handle c, int i) const + { + return side_of_power_circle(c, i, p, true) == ON_BOUNDED_SIDE; + } + + bool in_conflict_1(const Weighted_point& p, const Cell_handle c) const + { + return side_of_power_segment(c, p, true) == ON_BOUNDED_SIDE; + } + + bool in_conflict_0(const Weighted_point& p, const Cell_handle c) const + { + return power_test(c->vertex(0)->point(), p) == ON_POSITIVE_SIDE; + } + + bool in_conflict(const Weighted_point& p, const Cell_handle c) const + { + switch(dimension()) { - Self tmp; - Vertex_remover remover (tmp); - Tr_Base::remove_and_give_new_cells(v, remover, cit); - - CGAL_triangulation_expensive_postcondition(is_valid()); + case 0: + return in_conflict_0(p, c); + case 1: + return in_conflict_1(p, c); + case 2: + return in_conflict_2(p, c, 3); + case 3: + return in_conflict_3(p, c); } + return true; + } - - // DISPLACEMENT - Vertex_handle move_point(Vertex_handle v, const Weighted_point & p); - - // Displacement works only for regular triangulation - // without hidden points at any time - Vertex_handle move_if_no_collision(Vertex_handle v, const Weighted_point & p); - Vertex_handle move(Vertex_handle v, const Weighted_point & p); - - // REMOVE CLUSTER - works only when Regular has no hidden point at all - // "regular as Delaunay" - template < typename InputIterator > - size_type remove_cluster(InputIterator first, InputIterator beyond) - { - Self tmp; - Vertex_remover remover (tmp); - return Tr_Base::remove(first, beyond, remover); - } - - protected: - - Oriented_side - side_of_oriented_power_sphere(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p3, - const Weighted_point &p, - bool perturb = false) const; - - Oriented_side - side_of_oriented_power_circle(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p, - bool perturb = false) const; - - Bounded_side - side_of_bounded_power_circle(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p, - bool perturb = false) const; - - Bounded_side - side_of_bounded_power_segment(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p, - bool perturb = false) const; - + class Conflict_tester_3 + { + const Weighted_point& p; + const Self *t; public: - - // Queries - Bounded_side - side_of_power_sphere(Cell_handle c, const Weighted_point &p, - bool perturb = false) const; - - Bounded_side - side_of_power_circle(const Facet & f, const Weighted_point & p, - bool /* perturb */ = false) const + Conflict_tester_3(const Weighted_point& pt, const Self *tr) + : p(pt), t(tr) { - return side_of_power_circle(f.first, f.second, p); } - Bounded_side - side_of_power_circle(Cell_handle c, int i, const Weighted_point &p, - bool perturb = false) const; - - Bounded_side - side_of_power_segment(Cell_handle c, const Weighted_point &p, - bool perturb = false) const; - - // Undocumented, needed for Mesh_3 (because of Periodic_3_mesh_3) - bool greater_or_equal_power_distance(const Bare_point& p, - const Weighted_point& q, - const Weighted_point& r) const; - - double compute_power_distance_to_power_sphere(const Cell_handle c, - const int i) const; - double compute_power_distance_to_power_sphere(const Cell_handle c, - const Vertex_handle v) const; - - Vertex_handle - nearest_power_vertex_in_cell(const Bare_point& p, - Cell_handle c) const; - - Vertex_handle - nearest_power_vertex(const Bare_point& p, Cell_handle c = - Cell_handle()) const; - - bool is_Gabriel(Cell_handle c, int i) const; - bool is_Gabriel(Cell_handle c, int i, int j) const; - bool is_Gabriel(const Facet& f)const ; - bool is_Gabriel(const Edge& e) const; - bool is_Gabriel(Vertex_handle v) const; - - - // Dual functions - Bare_point dual(Cell_handle c) const; - Object dual(Cell_handle c, int i) const; - Object dual(const Facet& facet) const; - - void dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const; - void dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const; - void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const; - void dual_ray(Cell_handle c, int i, Ray& ray) const; - void dual_ray(const Facet& facet, Ray& ray) const; - void dual_ray_exact(const Facet& facet, Ray& ray) const; - - template < class Stream> - Stream& draw_dual(Stream & os) const; - - bool is_valid(bool verbose = false, int level = 0) const; - - protected: - bool - less_power_distance(const Bare_point &p, - const Weighted_point &q, - const Weighted_point &r) const + bool operator()(const Cell_handle c) const { - return - geom_traits().compare_power_distance_3_object()(p, q, r) == SMALLER; + return t->in_conflict_3(p, c); } - Bare_point - construct_weighted_circumcenter(const Weighted_point &p, - const Weighted_point &q, - const Weighted_point &r, - const Weighted_point &s) const + bool test_initial_cell(const Cell_handle c) const + { + return operator()(c); + } + Oriented_side compare_weight(const Weighted_point& wp1, + const Weighted_point& wp2) const + { + return t->power_test(wp1, wp2); + } + }; + + class Conflict_tester_2 + { + const Weighted_point& p; + const Self *t; + + public: + Conflict_tester_2(const Weighted_point& pt, const Self *tr) + : p(pt), t(tr) { - return geom_traits().construct_weighted_circumcenter_3_object()(p,q,r,s); } - Bare_point - construct_weighted_circumcenter(const Weighted_point &p, - const Weighted_point &q, - const Weighted_point &r) const + bool operator()(const Cell_handle c) const { - return geom_traits().construct_weighted_circumcenter_3_object()(p,q,r); + return t->in_conflict_2(p, c, 3); } - Segment - construct_segment(const Bare_point &p, const Bare_point &q) const + bool test_initial_cell(const Cell_handle c) const { - return geom_traits().construct_segment_3_object()(p, q); + return operator()(c); } - Line - construct_perpendicular_line(const Plane &pl, const Bare_point &p) const + Oriented_side compare_weight(const Weighted_point& wp1, + const Weighted_point& wp2) const + { + return t->power_test(wp1, wp2); + } + }; + + class Conflict_tester_1 + { + const Weighted_point& p; + const Self *t; + + public: + Conflict_tester_1(const Weighted_point& pt, const Self *tr) + : p(pt), t(tr) { - return geom_traits().construct_perpendicular_line_3_object()(pl, p); } - Plane - construct_plane(const Bare_point &p, const Bare_point &q, const Bare_point &r) const + bool operator()(const Cell_handle c) const { - return geom_traits().construct_plane_3_object()(p, q, r); + return t->in_conflict_1(p, c); } - Ray - construct_ray(const Bare_point &p, const Line &l) const + bool test_initial_cell(const Cell_handle c) const { - return geom_traits().construct_ray_3_object()(p, l); + return operator()(c); } - Object - construct_object(const Bare_point &p) const + Oriented_side compare_weight(const Weighted_point& wp1, + const Weighted_point& wp2) const + { + return t->power_test(wp1, wp2); + } + }; + + class Conflict_tester_0 + { + const Weighted_point& p; + const Self *t; + + public: + Conflict_tester_0(const Weighted_point& pt, const Self *tr) + : p(pt), t(tr) { - return geom_traits().construct_object_3_object()(p); } - Object - construct_object(const Segment &s) const + bool operator()(const Cell_handle c) const { - return geom_traits().construct_object_3_object()(s); + return t->in_conflict_0(p, c); } - Object - construct_object(const Ray &r) const + bool test_initial_cell(const Cell_handle c) const { - return geom_traits().construct_object_3_object()(r); + return operator()(c); } - Vertex_handle - nearest_power_vertex(const Bare_point &p, - Vertex_handle v, - Vertex_handle w) const + int compare_weight(const Weighted_point& wp1, + const Weighted_point& wp2) const { - // In case of equality, v is returned. - CGAL_triangulation_precondition(v != w); - if (is_infinite(v)) return w; - if (is_infinite(w)) return v; - return less_power_distance(p, w->point(), v->point()) ? w : v; + return t->power_test(wp1, wp2); } + }; - Oriented_side - power_test(const Weighted_point &p, const Weighted_point &q) const + // Sequential version + // "dummy" is here to allow the specialization (see below) + // See http://groups.google.com/group/comp.lang.c++.moderated/browse_thread/thread/285ab1eec49e1cb6 + template + class Hidden_point_visitor + { + Self *t; + mutable std::vector vertices; + mutable std::vector hidden_points; + + public: + Hidden_point_visitor(Self *tr) : t(tr) {} + + template + void process_cells_in_conflict(InputIterator start, InputIterator end) const { - CGAL_triangulation_precondition(this->equal(p, q)); - return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q); - } - - Oriented_side - power_test(const Weighted_point &p, const Weighted_point &q, - const Weighted_point &r) const - { - CGAL_triangulation_precondition(this->collinear(p, q, r)); - return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r); - } - - Oriented_side - power_test(const Weighted_point &p, const Weighted_point &q, - const Weighted_point &r, const Weighted_point &s) const - { - CGAL_triangulation_precondition(this->coplanar(p, q, r, s)); - return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r, s); - } - - Oriented_side - power_test(const Weighted_point &p, const Weighted_point &q, - const Weighted_point &r, const Weighted_point &s, - const Weighted_point &t) const - { - return geom_traits().power_side_of_oriented_power_sphere_3_object()(p, q, r, s, t); - } - - bool in_conflict_3(const Weighted_point &p, const Cell_handle c) const - { - return side_of_power_sphere(c, p, true) == ON_BOUNDED_SIDE; - } - - bool in_conflict_2(const Weighted_point &p, const Cell_handle c, int i) const - { - return side_of_power_circle(c, i, p, true) == ON_BOUNDED_SIDE; - } - - bool in_conflict_1(const Weighted_point &p, const Cell_handle c) const - { - return side_of_power_segment(c, p, true) == ON_BOUNDED_SIDE; - } - - bool in_conflict_0(const Weighted_point &p, const Cell_handle c) const - { - return power_test(c->vertex(0)->point(), p) == ON_POSITIVE_SIDE; - } - - bool in_conflict(const Weighted_point &p, const Cell_handle c) const - { - switch (dimension()) { - case 0: return in_conflict_0(p, c); - case 1: return in_conflict_1(p, c); - case 2: return in_conflict_2(p, c, 3); - case 3: return in_conflict_3(p, c); - } - return true; - } - - class Conflict_tester_3 - { - const Weighted_point &p; - const Self *t; - - public: - - Conflict_tester_3(const Weighted_point &pt, const Self *tr) - : p(pt), t(tr) {} - - bool operator()(const Cell_handle c) const { - return t->in_conflict_3(p, c); - } - - bool test_initial_cell(const Cell_handle c) const { - return operator()(c); - } - Oriented_side compare_weight(const Weighted_point &wp1, - const Weighted_point &wp2) const + int dim = t->dimension(); + while(start != end) { - return t->power_test (wp1, wp2); - } - }; + std::copy((*start)->hidden_points_begin(), + (*start)->hidden_points_end(), + std::back_inserter(hidden_points)); - class Conflict_tester_2 - { - const Weighted_point &p; - const Self *t; - public: - - Conflict_tester_2(const Weighted_point &pt, const Self *tr) - : p(pt), t(tr) {} - - bool operator()(const Cell_handle c) const - { - return t->in_conflict_2(p, c, 3); - } - bool test_initial_cell(const Cell_handle c) const { - return operator()(c); - } - Oriented_side compare_weight(const Weighted_point &wp1, - const Weighted_point &wp2) const - { - return t->power_test (wp1, wp2); - } - }; - - class Conflict_tester_1 - { - const Weighted_point &p; - const Self *t; - - public: - - Conflict_tester_1(const Weighted_point &pt, const Self *tr) - : p(pt), t(tr) {} - - bool operator()(const Cell_handle c) const - { - return t->in_conflict_1(p, c); - } - bool test_initial_cell(const Cell_handle c) const { - return operator()(c); - } - Oriented_side compare_weight(const Weighted_point &wp1, - const Weighted_point &wp2) const - { - return t->power_test (wp1, wp2); - } - }; - - class Conflict_tester_0 - { - const Weighted_point &p; - const Self *t; - - public: - - Conflict_tester_0(const Weighted_point &pt, const Self *tr) - : p(pt), t(tr) {} - - bool operator()(const Cell_handle c) const - { - return t->in_conflict_0(p, c); - } - bool test_initial_cell(const Cell_handle c) const { - return operator()(c); - } - int compare_weight(const Weighted_point &wp1, - const Weighted_point &wp2) const - { - return t->power_test (wp1, wp2); - } - }; - - // Sequential version - // "dummy" is here to allow the specialization (see below) - // See http://groups.google.com/group/comp.lang.c++.moderated/browse_thread/thread/285ab1eec49e1cb6 - template - class Hidden_point_visitor - { - Self *t; - mutable std::vector vertices; - mutable std::vector hidden_points; - - public: - - Hidden_point_visitor(Self *tr) : t(tr) {} - - template - void process_cells_in_conflict(InputIterator start, InputIterator end) const - { - int dim = t->dimension(); - while (start != end) { - std::copy((*start)->hidden_points_begin(), - (*start)->hidden_points_end(), - std::back_inserter(hidden_points)); - - for (int i=0; i<=dim; i++) { - Vertex_handle v = (*start)->vertex(i); - if (v->cell() != Cell_handle()) { - vertices.push_back(v); - v->set_cell(Cell_handle()); - } + for(int i=0; i<=dim; i++) + { + Vertex_handle v = (*start)->vertex(i); + if(v->cell() != Cell_handle()) + { + vertices.push_back(v); + v->set_cell(Cell_handle()); } - start ++; } + start ++; } - void reinsert_vertices(Vertex_handle v) { - Cell_handle hc = v->cell(); - for (typename std::vector::iterator - vi = vertices.begin(); vi != vertices.end(); ++vi) { - if ((*vi)->cell() != Cell_handle()) continue; - hc = t->locate ((*vi)->point(), hc); - hide_point(hc, (*vi)->point()); - t->tds().delete_vertex(*vi); - } - vertices.clear(); - for (typename std::vector::iterator - hp = hidden_points.begin(); hp != hidden_points.end(); ++hp) { - hc = t->locate (*hp, hc); - hide_point (hc, *hp); - } - hidden_points.clear(); + } + + void reinsert_vertices(Vertex_handle v) + { + Cell_handle hc = v->cell(); + for(typename std::vector::iterator + vi = vertices.begin(); vi != vertices.end(); ++vi) + { + if((*vi)->cell() != Cell_handle()) + continue; + + hc = t->locate((*vi)->point(), hc); + hide_point(hc, (*vi)->point()); + t->tds().delete_vertex(*vi); } - Vertex_handle replace_vertex(Cell_handle c, int index, - const Weighted_point &p) { - Vertex_handle v = c->vertex(index); - hide_point(c, v->point()); - v->set_point(p); - return v; + + vertices.clear(); + for(typename std::vector::iterator + hp = hidden_points.begin(); hp != hidden_points.end(); ++hp) + { + hc = t->locate(*hp, hc); + hide_point (hc, *hp); } - void hide_point(Cell_handle c, const Weighted_point &p) { - c->hide_point(p); - } - }; + hidden_points.clear(); + } + + Vertex_handle replace_vertex(Cell_handle c, int index, const Weighted_point& p) + { + Vertex_handle v = c->vertex(index); + hide_point(c, v->point()); + v->set_point(p); + return v; + } + + void hide_point(Cell_handle c, const Weighted_point& p) + { + c->hide_point(p); + } + }; #ifdef CGAL_LINKED_WITH_TBB - // Parallel version specialization - template - class Hidden_point_visitor + // Parallel version specialization + template + class Hidden_point_visitor + { + typedef Hidden_point_visitor HPV; + + Self *t; + mutable tbb::enumerable_thread_specific > vertices; + mutable tbb::enumerable_thread_specific > hidden_points; + + public: + Hidden_point_visitor(Self *tr) : t(tr) {} + + template + void process_cells_in_conflict(InputIterator start, InputIterator end) const { - typedef Hidden_point_visitor HPV; - - Self *t; - mutable tbb::enumerable_thread_specific > vertices; - mutable tbb::enumerable_thread_specific > hidden_points; - - public: - - Hidden_point_visitor(Self *tr) : t(tr) {} - - template - void process_cells_in_conflict(InputIterator start, InputIterator end) const + int dim = t->dimension(); + while(start != end) { - int dim = t->dimension(); - while (start != end) { - std::copy((*start)->hidden_points_begin(), - (*start)->hidden_points_end(), - std::back_inserter(hidden_points.local())); + std::copy((*start)->hidden_points_begin(), + (*start)->hidden_points_end(), + std::back_inserter(hidden_points.local())); - for (int i=0; i<=dim; i++) { - Vertex_handle v = (*start)->vertex(i); - if (v->cell() != Cell_handle()) { - vertices.local().push_back(v); - v->set_cell(Cell_handle()); - } + for(int i=0; i<=dim; i++) + { + Vertex_handle v = (*start)->vertex(i); + if(v->cell() != Cell_handle()) + { + vertices.local().push_back(v); + v->set_cell(Cell_handle()); } - start ++; } + start ++; } - void reinsert_vertices(Vertex_handle v) { - Cell_handle hc = v->cell(); - for (typename std::vector::iterator - vi = vertices.local().begin(); vi != vertices.local().end(); ++vi) { - if ((*vi)->cell() != Cell_handle()) continue; - hc = t->locate ((*vi)->point(), hc); - hide_point(hc, (*vi)->point()); - t->tds().delete_vertex(*vi); - } - vertices.local().clear(); - for (typename std::vector::iterator - hp = hidden_points.local().begin(); hp != hidden_points.local().end(); ++hp) { - hc = t->locate (*hp, hc); - hide_point (hc, *hp); - } - hidden_points.local().clear(); + } + + void reinsert_vertices(Vertex_handle v) + { + Cell_handle hc = v->cell(); + for(typename std::vector::iterator + vi = vertices.local().begin(); vi != vertices.local().end(); ++vi) + { + if((*vi)->cell() != Cell_handle()) + continue; + + hc = t->locate((*vi)->point(), hc); + hide_point(hc, (*vi)->point()); + t->tds().delete_vertex(*vi); } - Vertex_handle replace_vertex(Cell_handle c, int index, - const Weighted_point &p) { - Vertex_handle v = c->vertex(index); - hide_point(c, v->point()); - v->set_point(p); - return v; + + vertices.local().clear(); + for(typename std::vector::iterator + hp = hidden_points.local().begin(); hp != hidden_points.local().end(); ++hp) + { + hc = t->locate(*hp, hc); + hide_point (hc, *hp); } - void hide_point(Cell_handle c, const Weighted_point &p) { - c->hide_point(p); - } - }; + hidden_points.local().clear(); + } + + Vertex_handle replace_vertex(Cell_handle c, int index, const Weighted_point& p) + { + Vertex_handle v = c->vertex(index); + hide_point(c, v->point()); + v->set_point(p); + return v; + } + + void hide_point(Cell_handle c, const Weighted_point& p) + { + c->hide_point(p); + } + }; // Functor for parallel insert(begin, end) function template @@ -1313,51 +1340,50 @@ namespace CGAL { typedef typename RT::Weighted_point Weighted_point; typedef typename RT::Vertex_handle Vertex_handle; - RT & m_rt; - const std::vector & m_points; - tbb::enumerable_thread_specific & m_tls_hint; + RT& m_rt; + const std::vector& m_points; + tbb::enumerable_thread_specific& m_tls_hint; public: // Constructor - Insert_point(RT & rt, - const std::vector & points, - tbb::enumerable_thread_specific & tls_hint) - : m_rt(rt), m_points(points), m_tls_hint(tls_hint) + Insert_point(RT& rt, + const std::vector& points, + tbb::enumerable_thread_specific& tls_hint) + : m_rt(rt), m_points(points), m_tls_hint(tls_hint) {} // Constructor - Insert_point(const Insert_point &ip) - : m_rt(ip.m_rt), m_points(ip.m_points), m_tls_hint(ip.m_tls_hint) + Insert_point(const Insert_point& ip) + : m_rt(ip.m_rt), m_points(ip.m_points), m_tls_hint(ip.m_tls_hint) {} // operator() - void operator()( const tbb::blocked_range& r ) const + void operator()(const tbb::blocked_range& r) const { #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING static Profile_branch_counter_3 bcounter( - "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); + "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); #endif - Vertex_handle &hint = m_tls_hint.local(); - for( size_t i_point = r.begin() ; i_point != r.end() ; ++i_point) + Vertex_handle& hint = m_tls_hint.local(); + for(size_t i_point = r.begin() ; i_point != r.end() ; ++i_point) { bool success = false; - const Weighted_point &p = m_points[i_point]; + const Weighted_point& p = m_points[i_point]; while(!success) { - if (m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) + if(m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) { bool could_lock_zone; Locate_type lt; int li, lj; - Cell_handle c = m_rt.locate (p, lt, li, lj, hint->cell(), - &could_lock_zone); + Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(), &could_lock_zone); Vertex_handle v; - if (could_lock_zone) - v = m_rt.insert (p, lt, c, li, lj, &could_lock_zone); + if(could_lock_zone) + v = m_rt.insert(p, lt, c, li, lj, &could_lock_zone); - if (could_lock_zone) + if(could_lock_zone) { hint = (v == Vertex_handle() ? c->vertex(0) : v); m_rt.unlock_all_elements(); @@ -1394,27 +1420,27 @@ namespace CGAL { typedef typename RT::Vertex_handle Vertex_handle; typedef typename RT::Triangulation_data_structure::Vertex::Info Info; - RT & m_rt; - const std::vector & m_points; - const std::vector & m_infos; - const std::vector & m_indices; - tbb::enumerable_thread_specific & m_tls_hint; + RT& m_rt; + const std::vector& m_points; + const std::vector& m_infos; + const std::vector& m_indices; + tbb::enumerable_thread_specific& m_tls_hint; public: // Constructor - Insert_point_with_info(RT & rt, - const std::vector & points, - const std::vector & infos, - const std::vector & indices, - tbb::enumerable_thread_specific & tls_hint) - : m_rt(rt), m_points(points), m_infos(infos), m_indices(indices), - m_tls_hint(tls_hint) + Insert_point_with_info(RT& rt, + const std::vector& points, + const std::vector& infos, + const std::vector& indices, + tbb::enumerable_thread_specific& tls_hint) + : m_rt(rt), m_points(points), m_infos(infos), m_indices(indices), + m_tls_hint(tls_hint) {} // Constructor Insert_point_with_info(const Insert_point_with_info &ip) - : m_rt(ip.m_rt), m_points(ip.m_points), m_infos(ip.m_infos), - m_indices(ip.m_indices), m_tls_hint(ip.m_tls_hint) + : m_rt(ip.m_rt), m_points(ip.m_points), m_infos(ip.m_infos), + m_indices(ip.m_indices), m_tls_hint(ip.m_tls_hint) {} // operator() @@ -1422,32 +1448,32 @@ namespace CGAL { { #ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING static Profile_branch_counter_3 bcounter( - "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); + "early withdrawals / late withdrawals / successes [Delaunay_tri_3::insert]"); #endif - Vertex_handle &hint = m_tls_hint.local(); - for (size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) + Vertex_handle& hint = m_tls_hint.local(); + for(size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) { bool success = false; std::ptrdiff_t i_point = m_indices[i_idx]; - const Weighted_point &p = m_points[i_point]; - while (!success) + const Weighted_point& p = m_points[i_point]; + while(!success) { - if (m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) + if(m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) { bool could_lock_zone; Locate_type lt; int li, lj; Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(), - &could_lock_zone); + &could_lock_zone); Vertex_handle v; - if (could_lock_zone) + if(could_lock_zone) v = m_rt.insert(p, lt, c, li, lj, &could_lock_zone); - if (could_lock_zone) + if(could_lock_zone) { - if (v == Vertex_handle()) + if(v == Vertex_handle()) { hint = c->vertex(0); } @@ -1490,593 +1516,608 @@ namespace CGAL { typedef typename RT::Weighted_point Weighted_point; typedef typename RT::Vertex_handle Vertex_handle; - RT & m_rt; - const std::vector & m_vertices; - tbb::concurrent_vector & m_vertices_to_remove_sequentially; + RT& m_rt; + const std::vector& m_vertices; + tbb::concurrent_vector& m_vertices_to_remove_sequentially; public: // Constructor - Remove_point(RT & rt, - const std::vector & vertices, - tbb::concurrent_vector & - vertices_to_remove_sequentially) - : m_rt(rt), m_vertices(vertices), - m_vertices_to_remove_sequentially(vertices_to_remove_sequentially) + Remove_point(RT& rt, + const std::vector& vertices, + tbb::concurrent_vector& vertices_to_remove_sequentially) + : m_rt(rt), m_vertices(vertices), + m_vertices_to_remove_sequentially(vertices_to_remove_sequentially) {} // Constructor - Remove_point(const Remove_point &rp) - : m_rt(rp.m_rt), m_vertices(rp.m_vertices), - m_vertices_to_remove_sequentially(rp.m_vertices_to_remove_sequentially) + Remove_point(const Remove_point& rp) + : m_rt(rp.m_rt), m_vertices(rp.m_vertices), + m_vertices_to_remove_sequentially(rp.m_vertices_to_remove_sequentially) {} // operator() - void operator()( const tbb::blocked_range& r ) const + void operator()(const tbb::blocked_range& r) const { - for( size_t i_vertex = r.begin() ; i_vertex != r.end() ; ++i_vertex) + for(size_t i_vertex = r.begin() ; i_vertex != r.end() ; ++i_vertex) { Vertex_handle v = m_vertices[i_vertex]; bool could_lock_zone, needs_to_be_done_sequentially; do { - needs_to_be_done_sequentially = - !m_rt.remove(v, &could_lock_zone); + needs_to_be_done_sequentially = !m_rt.remove(v, &could_lock_zone); m_rt.unlock_all_elements(); - } while (!could_lock_zone); + } + while(!could_lock_zone); - if (needs_to_be_done_sequentially) + if(needs_to_be_done_sequentially) m_vertices_to_remove_sequentially.push_back(v); } } }; #endif // CGAL_LINKED_WITH_TBB - Hidden_point_visitor &get_hidden_point_visitor() - { - return hidden_point_visitor; - } - - template < class RegularTriangulation_3 > - class Vertex_remover; - - template < class RegularTriangulation_3 > - class Vertex_inserter; - - Hidden_point_visitor hidden_point_visitor; - }; - - - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - nearest_power_vertex_in_cell(const Bare_point& p, - Cell_handle c) const - // Returns the finite vertex of the cell c with smaller - // power distance to p. + Hidden_point_visitor& get_hidden_point_visitor() { - CGAL_triangulation_precondition(dimension() >= 1); - Vertex_handle nearest = nearest_power_vertex(p, - c->vertex(0), - c->vertex(1)); - if (dimension() >= 2) { - nearest = nearest_power_vertex(p, nearest, c->vertex(2)); - if (dimension() == 3) - nearest = nearest_power_vertex(p, nearest, c->vertex(3)); - } - return nearest; + return hidden_point_visitor; } + template < class RegularTriangulation_3 > + class Vertex_remover; - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - nearest_power_vertex(const Bare_point& p, Cell_handle start) const + template < class RegularTriangulation_3 > + class Vertex_inserter; + + Hidden_point_visitor hidden_point_visitor; +}; + + +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +nearest_power_vertex_in_cell(const Bare_point& p, Cell_handle c) const +// Returns the finite vertex of the cell c with smaller +// power distance to p. +{ + CGAL_triangulation_precondition(dimension() >= 1); + Vertex_handle nearest = nearest_power_vertex(p, c->vertex(0), c->vertex(1)); + if(dimension() >= 2) { - if (number_of_vertices() == 0) - return Vertex_handle(); + nearest = nearest_power_vertex(p, nearest, c->vertex(2)); + if(dimension() == 3) + nearest = nearest_power_vertex(p, nearest, c->vertex(3)); + } + return nearest; +} - // Use a brute-force algorithm if dimension < 3. - if (dimension() < 3) { - Finite_vertices_iterator vit = finite_vertices_begin(); - Vertex_handle res = vit; - ++vit; - for (Finite_vertices_iterator end = finite_vertices_end(); vit != end; ++vit) - res = nearest_power_vertex(p, res, vit); - return res; - } - Locate_type lt; - int li, lj; +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +nearest_power_vertex(const Bare_point& p, Cell_handle start) const +{ + if(number_of_vertices() == 0) + return Vertex_handle(); - typename Gt::Construct_weighted_point_3 p2wp = - geom_traits().construct_weighted_point_3_object(); - Cell_handle c = locate(p2wp(p), lt, li, lj, start); + // Use a brute-force algorithm if dimension < 3. + if(dimension() < 3) + { + Finite_vertices_iterator vit = finite_vertices_begin(); + Vertex_handle res = vit; + ++vit; + for(Finite_vertices_iterator end = finite_vertices_end(); vit != end; ++vit) + res = nearest_power_vertex(p, res, vit); - // - start with the closest vertex from the located cell. - // - repeatedly take the nearest of its incident vertices if any - // - if not, we're done. - Vertex_handle nearest = nearest_power_vertex_in_cell(p, c); - std::vector vs; - vs.reserve(32); - while (true) { - Vertex_handle tmp = nearest; - adjacent_vertices(nearest, std::back_inserter(vs)); - for (typename std::vector::const_iterator - vsit = vs.begin(); vsit != vs.end(); ++vsit) - tmp = nearest_power_vertex(p, tmp, *vsit); - if (tmp == nearest) - break; - vs.clear(); - nearest = tmp; - } - return nearest; + return res; } - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Bare_point - Regular_triangulation_3:: - dual(Cell_handle c) const - { - CGAL_triangulation_precondition(dimension()==3); - CGAL_triangulation_precondition( ! is_infinite(c) ); + Locate_type lt; + int li, lj; - return c->weighted_circumcenter(geom_traits()); + typename Gt::Construct_weighted_point_3 p2wp = geom_traits().construct_weighted_point_3_object(); + Cell_handle c = locate(p2wp(p), lt, li, lj, start); + + // - start with the closest vertex from the located cell. + // - repeatedly take the nearest of its incident vertices if any + // - if not, we're done. + Vertex_handle nearest = nearest_power_vertex_in_cell(p, c); + std::vector vs; + vs.reserve(32); + while(true) + { + Vertex_handle tmp = nearest; + adjacent_vertices(nearest, std::back_inserter(vs)); + for(typename std::vector::const_iterator + vsit = vs.begin(); vsit != vs.end(); ++vsit) + tmp = nearest_power_vertex(p, tmp, *vsit); + + if(tmp == nearest) + break; + + vs.clear(); + nearest = tmp; + } + return nearest; +} + +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Bare_point +Regular_triangulation_3:: +dual(Cell_handle c) const +{ + CGAL_triangulation_precondition(dimension()==3); + CGAL_triangulation_precondition(! is_infinite(c)); + + return c->weighted_circumcenter(geom_traits()); +} + +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const +{ + Cell_handle n = c->neighbor(i); + CGAL_assertion(! is_infinite(c) && ! is_infinite(n)); + + p = dual(c); + q = dual(n); +} + +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const +{ + return dual_segment(facet.first, facet.second, p, q); +} + +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_ray(Cell_handle c, int i, Ray& ray) const +{ + Cell_handle n = c->neighbor(i); + CGAL_triangulation_precondition((!is_infinite(c) != !is_infinite(n))); // xor + // either n or c is infinite + int in; + if(is_infinite(c)) + { + in = n->index(c); + } + else + { + n = c; + in = i; } - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const - { - Cell_handle n = c->neighbor(i); - CGAL_assertion( ! is_infinite(c) && ! is_infinite(n) ); + // n now denotes a finite cell, either c or c->neighbor(i) + int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; + if((in&1) == 1) + std::swap(ind[0], ind[1]); - p = dual(c); - q = dual(n); - } + const Weighted_point& p = n->vertex(ind[0])->point(); + const Weighted_point& q = n->vertex(ind[1])->point(); + const Weighted_point& r = n->vertex(ind[2])->point(); + const Bare_point& bp = construct_point(p); + const Bare_point& bq = construct_point(q); + const Bare_point& br = construct_point(r); - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const - { - return dual_segment(facet.first, facet.second, p, q); - } + Line l = construct_perpendicular_line(construct_plane(bp, bq, br), + construct_weighted_circumcenter(p,q,r)); - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_ray(Cell_handle c, int i, Ray& ray) const - { - Cell_handle n = c->neighbor(i); - CGAL_triangulation_precondition( (!is_infinite(c) != !is_infinite(n)) ); // xor - // either n or c is infinite - int in; - if ( is_infinite(c) ) - in = n->index(c); - else { - n = c; - in = i; - } - // n now denotes a finite cell, either c or c->neighbor(i) - int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; - if ( (in&1) == 1 ) - std::swap(ind[0], ind[1]); + ray = construct_ray(dual(n), l); +} - const Weighted_point& p = n->vertex(ind[0])->point(); - const Weighted_point& q = n->vertex(ind[1])->point(); - const Weighted_point& r = n->vertex(ind[2])->point(); - const Bare_point& bp = construct_point(p); - const Bare_point& bq = construct_point(q); - const Bare_point& br = construct_point(r); +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_ray(const Facet& facet, Ray& ray) const +{ + return dual_ray(facet.first, facet.second, ray); +} - Line l = construct_perpendicular_line(construct_plane(bp, bq, br), - construct_weighted_circumcenter(p,q,r)); +// Exact versions of dual_segment() and dual_ray() for Mesh_3. +// These functions are really dirty: they assume that the point type is nice enough +// such that EPECK can manipulate it (e.g. convert it to EPECK::Point_3) AND +// that the result of these manipulations will make sense. +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const +{ + typedef typename Gt::Kernel K; + typedef Exact_predicates_exact_constructions_kernel EK; + typedef Cartesian_converter To_exact; + typedef Cartesian_converter Back_from_exact; - ray = construct_ray(dual(n), l); - } + typedef EK Exact_Rt; - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_ray(const Facet& facet, Ray& ray) const - { - return dual_ray(facet.first, facet.second, ray); - } - - // Exact versions of dual_segment() and dual_ray() for Mesh_3. - // These functions are really dirty: they assume that the point type is nice enough - // such that EPECK can manipulate it (e.g. convert it to EPECK::Point_3) AND - // that the result of these manipulations will make sense. - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const - { - typedef typename Gt::Kernel K; - typedef Exact_predicates_exact_constructions_kernel EK; - typedef Cartesian_converter To_exact; - typedef Cartesian_converter Back_from_exact; - - typedef EK Exact_Rt; - - To_exact to_exact; - Back_from_exact back_from_exact; - Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = + To_exact to_exact; + Back_from_exact back_from_exact; + Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = Exact_Rt().construct_weighted_circumcenter_3_object(); - Cell_handle c = facet.first; - int i = facet.second; - Cell_handle n = c->neighbor(i); + Cell_handle c = facet.first; + int i = facet.second; + Cell_handle n = c->neighbor(i); - const typename Exact_Rt::Weighted_point_3& cp = to_exact(c->vertex(0)->point()); - const typename Exact_Rt::Weighted_point_3& cq = to_exact(c->vertex(1)->point()); - const typename Exact_Rt::Weighted_point_3& cr = to_exact(c->vertex(2)->point()); - const typename Exact_Rt::Weighted_point_3& cs = to_exact(c->vertex(3)->point()); + const typename Exact_Rt::Weighted_point_3& cp = to_exact(c->vertex(0)->point()); + const typename Exact_Rt::Weighted_point_3& cq = to_exact(c->vertex(1)->point()); + const typename Exact_Rt::Weighted_point_3& cr = to_exact(c->vertex(2)->point()); + const typename Exact_Rt::Weighted_point_3& cs = to_exact(c->vertex(3)->point()); - const typename Exact_Rt::Weighted_point_3& np = to_exact(n->vertex(0)->point()); - const typename Exact_Rt::Weighted_point_3& nq = to_exact(n->vertex(1)->point()); - const typename Exact_Rt::Weighted_point_3& nr = to_exact(n->vertex(2)->point()); - const typename Exact_Rt::Weighted_point_3& ns = to_exact(n->vertex(3)->point()); + const typename Exact_Rt::Weighted_point_3& np = to_exact(n->vertex(0)->point()); + const typename Exact_Rt::Weighted_point_3& nq = to_exact(n->vertex(1)->point()); + const typename Exact_Rt::Weighted_point_3& nr = to_exact(n->vertex(2)->point()); + const typename Exact_Rt::Weighted_point_3& ns = to_exact(n->vertex(3)->point()); - p = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs)); - q = back_from_exact(exact_weighted_circumcenter(np, nq, nr, ns)); + p = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs)); + q = back_from_exact(exact_weighted_circumcenter(np, nq, nr, ns)); +} + +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +dual_ray_exact(const Facet& facet, Ray& ray) const +{ + Cell_handle c = facet.first; + int i = facet.second; + Cell_handle n = c->neighbor(i); + CGAL_triangulation_precondition(!is_infinite(c) != !is_infinite(n)); // xor + // either n or c is infinite + int in; + if(is_infinite(c)) + { + in = n->index(c); + } + else + { + n = c; + in = i; } - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - dual_ray_exact(const Facet& facet, Ray& ray) const - { - Cell_handle c = facet.first; - int i = facet.second; - Cell_handle n = c->neighbor(i); - CGAL_triangulation_precondition( !is_infinite(c) != !is_infinite(n) ); // xor - // either n or c is infinite - int in; - if ( is_infinite(c) ) - in = n->index(c); - else { - n = c; - in = i; - } - // n now denotes a finite cell, either c or c->neighbor(i) - int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; - if ( (in&1) == 1 ) - std::swap(ind[0], ind[1]); + // n now denotes a finite cell, either c or c->neighbor(i) + int ind[3] = {(in+1)&3,(in+2)&3,(in+3)&3}; + if((in&1) == 1) + std::swap(ind[0], ind[1]); - // exact part - typedef typename Gt::Kernel K; - typedef Exact_predicates_exact_constructions_kernel EK; - typedef Cartesian_converter To_exact; - typedef Cartesian_converter Back_from_exact; + // exact part + typedef typename Gt::Kernel K; + typedef Exact_predicates_exact_constructions_kernel EK; + typedef Cartesian_converter To_exact; + typedef Cartesian_converter Back_from_exact; - typedef EK Exact_Rt; + typedef EK Exact_Rt; - To_exact to_exact; - Back_from_exact back_from_exact; + To_exact to_exact; + Back_from_exact back_from_exact; - Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = + Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = Exact_Rt().construct_weighted_circumcenter_3_object(); - Exact_Rt::Construct_perpendicular_line_3 exact_perpendicular_line = + Exact_Rt::Construct_perpendicular_line_3 exact_perpendicular_line = Exact_Rt().construct_perpendicular_line_3_object(); - Exact_Rt::Construct_plane_3 exact_plane_3 = Exact_Rt().construct_plane_3_object(); - Exact_Rt::Construct_ray_3 exact_ray_3 = Exact_Rt().construct_ray_3_object(); - Exact_Rt::Construct_point_3 exact_point_3 = Exact_Rt().construct_point_3_object(); + Exact_Rt::Construct_plane_3 exact_plane_3 = Exact_Rt().construct_plane_3_object(); + Exact_Rt::Construct_ray_3 exact_ray_3 = Exact_Rt().construct_ray_3_object(); + Exact_Rt::Construct_point_3 exact_point_3 = Exact_Rt().construct_point_3_object(); - const typename Exact_Rt::Weighted_point_3& p = to_exact(n->vertex(ind[0])->point()); - const typename Exact_Rt::Weighted_point_3& q = to_exact(n->vertex(ind[1])->point()); - const typename Exact_Rt::Weighted_point_3& r = to_exact(n->vertex(ind[2])->point()); - const typename Exact_Rt::Weighted_point_3& s = to_exact(n->vertex(in)->point()); + const typename Exact_Rt::Weighted_point_3& p = to_exact(n->vertex(ind[0])->point()); + const typename Exact_Rt::Weighted_point_3& q = to_exact(n->vertex(ind[1])->point()); + const typename Exact_Rt::Weighted_point_3& r = to_exact(n->vertex(ind[2])->point()); + const typename Exact_Rt::Weighted_point_3& s = to_exact(n->vertex(in)->point()); - const typename Exact_Rt::Point_3& bp = exact_point_3(p); - const typename Exact_Rt::Point_3& bq = exact_point_3(q); - const typename Exact_Rt::Point_3& br = exact_point_3(r); + const typename Exact_Rt::Point_3& bp = exact_point_3(p); + const typename Exact_Rt::Point_3& bq = exact_point_3(q); + const typename Exact_Rt::Point_3& br = exact_point_3(r); - typename Exact_Rt::Line_3 l = exact_perpendicular_line( - exact_plane_3(bp, bq, br), - exact_weighted_circumcenter(p, q, r)); + typename Exact_Rt::Line_3 l = exact_perpendicular_line( + exact_plane_3(bp, bq, br), + exact_weighted_circumcenter(p, q, r)); - ray = back_from_exact( - exact_ray_3( - exact_weighted_circumcenter(p, q, r, s), l)); + ray = back_from_exact(exact_ray_3(exact_weighted_circumcenter(p, q, r, s), l)); +} + +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Object +Regular_triangulation_3:: +dual(Cell_handle c, int i) const +{ + CGAL_triangulation_precondition(dimension()>=2); + CGAL_triangulation_precondition(! is_infinite(c,i)); + + if(dimension() == 2) + { + CGAL_triangulation_precondition(i == 3); + return construct_object(construct_weighted_circumcenter(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point())); } - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Object - Regular_triangulation_3:: - dual(Cell_handle c, int i) const + // dimension() == 3 + Cell_handle n = c->neighbor(i); + if(! is_infinite(c) && ! is_infinite(n)) { - CGAL_triangulation_precondition(dimension()>=2); - CGAL_triangulation_precondition( ! is_infinite(c,i) ); - - if ( dimension() == 2 ) { - CGAL_triangulation_precondition( i == 3 ); - return construct_object( - construct_weighted_circumcenter(c->vertex(0)->point(), - c->vertex(1)->point(), - c->vertex(2)->point()) ); - } - - // dimension() == 3 - Cell_handle n = c->neighbor(i); - if ( ! is_infinite(c) && ! is_infinite(n) ){ - // dual is a segment - Bare_point bp = dual(c); - Bare_point np = dual(n); - return construct_object(construct_segment(bp, np)); - } - - // either n or c is infinite, dual is a ray - Ray r; - dual_ray(c, i, r); - return construct_object(r); + // dual is a segment + Bare_point bp = dual(c); + Bare_point np = dual(n); + return construct_object(construct_segment(bp, np)); } - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Object - Regular_triangulation_3:: - dual(const Facet& facet) const - { - return dual(facet.first, facet.second); - } + // either n or c is infinite, dual is a ray + Ray r; + dual_ray(c, i, r); + return construct_object(r); +} - template < class Gt, class Tds, class Lds > - template < class Stream> - Stream& - Regular_triangulation_3:: - draw_dual(Stream & os) const +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Object +Regular_triangulation_3:: +dual(const Facet& facet) const +{ + return dual(facet.first, facet.second); +} + +template < class Gt, class Tds, class Lds > +template < class Stream> +Stream& +Regular_triangulation_3:: +draw_dual(Stream& os) const +{ + for(Finite_facets_iterator fit = finite_facets_begin(), end = finite_facets_end(); + fit != end; ++fit) { - for (Finite_facets_iterator fit = finite_facets_begin(), - end = finite_facets_end(); - fit != end; ++fit) { - Object o = dual(*fit); - if (const Segment *s = object_cast(&o)) os << *s; - else if (const Ray *r = object_cast(&o)) os << *r; - else if (const Bare_point *p = object_cast(&o)) os << *p; - } + Object o = dual(*fit); + if(const Segment* s = object_cast(&o)) + os << *s; + else if(const Ray* r = object_cast(&o)) + os << *r; + else if(const Bare_point* p = object_cast(&o)) + os << *p; + } + return os; +} + +template < class Gt, class Tds, class Lds > +Oriented_side +Regular_triangulation_3:: +side_of_oriented_power_sphere(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p3, + const Weighted_point& p, bool perturb) const +{ + CGAL_triangulation_precondition(orientation(p0, p1, p2, p3) == POSITIVE); + + using namespace boost; + + Oriented_side os = power_test(p0, p1, p2, p3, p); + + if(os != ON_ORIENTED_BOUNDARY || !perturb) return os; + + // We are now in a degenerate case => we do a symbolic perturbation. + + // We sort the points lexicographically. + const Weighted_point * points[5] = {&p0, &p1, &p2, &p3, &p}; + + std::sort(points, points + 5, typename Tr_Base::Perturbation_order(this)); + + // We successively look whether the leading monomial, then 2nd monomial + // of the determinant has non null coefficient. + for(int i=4; i>1; --i) + { + if(points[i] == &p) + return ON_NEGATIVE_SIDE; // since p0 p1 p2 p3 are non coplanar and positively oriented + Orientation o; + if(points[i] == &p3 && (o = orientation(p0,p1,p2,p)) != COPLANAR) + return o; + if(points[i] == &p2 && (o = orientation(p0,p1,p,p3)) != COPLANAR) + return o; + if(points[i] == &p1 && (o = orientation(p0,p,p2,p3)) != COPLANAR) + return o; + if(points[i] == &p0 && (o = orientation(p,p1,p2,p3)) != COPLANAR) + return o; } - template < class Gt, class Tds, class Lds > - Oriented_side - Regular_triangulation_3:: - side_of_oriented_power_sphere(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p3, - const Weighted_point &p, bool perturb) const + CGAL_triangulation_assertion(false); + return ON_NEGATIVE_SIDE; +} + +template < class Gt, class Tds, class Lds > +Bounded_side +Regular_triangulation_3:: +side_of_power_sphere(Cell_handle c, const Weighted_point& p, bool perturb) const +{ + CGAL_triangulation_precondition(dimension() == 3); + int i3; + if(! c->has_vertex(infinite_vertex(), i3)) { - CGAL_triangulation_precondition( orientation(p0, p1, p2, p3) == POSITIVE ); - - using namespace boost; - - Oriented_side os = power_test(p0, p1, p2, p3, p); - - if (os != ON_ORIENTED_BOUNDARY || !perturb) - return os; - - // We are now in a degenerate case => we do a symbolic perturbation. - - // We sort the points lexicographically. - const Weighted_point * points[5] = {&p0, &p1, &p2, &p3, &p}; - - std::sort(points, points + 5, typename Tr_Base::Perturbation_order(this)); - - // We successively look whether the leading monomial, then 2nd monomial - // of the determinant has non null coefficient. - for (int i=4; i>1; --i) { - if (points[i] == &p) - return ON_NEGATIVE_SIDE; // since p0 p1 p2 p3 are non coplanar - // and positively oriented - Orientation o; - if (points[i] == &p3 && (o = orientation(p0,p1,p2,p)) != COPLANAR ) - return o; - if (points[i] == &p2 && (o = orientation(p0,p1,p,p3)) != COPLANAR ) - return o; - if (points[i] == &p1 && (o = orientation(p0,p,p2,p3)) != COPLANAR ) - return o; - if (points[i] == &p0 && (o = orientation(p,p1,p2,p3)) != COPLANAR ) - return o; - } - - CGAL_triangulation_assertion(false); - return ON_NEGATIVE_SIDE; + return Bounded_side(side_of_oriented_power_sphere(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + c->vertex(3)->point(), + p, perturb)); } - - template < class Gt, class Tds, class Lds > - Bounded_side - Regular_triangulation_3:: - side_of_power_sphere(Cell_handle c, const Weighted_point &p, - bool perturb) const + // else infinite cell : + int i0,i1,i2; + if((i3%2) == 1) { - CGAL_triangulation_precondition( dimension() == 3 ); - int i3; - if ( ! c->has_vertex( infinite_vertex(), i3 ) ) { - return Bounded_side( side_of_oriented_power_sphere(c->vertex(0)->point(), - c->vertex(1)->point(), - c->vertex(2)->point(), - c->vertex(3)->point(), - p, perturb) ); - } - // else infinite cell : - int i0,i1,i2; - if ( (i3%2) == 1 ) { - i0 = (i3+1)&3; - i1 = (i3+2)&3; - i2 = (i3+3)&3; - } - else { - i0 = (i3+2)&3; - i1 = (i3+1)&3; - i2 = (i3+3)&3; - } - - // general case - Orientation o = orientation(c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), p); - if (o != ZERO) - return Bounded_side(o); - - // else p coplanar with i0,i1,i2 - return side_of_bounded_power_circle(c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), - p, perturb); + i0 = (i3+1)&3; + i1 = (i3+2)&3; + i2 = (i3+3)&3; + } + else + { + i0 = (i3+2)&3; + i1 = (i3+1)&3; + i2 = (i3+3)&3; } + // general case + Orientation o = orientation(c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), p); + if(o != ZERO) + return Bounded_side(o); - template < class Gt, class Tds, class Lds > - Bounded_side - Regular_triangulation_3:: - side_of_bounded_power_circle(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p, bool perturb) const + // else p coplanar with i0,i1,i2 + return side_of_bounded_power_circle(c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), + p, perturb); +} + +template < class Gt, class Tds, class Lds > +Bounded_side +Regular_triangulation_3:: +side_of_bounded_power_circle(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p, bool perturb) const +{ + CGAL_triangulation_precondition(coplanar_orientation(p0, p1, p2) != 0); + if(coplanar_orientation(p0, p1, p2) == POSITIVE) + return Bounded_side (side_of_oriented_power_circle(p0, p1, p2, p, perturb)); + + // Wrong because the low level power test already does a coplanar orientation test. + // return Bounded_side (- side_of_oriented_power_circle (p0, p2, p1, p, perturb)); + + return Bounded_side (side_of_oriented_power_circle(p0, p2, p1, p, perturb)); +} + + +template < class Gt, class Tds, class Lds > +Oriented_side +Regular_triangulation_3:: +side_of_oriented_power_circle(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p2, + const Weighted_point& p, bool perturb) const +{ + CGAL_triangulation_precondition(coplanar_orientation(p0, p1, p2) == POSITIVE); + + using namespace boost; + + Oriented_side os = power_test(p0, p1, p2, p); + + if(os != ON_ORIENTED_BOUNDARY || !perturb) + return os; + + // We are now in a degenerate case => we do a symbolic perturbation. + + // We sort the points lexicographically. + const Weighted_point * points[4] = {&p0, &p1, &p2, &p}; + + std::sort(points, points + 4, typename Tr_Base::Perturbation_order(this)); + + // We successively look whether the leading monomial, then 2nd monomial + // of the determinant has non null coefficient. + // 2 iterations are enough (cf paper) + for(int i=3; i>1; --i) { - CGAL_triangulation_precondition(coplanar_orientation(p0, p1, p2) != 0); - if (coplanar_orientation(p0, p1, p2) == POSITIVE) - return Bounded_side (side_of_oriented_power_circle(p0, p1, p2, p, perturb)); - // Wrong because the low level power test already does a coplanar orientation - // test. - // return Bounded_side (- side_of_oriented_power_circle (p0, p2, p1, p, - // perturb)); - return Bounded_side (side_of_oriented_power_circle(p0, p2, p1, p, perturb)); + if(points[i] == &p) + return ON_NEGATIVE_SIDE; // since p0 p1 p2 are non collinear and positively oriented + Orientation o; + if(points[i] == &p2 && (o = coplanar_orientation(p0,p1,p)) != COPLANAR) + return o; + if(points[i] == &p1 && (o = coplanar_orientation(p0,p,p2)) != COPLANAR) + return o; + if(points[i] == &p0 && (o = coplanar_orientation(p,p1,p2)) != COPLANAR) + return o; } + CGAL_triangulation_assertion(false); + return ON_NEGATIVE_SIDE; +} - template < class Gt, class Tds, class Lds > - Oriented_side - Regular_triangulation_3:: - side_of_oriented_power_circle(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p2, - const Weighted_point &p, bool perturb) const +template < class Gt, class Tds, class Lds > +Bounded_side +Regular_triangulation_3:: +side_of_power_circle(Cell_handle c, int i, const Weighted_point& p, + bool perturb) const +{ + CGAL_triangulation_precondition(dimension() >= 2); + int i3 = 5; + if(dimension() == 2) { - CGAL_triangulation_precondition( coplanar_orientation(p0, p1, p2) == POSITIVE ); - - using namespace boost; - - Oriented_side os = power_test(p0, p1, p2, p); - - if (os != ON_ORIENTED_BOUNDARY || !perturb) - return os; - - // We are now in a degenerate case => we do a symbolic perturbation. - - // We sort the points lexicographically. - const Weighted_point * points[4] = {&p0, &p1, &p2, &p}; - - std::sort(points, points + 4, typename Tr_Base::Perturbation_order(this)); - - // We successively look whether the leading monomial, then 2nd monomial - // of the determinant has non null coefficient. - // 2 iterations are enough (cf paper) - for (int i=3; i>1; --i) { - if (points[i] == &p) - return ON_NEGATIVE_SIDE; // since p0 p1 p2 are non collinear - // and positively oriented - Orientation o; - if (points[i] == &p2 && (o = coplanar_orientation(p0,p1,p)) != COPLANAR ) - return o; - if (points[i] == &p1 && (o = coplanar_orientation(p0,p,p2)) != COPLANAR ) - return o; - if (points[i] == &p0 && (o = coplanar_orientation(p,p1,p2)) != COPLANAR ) - return o; - } - - CGAL_triangulation_assertion(false); - return ON_NEGATIVE_SIDE; - } - - - template < class Gt, class Tds, class Lds > - Bounded_side - Regular_triangulation_3:: - side_of_power_circle(Cell_handle c, int i, const Weighted_point &p, - bool perturb) const - { - CGAL_triangulation_precondition( dimension() >= 2 ); - int i3 = 5; - if ( dimension() == 2 ) { - CGAL_triangulation_precondition( i == 3 ); - // the triangulation is supposed to be valid, ie the facet - // with vertices 0 1 2 in this order is positively oriented - if ( ! c->has_vertex( infinite_vertex(), i3 ) ) - return Bounded_side( side_of_oriented_power_circle(c->vertex(0)->point(), - c->vertex(1)->point(), - c->vertex(2)->point(), - p, perturb) ); - // else infinite facet - // v1, v2 finite vertices of the facet such that v1,v2,infinite - // is positively oriented - Vertex_handle v1 = c->vertex( ccw(i3) ), - v2 = c->vertex( cw(i3) ); - CGAL_triangulation_assertion(coplanar_orientation(v1->point(), v2->point(), - mirror_vertex(c, i3)->point()) == NEGATIVE); - Orientation o = coplanar_orientation(v1->point(), v2->point(), p); - if ( o != ZERO ) - return Bounded_side( o ); - // case when p collinear with v1v2 - return side_of_bounded_power_segment(v1->point(), - v2->point(), - p, perturb); - }// dim 2 - - // else dimension == 3 - CGAL_triangulation_precondition( (i >= 0) && (i < 4) ); - if ( ( ! c->has_vertex(infinite_vertex(),i3) ) || ( i3 != i ) ) { - // finite facet - // initialization of i0 i1 i2, vertices of the facet positively - // oriented (if the triangulation is valid) - int i0 = (i>0) ? 0 : 1; - int i1 = (i>1) ? 1 : 2; - int i2 = (i>2) ? 2 : 3; - CGAL_triangulation_precondition(this->coplanar(c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), p)); - return side_of_bounded_power_circle(c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), - p, perturb); - } - //else infinite facet + CGAL_triangulation_precondition(i == 3); + // the triangulation is supposed to be valid, ie the facet + // with vertices 0 1 2 in this order is positively oriented + if(! c->has_vertex(infinite_vertex(), i3)) + return Bounded_side(side_of_oriented_power_circle(c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + p, perturb)); + // else infinite facet // v1, v2 finite vertices of the facet such that v1,v2,infinite // is positively oriented - Vertex_handle v1 = c->vertex( next_around_edge(i3,i) ), - v2 = c->vertex( next_around_edge(i,i3) ); - Orientation o = (Orientation) - (coplanar_orientation( v1->point(), v2->point(), - c->vertex(i)->point()) * - coplanar_orientation( v1->point(), v2->point(), p)); - // then the code is duplicated from 2d case - if ( o != ZERO ) - return Bounded_side( -o ); - // because p is in f iff - // it is not on the same side of v1v2 as c->vertex(i) - // case when p collinear with v1v2 : + Vertex_handle v1 = c->vertex(ccw(i3)), + v2 = c->vertex(cw(i3)); + CGAL_triangulation_assertion( + coplanar_orientation(v1->point(), v2->point(), mirror_vertex(c, i3)->point()) == NEGATIVE); + + Orientation o = coplanar_orientation(v1->point(), v2->point(), p); + if(o != ZERO) + return Bounded_side(o); + + // case when p collinear with v1v2 return side_of_bounded_power_segment(v1->point(), - v2->point(), - p, perturb); - } + v2->point(), + p, perturb); + } // dim 2 - template < class Gt, class Tds, class Lds > - Bounded_side - Regular_triangulation_3:: - side_of_bounded_power_segment(const Weighted_point &p0, - const Weighted_point &p1, - const Weighted_point &p, bool perturb) const + // else dimension == 3 + CGAL_triangulation_precondition((i >= 0) && (i < 4)); + if((! c->has_vertex(infinite_vertex(),i3)) || (i3 != i)) { - Oriented_side os = power_test(p0, p1, p); + // finite facet + // initialization of i0 i1 i2, vertices of the facet positively + // oriented (if the triangulation is valid) + int i0 = (i>0) ? 0 : 1; + int i1 = (i>1) ? 1 : 2; + int i2 = (i>2) ? 2 : 3; + CGAL_triangulation_precondition(this->coplanar(c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), p)); + return side_of_bounded_power_circle(c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), + p, perturb); + } + //else infinite facet + // v1, v2 finite vertices of the facet such that v1,v2,infinite + // is positively oriented + Vertex_handle v1 = c->vertex(next_around_edge(i3,i)), + v2 = c->vertex(next_around_edge(i,i3)); + Orientation o = (Orientation) + (coplanar_orientation(v1->point(), v2->point(), + c->vertex(i)->point()) * + coplanar_orientation(v1->point(), v2->point(), p)); + // then the code is duplicated from 2d case + if(o != ZERO) + return Bounded_side(-o); + // because p is in f iff + // it is not on the same side of v1v2 as c->vertex(i) + // case when p collinear with v1v2 : + return side_of_bounded_power_segment(v1->point(), + v2->point(), + p, perturb); +} - if (os != ON_ORIENTED_BOUNDARY || !perturb) - return Bounded_side(os); +template < class Gt, class Tds, class Lds > +Bounded_side +Regular_triangulation_3:: +side_of_bounded_power_segment(const Weighted_point& p0, + const Weighted_point& p1, + const Weighted_point& p, bool perturb) const +{ + Oriented_side os = power_test(p0, p1, p); - // We are now in a degenerate case => we do a symbolic perturbation. + if(os != ON_ORIENTED_BOUNDARY || !perturb) + return Bounded_side(os); - switch (this->collinear_position(p0, p, p1)) { + // We are now in a degenerate case => we do a symbolic perturbation. + switch (this->collinear_position(p0, p, p1)) + { case Tr_Base::BEFORE: case Tr_Base::AFTER: return ON_UNBOUNDED_SIDE; @@ -2084,543 +2125,568 @@ namespace CGAL { return ON_BOUNDED_SIDE; default: ; - } - - CGAL_triangulation_assertion(false); - return ON_UNBOUNDED_SIDE; } - template < class Gt, class Tds, class Lds > - Bounded_side - Regular_triangulation_3:: - side_of_power_segment(Cell_handle c, const Weighted_point &p, - bool perturb) const - { - CGAL_triangulation_precondition( dimension() == 1 ); - if ( ! is_infinite(c,0,1) ) - return side_of_bounded_power_segment(c->vertex(0)->point(), - c->vertex(1)->point(), - p, perturb); - Locate_type lt; int i; - Bounded_side soe = side_of_edge( p, c, lt, i ); - if (soe != ON_BOUNDARY) - return soe; - // Either we compare weights, or we use the finite neighboring edge - Cell_handle finite_neighbor = c->neighbor(c->index(infinite_vertex())); - CGAL_triangulation_assertion(!is_infinite(finite_neighbor,0,1)); - return side_of_bounded_power_segment(finite_neighbor->vertex(0)->point(), - finite_neighbor->vertex(1)->point(), - p, perturb); - } + CGAL_triangulation_assertion(false); + return ON_UNBOUNDED_SIDE; +} - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - greater_or_equal_power_distance(const Bare_point& p, - const Weighted_point& q, - const Weighted_point& r) const - { - return ! less_power_distance(p, q, r); - } +template < class Gt, class Tds, class Lds > +Bounded_side +Regular_triangulation_3:: +side_of_power_segment(Cell_handle c, const Weighted_point& p, bool perturb) const +{ + CGAL_triangulation_precondition(dimension() == 1); + if(! is_infinite(c,0,1)) + return side_of_bounded_power_segment(c->vertex(0)->point(), + c->vertex(1)->point(), + p, perturb); - // @fixme, below can't be defined there because RegularTriangulationTraits_3 - // does not require the following construction + Locate_type lt; int i; + Bounded_side soe = side_of_edge(p, c, lt, i); + if(soe != ON_BOUNDARY) + return soe; - // Undocumented, needed for Mesh_3 (because of Periodic_3_mesh_3) - // \pre c->neighbor(i) is finite - template < class Gt, class Tds, class Lds > - double - Regular_triangulation_3:: - compute_power_distance_to_power_sphere(const Cell_handle c, - const int i) const - { - Cell_handle nc = c->neighbor(i); - CGAL_precondition(!is_infinite(nc)); - Vertex_handle v = nc->vertex(nc->index(c)); + // Either we compare weights, or we use the finite neighboring edge + Cell_handle finite_neighbor = c->neighbor(c->index(infinite_vertex())); + CGAL_triangulation_assertion(!is_infinite(finite_neighbor,0,1)); + return side_of_bounded_power_segment(finite_neighbor->vertex(0)->point(), + finite_neighbor->vertex(1)->point(), + p, perturb); +} - return compute_power_distance_to_power_sphere(c, v); - } +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +greater_or_equal_power_distance(const Bare_point& p, + const Weighted_point& q, + const Weighted_point& r) const +{ + return ! less_power_distance(p, q, r); +} - template < class Gt, class Tds, class Lds > - double - Regular_triangulation_3:: - compute_power_distance_to_power_sphere(const Cell_handle c, - const Vertex_handle v) const - { - typedef typename Geom_traits::Compute_power_distance_to_power_sphere_3 Critical_radius; +// @fixme, below can't be defined there because RegularTriangulationTraits_3 +// does not require the following construction - Critical_radius critical_radius = +// Undocumented, needed for Mesh_3 (because of Periodic_3_mesh_3) +// \pre c->neighbor(i) is finite +template < class Gt, class Tds, class Lds > +double +Regular_triangulation_3:: +compute_power_distance_to_power_sphere(const Cell_handle c, + const int i) const +{ + Cell_handle nc = c->neighbor(i); + CGAL_precondition(!is_infinite(nc)); + Vertex_handle v = nc->vertex(nc->index(c)); + + return compute_power_distance_to_power_sphere(c, v); +} + +template < class Gt, class Tds, class Lds > +double +Regular_triangulation_3:: +compute_power_distance_to_power_sphere(const Cell_handle c, + const Vertex_handle v) const +{ + typedef typename Geom_traits::Compute_power_distance_to_power_sphere_3 Critical_radius; + + Critical_radius critical_radius = geom_traits().compute_power_distance_to_power_sphere_3_object(); - return CGAL::to_double(critical_radius(point(c, 0), point(c, 1), - point(c, 2), point(c, 3), - point(v))); - } + return CGAL::to_double(critical_radius(point(c, 0), point(c, 1), + point(c, 2), point(c, 3), + point(v))); +} - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_Gabriel(const Facet& f) const - { - return is_Gabriel(f.first, f.second); - } +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_Gabriel(const Facet& f) const +{ + return is_Gabriel(f.first, f.second); +} - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_Gabriel(Cell_handle c, int i) const - { - CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i)); - typename Geom_traits::Power_side_of_bounded_power_sphere_3 +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_Gabriel(Cell_handle c, int i) const +{ + CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i)); + typename Geom_traits::Power_side_of_bounded_power_sphere_3 side_of_bounded_orthogonal_sphere = geom_traits().power_side_of_bounded_power_sphere_3_object(); - if ((!is_infinite(c->vertex(i))) && - side_of_bounded_orthogonal_sphere( - c->vertex(vertex_triple_index(i,0))->point(), - c->vertex(vertex_triple_index(i,1))->point(), - c->vertex(vertex_triple_index(i,2))->point(), - c->vertex(i)->point()) == ON_BOUNDED_SIDE ) return false; + if((!is_infinite(c->vertex(i))) && + side_of_bounded_orthogonal_sphere(c->vertex(vertex_triple_index(i,0))->point(), + c->vertex(vertex_triple_index(i,1))->point(), + c->vertex(vertex_triple_index(i,2))->point(), + c->vertex(i)->point()) == ON_BOUNDED_SIDE) + return false; - Cell_handle neighbor = c->neighbor(i); - int in = neighbor->index(c); + Cell_handle neighbor = c->neighbor(i); + int in = neighbor->index(c); - if ((!is_infinite(neighbor->vertex(in))) && - side_of_bounded_orthogonal_sphere( - c->vertex(vertex_triple_index(i,0))->point(), - c->vertex(vertex_triple_index(i,1))->point(), - c->vertex(vertex_triple_index(i,2))->point(), - neighbor->vertex(in)->point()) == ON_BOUNDED_SIDE ) return false; + if((!is_infinite(neighbor->vertex(in))) && + side_of_bounded_orthogonal_sphere(c->vertex(vertex_triple_index(i,0))->point(), + c->vertex(vertex_triple_index(i,1))->point(), + c->vertex(vertex_triple_index(i,2))->point(), + neighbor->vertex(in)->point()) == ON_BOUNDED_SIDE) + return false; - return true; - } + return true; +} +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_Gabriel(const Edge& e) const +{ + return is_Gabriel(e.first, e.second, e.third); +} - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_Gabriel(const Edge& e) const - { - return is_Gabriel(e.first, e.second, e.third); - } - - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_Gabriel(Cell_handle c, int i, int j) const - { - CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i,j)); - typename Geom_traits::Power_side_of_bounded_power_sphere_3 +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_Gabriel(Cell_handle c, int i, int j) const +{ + CGAL_triangulation_precondition(dimension() == 3 && !is_infinite(c,i,j)); + typename Geom_traits::Power_side_of_bounded_power_sphere_3 side_of_bounded_orthogonal_sphere = geom_traits().power_side_of_bounded_power_sphere_3_object(); - Facet_circulator fcirc = incident_facets(c,i,j), - fdone(fcirc); - Vertex_handle v1 = c->vertex(i); - Vertex_handle v2 = c->vertex(j); - do { - // test whether the vertex of cc opposite to *fcirc - // is inside the sphere defined by the edge e = (s, i,j) - Cell_handle cc = (*fcirc).first; - int ii = (*fcirc).second; - if (!is_infinite(cc->vertex(ii)) && - side_of_bounded_orthogonal_sphere( v1->point(), - v2->point(), - cc->vertex(ii)->point()) - == ON_BOUNDED_SIDE ) return false; - } while(++fcirc != fdone); - return true; - } - - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_Gabriel(Vertex_handle v) const + Facet_circulator fcirc = incident_facets(c,i,j), fdone(fcirc); + Vertex_handle v1 = c->vertex(i); + Vertex_handle v2 = c->vertex(j); + do { - return nearest_power_vertex( - geom_traits().construct_point_3_object()(v->point()), v->cell()) == v; + // test whether the vertex of cc opposite to *fcirc + // is inside the sphere defined by the edge e = (s, i,j) + Cell_handle cc = (*fcirc).first; + int ii = (*fcirc).second; + if(!is_infinite(cc->vertex(ii)) && + side_of_bounded_orthogonal_sphere(v1->point(), + v2->point(), + cc->vertex(ii)->point()) == ON_BOUNDED_SIDE) + return false; } + while(++fcirc != fdone); - // Returns - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - insert(const Weighted_point & p, Cell_handle start, bool *could_lock_zone) + return true; +} + +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_Gabriel(Vertex_handle v) const +{ + return nearest_power_vertex( + geom_traits().construct_point_3_object()(v->point()), v->cell()) == v; +} + +// Returns +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert(const Weighted_point& p, Cell_handle start, bool *could_lock_zone) +{ + Locate_type lt; + int li, lj; + + // Parallel + if(could_lock_zone) { - Locate_type lt; - int li, lj; - - // Parallel - if (could_lock_zone) - { - Cell_handle c = locate(p, lt, li, lj, start, could_lock_zone); - if (*could_lock_zone) - return insert(p, lt, c, li, lj, could_lock_zone); - else - return Vertex_handle(); - } - // Sequential + Cell_handle c = locate(p, lt, li, lj, start, could_lock_zone); + if(*could_lock_zone) + return insert(p, lt, c, li, lj, could_lock_zone); else - { - Cell_handle c = locate(p, lt, li, lj, start); - return insert(p, lt, c, li, lj); - } + return Vertex_handle(); } - - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - insert(const Weighted_point & p, Locate_type lt, Cell_handle c, - int li, int lj, bool *could_lock_zone) + // Sequential + else + { + Cell_handle c = locate(p, lt, li, lj, start); + return insert(p, lt, c, li, lj); + } +} + +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert(const Weighted_point& p, Locate_type lt, Cell_handle c, + int li, int lj, bool *could_lock_zone) +{ + switch(dimension()) { - switch (dimension()) { case 3: - { - Conflict_tester_3 tester (p, this); - return insert_in_conflict(p, lt,c,li,lj, tester, - get_hidden_point_visitor(), - could_lock_zone); - } + { + Conflict_tester_3 tester(p, this); + return insert_in_conflict(p, lt,c,li,lj, tester, + get_hidden_point_visitor(), + could_lock_zone); + } case 2: - { - Conflict_tester_2 tester (p, this); - return insert_in_conflict(p, lt,c,li,lj, tester, - get_hidden_point_visitor(), - could_lock_zone); - } + { + Conflict_tester_2 tester(p, this); + return insert_in_conflict(p, lt,c,li,lj, tester, + get_hidden_point_visitor(), + could_lock_zone); + } case 1: - { - Conflict_tester_1 tester (p, this); - return insert_in_conflict(p, lt,c,li,lj, tester, - get_hidden_point_visitor(), - could_lock_zone); - } + { + Conflict_tester_1 tester(p, this); + return insert_in_conflict(p, lt,c,li,lj, tester, + get_hidden_point_visitor(), + could_lock_zone); } - - Conflict_tester_0 tester (p, this); - return insert_in_conflict(p, lt,c,li,lj, tester, - get_hidden_point_visitor(), - could_lock_zone); } + Conflict_tester_0 tester(p, this); + return insert_in_conflict(p, lt,c,li,lj, tester, + get_hidden_point_visitor(), + could_lock_zone); +} - template < class Gt, class Tds, class Lds > - template - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - insert_in_hole(const Weighted_point & p, CellIt cell_begin, CellIt cell_end, - Cell_handle begin, int i) +template < class Gt, class Tds, class Lds > +template +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert_in_hole(const Weighted_point& p, CellIt cell_begin, CellIt cell_end, + Cell_handle begin, int i) +{ + CGAL_triangulation_precondition(cell_begin != cell_end); + + get_hidden_point_visitor().process_cells_in_conflict(cell_begin,cell_end); + + Vertex_handle v = Tr_Base::insert_in_hole(p, cell_begin, cell_end, begin, i); + + // Store the hidden points in their new cells and hide vertices that + // have to be hidden + get_hidden_point_visitor().reinsert_vertices(v); + return v; +} + +template < class Gt, class Tds, class Lds > +template +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert_in_hole(const Weighted_point& p, CellIt cell_begin, CellIt cell_end, + Cell_handle begin, int i, Vertex_handle newv) +{ + CGAL_triangulation_precondition(cell_begin != cell_end); + + get_hidden_point_visitor().process_cells_in_conflict(cell_begin,cell_end); + + Vertex_handle v = Tr_Base::insert_in_hole(p, cell_begin, cell_end, begin, i, newv); + + // Store the hidden points in their new cells and hide vertices that + // have to be hidden + get_hidden_point_visitor().reinsert_vertices(v); + return v; +} + +template +template +class Regular_triangulation_3::Vertex_remover +{ + typedef RegularTriangulation_3 Regular; + typedef typename Gt::Weighted_point_3 Weighted_point; + +public: + typedef typename std::vector::iterator Hidden_points_iterator; + + Vertex_remover(Regular &tmp_) : tmp(tmp_) {} + + Regular& tmp; + + void add_hidden_points(Cell_handle ch) { - CGAL_triangulation_precondition(cell_begin != cell_end); - - get_hidden_point_visitor().process_cells_in_conflict(cell_begin,cell_end); - - Vertex_handle v = - Tr_Base::insert_in_hole(p, cell_begin, cell_end, begin, i); - - // Store the hidden points in their new cells and hide vertices that - // have to be hidden - get_hidden_point_visitor().reinsert_vertices(v); - return v; + std::copy(ch->hidden_points_begin(), ch->hidden_points_end(), + std::back_inserter(hidden)); } + Hidden_points_iterator hidden_points_begin() { return hidden.begin(); } + Hidden_points_iterator hidden_points_end() { return hidden.end(); } - template < class Gt, class Tds, class Lds > - template - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - insert_in_hole(const Weighted_point & p, CellIt cell_begin, CellIt cell_end, - Cell_handle begin, int i, Vertex_handle newv) + Bounded_side side_of_bounded_circle(const Weighted_point& p, const Weighted_point& q, + const Weighted_point& r, const Weighted_point& s, bool perturb = false) const { - CGAL_triangulation_precondition(cell_begin != cell_end); - - get_hidden_point_visitor().process_cells_in_conflict(cell_begin,cell_end); - - Vertex_handle v = - Tr_Base::insert_in_hole(p, cell_begin, cell_end, begin, i, newv); - - // Store the hidden points in their new cells and hide vertices that - // have to be hidden - get_hidden_point_visitor().reinsert_vertices(v); - return v; + return tmp.side_of_bounded_power_circle(p,q,r,s,perturb); } - template - template - class Regular_triangulation_3::Vertex_remover { - typedef RegularTriangulation_3 Regular; - typedef typename Gt::Weighted_point_3 Weighted_point; - public: - typedef typename std::vector::iterator - Hidden_points_iterator; +private: + // The removal of v may un-hide some points, + // Space functions output them. + std::vector hidden; +}; - Vertex_remover(Regular &tmp_) : tmp(tmp_) {} +// The displacement method works only +// on regular triangulation without hidden points at any time +// the vertex inserter is used only +// for the purpose of displacements +template +template +class Regular_triangulation_3::Vertex_inserter +{ + typedef RegularTriangulation_3 Regular; - Regular &tmp; +public: + typedef Nullptr_t Hidden_points_iterator; - void add_hidden_points(Cell_handle ch) { - std::copy (ch->hidden_points_begin(), ch->hidden_points_end(), - std::back_inserter(hidden)); - } + Vertex_inserter(Regular &tmp_) : tmp(tmp_) {} - Hidden_points_iterator hidden_points_begin() { - return hidden.begin(); - } - Hidden_points_iterator hidden_points_end() { - return hidden.end(); - } + Regular& tmp; - Bounded_side side_of_bounded_circle(const Weighted_point &p, const Weighted_point &q, - const Weighted_point &r, const Weighted_point &s, bool perturb = false) const { - return tmp.side_of_bounded_power_circle(p,q,r,s,perturb); - } + void add_hidden_points(Cell_handle) {} + Hidden_points_iterator hidden_points_begin() { return NULL; } + Hidden_points_iterator hidden_points_end() { return NULL; } - private: - // The removal of v may un-hide some points, - // Space functions output them. - std::vector hidden; - }; - - // The displacement method works only - // on regular triangulation without hidden points at any time - // the vertex inserter is used only - // for the purpose of displacements - template - template - class Regular_triangulation_3::Vertex_inserter { - typedef RegularTriangulation_3 Regular; - public: - typedef Nullptr_t Hidden_points_iterator; - - Vertex_inserter(Regular &tmp_) : tmp(tmp_) {} - - Regular &tmp; - - void add_hidden_points(Cell_handle) {} - Hidden_points_iterator hidden_points_begin() { return NULL; } - Hidden_points_iterator hidden_points_end() { return NULL; } - - Vertex_handle insert(const Weighted_point& p, - Locate_type lt, Cell_handle c, int li, int lj) { - return tmp.insert(p, lt, c, li, lj); - } - - Vertex_handle insert(const Weighted_point& p, Cell_handle c) { - return tmp.insert(p, c); - } - - Vertex_handle insert(const Weighted_point& p) { - return tmp.insert(p); - } - }; - - template < class Gt, class Tds, class Lds > - void - Regular_triangulation_3:: - remove(Vertex_handle v) + Vertex_handle insert(const Weighted_point& p, + Locate_type lt, Cell_handle c, int li, int lj) { - Cell_handle c; - if (dimension() > 0) - c = v->cell()->neighbor(v->cell()->index(v)); + return tmp.insert(p, lt, c, li, lj); + } + + Vertex_handle insert(const Weighted_point& p, Cell_handle c) + { + return tmp.insert(p, c); + } + + Vertex_handle insert(const Weighted_point& p) + { + return tmp.insert(p); + } +}; + +template < class Gt, class Tds, class Lds > +void +Regular_triangulation_3:: +remove(Vertex_handle v) +{ + Cell_handle c; + if(dimension() > 0) + c = v->cell()->neighbor(v->cell()->index(v)); + + Self tmp; + Vertex_remover remover(tmp); + Tr_Base::remove(v,remover); + + // Re-insert the points that v was hiding. + for(typename Vertex_remover::Hidden_points_iterator + hi = remover.hidden_points_begin(); + hi != remover.hidden_points_end(); ++hi) + { + Vertex_handle hv = insert(*hi, c); + if(hv != Vertex_handle()) + c = hv->cell(); + } + CGAL_triangulation_expensive_postcondition(is_valid()); +} + +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +remove(Vertex_handle v, bool *could_lock_zone) +{ + bool removed = true; + + // Locking vertex v... + if(!this->try_lock_vertex(v)) + { + *could_lock_zone = false; + } + else + { + Vertex_handle hint = (v->cell()->vertex(0) == v ? + v->cell()->vertex(1) : v->cell()->vertex(0)); Self tmp; Vertex_remover remover(tmp); - Tr_Base::remove(v,remover); + removed = Tr_Base::remove(v, remover, could_lock_zone); - // Re-insert the points that v was hiding. - for (typename Vertex_remover::Hidden_points_iterator - hi = remover.hidden_points_begin(); - hi != remover.hidden_points_end(); ++hi) { - Vertex_handle hv = insert (*hi, c); - if (hv != Vertex_handle()) c = hv->cell(); - } - CGAL_triangulation_expensive_postcondition (is_valid()); - } - - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - remove(Vertex_handle v, bool *could_lock_zone) - { - bool removed = true; - - // Locking vertex v... - if (!this->try_lock_vertex(v)) + if(*could_lock_zone && removed) { - *could_lock_zone = false; - } - else - { - Vertex_handle hint = (v->cell()->vertex(0) == v ? - v->cell()->vertex(1) : v->cell()->vertex(0)); - - Self tmp; - Vertex_remover remover(tmp); - removed = Tr_Base::remove(v, remover, could_lock_zone); - - if (*could_lock_zone && removed) + // Re-insert the points that v was hiding. + for(typename Vertex_remover::Hidden_points_iterator + hi = remover.hidden_points_begin(); + hi != remover.hidden_points_end(); ++hi) { - // Re-insert the points that v was hiding. - for (typename Vertex_remover::Hidden_points_iterator - hi = remover.hidden_points_begin(); - hi != remover.hidden_points_end(); ++hi) + bool could_lock_zone = false; + Vertex_handle hv; + while(!could_lock_zone) { - bool could_lock_zone = false; - Vertex_handle hv; - while (!could_lock_zone) - { - hv = insert (*hi, hint, &could_lock_zone); - } - if (hv != Vertex_handle()) - hint = hv; + hv = insert(*hi, hint, &could_lock_zone); } - CGAL_triangulation_expensive_postcondition (is_valid()); + + if(hv != Vertex_handle()) + hint = hv; } + CGAL_triangulation_expensive_postcondition(is_valid()); } - - return removed; } - // Again, verbatim copy from Delaunay. - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - move_point(Vertex_handle v, const Weighted_point & p) + return removed; +} + +// Again, verbatim copy from Delaunay. +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +move_point(Vertex_handle v, const Weighted_point& p) +{ + CGAL_triangulation_precondition(! is_infinite(v)); + CGAL_triangulation_expensive_precondition(is_vertex(v)); + + // Dummy implementation for a start. + + // Remember an incident vertex to restart + // the point location after the removal. + Cell_handle c = v->cell(); + Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0); + CGAL_triangulation_assertion(old_neighbor != v); + + remove(v); + + if(dimension() <= 0) + return insert(p); + + return insert(p, old_neighbor->cell()); +} + +// Displacement works only for regular triangulation +// without hidden points at any time +template < class Gt, class Tds, class Lds > +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +move_if_no_collision(Vertex_handle v, const Weighted_point& p) +{ + Self tmp; + Vertex_remover remover(tmp); + Vertex_inserter inserter(*this); + Vertex_handle res = Tr_Base::move_if_no_collision(v,p,remover,inserter); + + CGAL_triangulation_expensive_postcondition(is_valid()); + return res; +} + +template +typename Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +move(Vertex_handle v, const Weighted_point& p) +{ + CGAL_triangulation_precondition(!is_infinite(v)); + if(v->point() == p) + return v; + + Self tmp; + Vertex_remover remover(tmp); + Vertex_inserter inserter(*this); + return Tr_Base::move(v,p,remover,inserter); +} + +template < class Gt, class Tds, class Lds > +bool +Regular_triangulation_3:: +is_valid(bool verbose, int level) const +{ + if(! Tr_Base::is_valid(verbose,level)) { - CGAL_triangulation_precondition(! is_infinite(v)); - CGAL_triangulation_expensive_precondition(is_vertex(v)); + if(verbose) + std::cerr << "invalid base triangulation" << std::endl; - // Dummy implementation for a start. - - // Remember an incident vertex to restart - // the point location after the removal. - Cell_handle c = v->cell(); - Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0); - CGAL_triangulation_assertion(old_neighbor != v); - - remove(v); - - if (dimension() <= 0) - return insert(p); - return insert(p, old_neighbor->cell()); + CGAL_triangulation_assertion(false); + return false; } - // Displacement works only for regular triangulation - // without hidden points at any time - template < class Gt, class Tds, class Lds > - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - move_if_no_collision(Vertex_handle v, const Weighted_point &p) + switch(dimension()) { - Self tmp; - Vertex_remover remover (tmp); - Vertex_inserter inserter (*this); - Vertex_handle res = Tr_Base::move_if_no_collision(v,p,remover,inserter); - - CGAL_triangulation_expensive_postcondition(is_valid()); - return res; - } - - template - typename Regular_triangulation_3::Vertex_handle - Regular_triangulation_3:: - move(Vertex_handle v, const Weighted_point &p) { - CGAL_triangulation_precondition(!is_infinite(v)); - if(v->point() == p) return v; - Self tmp; - Vertex_remover remover (tmp); - Vertex_inserter inserter (*this); - return Tr_Base::move(v,p,remover,inserter); - } - - template < class Gt, class Tds, class Lds > - bool - Regular_triangulation_3:: - is_valid(bool verbose, int level) const - { - if ( ! Tr_Base::is_valid(verbose,level) ) { - if (verbose) - std::cerr << "invalid base triangulation" << std::endl; - CGAL_triangulation_assertion(false); - return false; - } - - switch ( dimension() ) { case 3: + { + for(Finite_cells_iterator it = finite_cells_begin(), + end = finite_cells_end(); it != end; ++it) { - for(Finite_cells_iterator it = finite_cells_begin(), end = finite_cells_end(); it != end; ++it) { - is_valid_finite(it, verbose, level); - for(int i=0; i<4; i++) { - if ( !is_infinite - (it->neighbor(i)->vertex(it->neighbor(i)->index(it))) ) { - if ( side_of_power_sphere - (it, - it->neighbor(i)->vertex(it->neighbor(i)->index(it))->point()) - == ON_BOUNDED_SIDE ) { - if (verbose) - std::cerr << "non-empty sphere " << std::endl; - CGAL_triangulation_assertion(false); - return false; - } + is_valid_finite(it, verbose, level); + for(int i=0; i<4; i++) + { + if(!is_infinite(it->neighbor(i)->vertex(it->neighbor(i)->index(it)))) + { + if(side_of_power_sphere(it, + it->neighbor(i)->vertex( + it->neighbor(i)->index(it))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) + std::cerr << "non-empty sphere " << std::endl; + + CGAL_triangulation_assertion(false); + return false; } } } - break; - } - case 2: - { - for(Finite_facets_iterator it = finite_facets_begin(), end = finite_facets_end(); it!= end; ++it) { - is_valid_finite((*it).first, verbose, level); - for(int i=0; i<3; i++) { - if( !is_infinite - ((*it).first->neighbor(i)->vertex( (((*it).first)->neighbor(i)) - ->index((*it).first))) ) { - if ( side_of_power_circle - ( (*it).first, 3, - (*it).first->neighbor(i)-> - vertex( (((*it).first)->neighbor(i)) - ->index((*it).first) )->point() ) - == ON_BOUNDED_SIDE ) { - if (verbose) - std::cerr << "non-empty circle " << std::endl; - CGAL_triangulation_assertion(false); - return false; - } - } - } - } - break; - } - case 1: - { - for(Finite_edges_iterator it = finite_edges_begin(), end = finite_edges_end(); it != end; ++it) { - is_valid_finite((*it).first, verbose, level); - for(int i=0; i<2; i++) { - if( !is_infinite - ((*it).first->neighbor(i)->vertex( (((*it).first)->neighbor(i)) - ->index((*it).first))) ) { - if ( side_of_power_segment - ( (*it).first, - (*it).first->neighbor(i)-> - vertex( (((*it).first)->neighbor(i)) - ->index((*it).first) )->point() ) - == ON_BOUNDED_SIDE ) { - if (verbose) - std::cerr << "non-empty edge " << std::endl; - CGAL_triangulation_assertion(false); - return false; - } - } - } - } - break; } + break; + } + case 2: + { + for(Finite_facets_iterator it = finite_facets_begin(), end = finite_facets_end(); it!= end; ++it) + { + is_valid_finite((*it).first, verbose, level); + for(int i=0; i<3; i++) + { + if(!is_infinite((*it).first->neighbor(i)->vertex( + (((*it).first)->neighbor(i))->index((*it).first)))) + { + if(side_of_power_circle((*it).first, 3, + (*it).first->neighbor(i)-> + vertex((((*it).first)->neighbor(i)) + ->index((*it).first))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) + std::cerr << "non-empty circle " << std::endl; + + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + break; + } + case 1: + { + for(Finite_edges_iterator it = finite_edges_begin(), + end = finite_edges_end(); it != end; ++it) + { + is_valid_finite((*it).first, verbose, level); + for(int i=0; i<2; i++) + { + if(!is_infinite + ((*it).first->neighbor(i)->vertex( + (((*it).first)->neighbor(i))->index((*it).first)))) + { + if(side_of_power_segment((*it).first, + (*it).first->neighbor(i)->vertex( + (((*it).first)->neighbor(i))->index( + (*it).first))->point()) == ON_BOUNDED_SIDE) + { + if(verbose) + std::cerr << "non-empty edge " << std::endl; + + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + break; } - if (verbose) - std::cerr << "valid regular triangulation" << std::endl; - return true; } + if(verbose) + std::cerr << "valid regular triangulation" << std::endl; + + return true; +} + } //namespace CGAL #if defined(BOOST_MSVC) diff --git a/Triangulation_3/include/CGAL/Triangulation_3.h b/Triangulation_3/include/CGAL/Triangulation_3.h index 68c6d7cc8d5..d2299eb694e 100644 --- a/Triangulation_3/include/CGAL/Triangulation_3.h +++ b/Triangulation_3/include/CGAL/Triangulation_3.h @@ -201,8 +201,8 @@ class Triangulation_3_base public: // If Lock_data_structure_ = Default => use Spatial_lock_grid_3 typedef typename Default::Get< - Lock_data_structure_, Spatial_lock_grid_3 > - ::type Lock_data_structure; + Lock_data_structure_, + Spatial_lock_grid_3 >::type Lock_data_structure; protected: Triangulation_3_base() @@ -237,9 +237,9 @@ protected: struct Vertex_handle_unique_hash_map_generator { typedef Unique_hash_map > type; + Vertex_handle, + Handle_hash_function, + tbb::scalable_allocator > type; }; public: