diff --git a/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h b/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h index 9d6d449cd41..da6aa2538cf 100644 --- a/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h +++ b/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h @@ -238,6 +238,18 @@ public: public: Vertex_handle insert(const Point& p, Face_handle f = Face_handle()); + + // For convenience, when P3 != PoS2 + template + Vertex_handle insert(const P& p, + Face_handle f = Face_handle(), + typename std::enable_if::value>::type* = nullptr) + { + CGAL_triangulation_assertion((std::is_same::value)); + + return insert(geom_traits().construct_point_on_sphere_2_object()(p), f); + } + Vertex_handle push_back(const Point& p, Face_handle f = Face_handle()) { return insert(p, f); } Vertex_handle insert(const Point& p, Locate_type lt, Face_handle loc, int li); Vertex_handle insert_first(const Point& p); @@ -246,37 +258,25 @@ public: Vertex_handle insert_outside_affine_hull_regular(const Point& p); Vertex_handle insert_cocircular(const Point& p, Locate_type lt, Face_handle loc); - template - int insert(InputIterator first, InputIterator beyond, - typename std::enable_if< - std::is_same::value_type, - Point>::value>::type* = nullptr); - - // For convenience when P3 != PoS2 - template - Vertex_handle insert(const P& p, - Face_handle f = Face_handle(), - typename std::enable_if::value>::type* = nullptr) - { - return insert(geom_traits().construct_point_on_sphere_2_object()(p), f); - } + // Range insertion + // Input range has value type Point, with Point != Point_3 template int insert(InputIterator first, InputIterator beyond, typename std::enable_if< !std::is_same::value_type, - Point>::value>::type* = nullptr) - { - typename Geom_traits::Construct_point_on_sphere_2 cst = geom_traits().construct_point_on_sphere_2_object(); + Point_3>::value>::type* = nullptr); - return insert(boost::make_transform_iterator(first, cst), - boost::make_transform_iterator(beyond, cst)); - } + // Input range has value type Point_3, possibly with Point == Point_3 + template + int insert(InputIterator first, InputIterator beyond, + typename std::enable_if< + std::is_same::value_type, + Point_3>::value>::type* = nullptr); - bool update_ghost_faces(Vertex_handle v = Vertex_handle(), bool first = false); + bool update_ghost_faces(Vertex_handle v, bool first = false); - //REMOVAL - void remove_degree_3(Vertex_handle v, Face_handle f = Face_handle()); + // Removal void remove_1D(Vertex_handle v); void remove_2D(Vertex_handle v); void remove(Vertex_handle v); @@ -367,36 +367,6 @@ test_conflict(const Point& p, Face_handle fh) const } // ------------------------ INSERTION --------------------------------// -template -typename Delaunay_triangulation_on_sphere_2::Vertex_handle -Delaunay_triangulation_on_sphere_2:: -insert(const Point& p, Face_handle start) -{ - Locate_type lt; - int li; - Face_handle loc = Base::locate(p, lt, li, start); - - switch(lt) - { - case NOT_ON_SPHERE: - return Vertex_handle(); - case TOO_CLOSE: - { - CGAL_assertion(loc != Face_handle()); - return loc->vertex(li); // @todo test this for dim < 2 - } - case VERTEX: - { - if(number_of_vertices() == 1) - return vertices_begin(); - - return (loc->vertex(li)); - } - default: // the point can be inserted - return insert(p, lt, loc, li); - } -} - // inserts a new point to a 1d triangulation, the new point is also coplanar with the existing points. template typename Delaunay_triangulation_on_sphere_2::Vertex_handle @@ -426,7 +396,8 @@ insert_cocircular(const Point& p, tds().delete_face(loc); - update_ghost_faces(v); + CGAL_postcondition(dimension() == 1); + return v; } @@ -484,7 +455,6 @@ insert_third(const Point& p) point(f, 1), point(f->neighbor(0), 1)) != RIGHT_TURN); - update_ghost_faces(nv); return nv; } @@ -579,6 +549,128 @@ insert(const Point& p, Locate_type lt, Face_handle loc, int /*li*/) return v; } +template +typename Delaunay_triangulation_on_sphere_2::Vertex_handle +Delaunay_triangulation_on_sphere_2:: +insert(const Point& p, Face_handle start) +{ + Locate_type lt; + int li; + Face_handle loc = Base::locate(p, lt, li, start); + + switch(lt) + { + case NOT_ON_SPHERE: + return Vertex_handle(); + case TOO_CLOSE: + { + CGAL_triangulation_assertion(loc != Face_handle()); + return loc->vertex(li); + } + case VERTEX: + { + if(number_of_vertices() == 1) + return vertices_begin(); + + return loc->vertex(li); + } + default: // the point can be inserted + return insert(p, lt, loc, li); + } +} + +template +template +int +Delaunay_triangulation_on_sphere_2:: +insert(InputIterator first, InputIterator beyond, + typename std::enable_if< + !std::is_same::value_type, + Point_3>::value>::type*) +{ + CGAL_static_assertion((std::is_same::value_type, Point>::value)); + CGAL_static_assertion(!(std::is_same::value)); + + const int n = number_of_vertices(); + + // On paper, Spatial_sort_traits_adapter_3 should be used. However, it is not compatible + // with spatial_sort_on_sphere() because of the way Transform_coordinates_traits_3 were written: + // there are hard calls to x(), y(), and z() whereas it should have been calls to Compute_x_3 & such. + + std::vector input_points(first, beyond); + + struct Point_3_with_iterator + : public Point_3 + { + Point_3_with_iterator(const Point& p, Self const * const tr) + : Point_3(tr->geom_traits().construct_point_3_object()(p)), input_point_ptr(&p) + { } + + const Point* input_point_ptr; + }; + + std::vector p3_points; + p3_points.reserve(std::distance(first, beyond)); + for(const Point& p : input_points) + p3_points.push_back(Point_3_with_iterator(p, this)); + + CGAL::cpp98::random_shuffle(p3_points.begin(), p3_points.end()); + + // @todo: + // - bench this + // - should points not on the sphere already be filtered because it might mess up the sort? + spatial_sort_on_sphere(p3_points.begin(), p3_points.end(), + geom_traits(), square(geom_traits().radius()), geom_traits().center()); + + // Reordering could be done in-place, but not worth the hassle for now + std::vector points; + points.reserve(input_points.size()); + for(const Point_3_with_iterator& p3wi : p3_points) + { + CGAL_triangulation_assertion(p3wi.input_point_ptr != nullptr); + points.push_back(*(p3wi.input_point_ptr)); + } + + Face_handle hint; + Vertex_handle v; + for(const Point& p : points) + { + v = insert(p, hint); + if(v != Vertex_handle()) // can happen e.g. if the point is not on the sphere + hint = v->face(); + } + + return number_of_vertices() - n; +} + +template +template +int +Delaunay_triangulation_on_sphere_2:: +insert(InputIterator first, InputIterator beyond, + typename std::enable_if< + std::is_same::value_type, + Point_3>::value>::type*) +{ + const int n = number_of_vertices(); + + std::vector points(first, beyond); + CGAL::cpp98::random_shuffle(points.begin(), points.end()); + + spatial_sort_on_sphere(points.begin(), points.end(), geom_traits(), + square(geom_traits().radius()), geom_traits().center()); + + Face_handle hint; + for(const Point_3& p : points) + { + Vertex_handle v = insert(p, hint); + if(v != Vertex_handle()) // could happen if the point is not on the sphere + hint = v->face(); + } + + return number_of_vertices() - n; +} + /* * method to test and mark faces incident to v as ghost-faces or solid-faces. * the Boolean 'first' defines whether the dimension of the triangulation increased @@ -589,152 +681,57 @@ bool Delaunay_triangulation_on_sphere_2:: update_ghost_faces(Vertex_handle v, bool first) { - if(number_of_vertices() < 3) - return false; - - CGAL_assertion(dimension() >= 1); + CGAL_triangulation_assertion(dimension() == 2); bool ghost_found = false; - if(dimension() == 1) + + if(first) // first time dimension 2 { - All_edges_iterator eit = all_edges_begin(); - for(; eit!=all_edges_end(); ++eit) + for(All_faces_iterator fi=all_faces_begin(); fi!=all_faces_end(); ++fi) { - Face_handle f = eit->first; - Face_handle fn = f->neighbor(0); - const Point& q = point(fn, 1); - if(this->collinear_between(point(f, 0), point(f, 1), q)) + if(orientation_on_sphere(fi) != POSITIVE) { - f->set_ghost(true); + fi->set_ghost(true); ghost_found = true; } else { - f->set_ghost(false); + fi->set_ghost(false); } } } - else // dimension==2 + else // not first { - if(first) // first time dimension 2 + CGAL_triangulation_assertion(v != Vertex_handle()); + Face_circulator fc = this->incident_faces(v, v->face()); + Face_circulator done(fc); + do { - All_faces_iterator fi = all_faces_begin(); - for(; fi!=all_faces_end(); fi++) + if(orientation_on_sphere(fc) != POSITIVE) { - if(orientation_on_sphere(fi) != POSITIVE) - { - fi->set_ghost(true); - ghost_found = true; - } - else - { - fi->set_ghost(false); - } + fc->set_ghost(true); + ghost_found = true; + } + else + { + fc->set_ghost(false); } } - else // not first - { - CGAL_assertion(v != Vertex_handle()); - Face_circulator fc = this->incident_faces(v, v->face()); - Face_circulator done(fc); - do - { - if(orientation_on_sphere(fc) != POSITIVE) - { - fc->set_ghost(true); - ghost_found = true; - } - else - { - fc->set_ghost(false); - } - } - while(++fc != done); - } + while(++fc != done); } return ghost_found; } -template -template -int -Delaunay_triangulation_on_sphere_2:: -insert(InputIterator first, InputIterator beyond, - typename std::enable_if< - std::is_same::value_type, Point>::value>::type*) -{ - int n = number_of_vertices(); - - std::vector points(first, beyond); - CGAL::cpp98::random_shuffle(points.begin(), points.end()); - - typename Geom_traits::Construct_point_3 cp3 = geom_traits().construct_point_3_object(); - - typedef typename Geom_traits::Construct_point_3 Construct_point_3; - typedef typename boost::result_of::type Ret; - typedef boost::function_property_map fpmap; - typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; - - // @todo: - // - bench this - // - should points not on the sphere already be filtered because it might mess up the sort? - spatial_sort_on_sphere(points.begin(), points.end(), - Search_traits_3( - boost::make_function_property_map( - geom_traits().construct_point_3_object()), geom_traits()), - square(geom_traits().radius()), geom_traits().center()); - - Face_handle hint; - Vertex_handle v; - for(typename std::vector::const_iterator p=points.begin(), end=points.end(); p!=end; ++p) - { - v = insert(*p, hint); - if(v != Vertex_handle()) // could happen if the point is not on the sphere - hint = v->face(); - } - - return number_of_vertices() - n; -} - //-------------------------------------REMOVAL----------------------------------------------------// -template -void -Delaunay_triangulation_on_sphere_2:: -remove_degree_3(Vertex_handle v, Face_handle f) -{ - CGAL_precondition(v != Vertex_handle()); - - if(f == Face_handle()) - f = v->face(); - - tds().remove_degree_3(v, f); -} - -template -void -Delaunay_triangulation_on_sphere_2:: -remove(Vertex_handle v) -{ - CGAL_precondition(v != Vertex_handle()); - - if(number_of_vertices() <= 3) - tds().remove_dim_down(v); - else if(dimension() == 2) - remove_2D(v); - else - remove_1D(v); -} - template void Delaunay_triangulation_on_sphere_2:: remove_1D(Vertex_handle v) { - CGAL_precondition(v != Vertex_handle()); + CGAL_triangulation_precondition(v != Vertex_handle()); tds().remove_1D(v); - update_ghost_faces(); } template @@ -747,7 +744,6 @@ remove_2D(Vertex_handle v) if(test_dim_down(v)) // resulting triangulation has dim 1 { tds().remove_dim_down(v); - update_ghost_faces(); //1d triangulation, no vertex needed to update ghost-faces } else { @@ -757,6 +753,21 @@ remove_2D(Vertex_handle v) } } +template +void +Delaunay_triangulation_on_sphere_2:: +remove(Vertex_handle v) +{ + CGAL_triangulation_precondition(v != Vertex_handle()); + + if(number_of_vertices() <= 3) + tds().remove_dim_down(v); + else if(dimension() == 2) + remove_2D(v); + else + remove_1D(v); +} + // tests whether the dimension will decrease when removing v from the triangulation. template bool diff --git a/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h b/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h index 6c76587f38e..e8964f6b83e 100644 --- a/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h +++ b/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h @@ -695,22 +695,99 @@ test_distance(const Point& p, } } +/* + * Location for degenerated cases: locates the conflicting edge in a 1-dimensional triangulation. + * This method is used when the new point is coplanar with the existing vertices. + * The Boolean 'on_diametral_plane' indicates whether the points are also coplanar with the + * center of the sphere (true) or not (false). + */ +template +typename Triangulation_on_sphere_2::Face_handle +Triangulation_on_sphere_2:: +locate_edge(const Point& p, + Locate_type& lt, + int& li, + const bool on_diametral_plane) const +{ + CGAL_precondition(dimension() == 1); + + Face_handle loc; + if(on_diametral_plane) + { + All_edges_iterator eit; + for(eit=all_edges_begin(); eit!=all_edges_end(); ++eit) + { + const Face_handle f = eit->first; + const Face_handle fn = f->neighbor(0); + const Point& q = point(fn, 1); + + // Check if *eit is an a "ghost" edge, that is if in the diametral plane, the center + // of the sphere is to the right of (eit->source(); eit->target()). + // Note that this doesn't check if p is actually "in conflict" with that edge, + // so that edge is just kept in memory, and if no solid edge is found, then it has to be this edge + if(collinear_between(point(f, 0), point(f, 1), q)) + { + // The new point is on the same 3D plane as the existing vertices, but not on an existing edge + loc = eit->first; + } + else + { + // not ghost, check the cone + if(collinear_between(point(eit->first, 0), point(eit->first, 1), p)) + { + loc = eit->first; + lt = EDGE; + li = 2; + test_distance(p, loc, lt, li); + return loc; + } + } + } + + // Couldn't find a solid edge that is in conflict with p + lt = OUTSIDE_CONVEX_HULL; + li = 4; + + test_distance(p, loc, lt, li); + return loc; + } + else // not coplanar with the center of the sphere + { + for(All_edges_iterator eit = all_edges_begin(); eit!=all_edges_end(); ++eit) + { + // The plane going through the center of the sphere, v1, and v2 cuts the non-diametral circle + // in two, with the negative side of the plane corresponding to the shorter part, so + // if p is coplanar with all the points on that circle and on the negative side, it splits + // the edge v1-v2 + if(orientation_on_sphere(point(eit->first, 0), point(eit->first, 1), p) == ON_NEGATIVE_SIDE) + { + loc = eit->first; + lt = EDGE; + li = 2; + test_distance(p, loc, lt, li); + return loc; + } + } + + CGAL_assertion(false); + return loc; + } +} + template typename Triangulation_on_sphere_2::Face_handle Triangulation_on_sphere_2:: march_locate_1D(const Point& p, Locate_type& lt, int& li) const { - CGAL_assertion(dimension() == 1); + CGAL_triangulation_assertion(dimension() == 1); - - // check if p is coplanar with existing points - // first three points of triangulation + // Check if p is coplanar with the existing points first three points of the triangulation Face_handle f = all_edges_begin()->first; - Vertex_handle v1 = f->vertex(0); - Vertex_handle v2 = f->vertex(1); - Vertex_handle v3 = f->neighbor(0)->vertex(1); + const Vertex_handle v1 = f->vertex(0); + const Vertex_handle v2 = f->vertex(1); + const Vertex_handle v3 = f->neighbor(0)->vertex(1); - Orientation orient = orientation(point(v1), point(v2), point(v3), p); + const Orientation orient = orientation(point(v1), point(v2), point(v3), p); if(orient != ON_ORIENTED_BOUNDARY) { lt = OUTSIDE_AFFINE_HULL; @@ -719,22 +796,23 @@ march_locate_1D(const Point& p, Locate_type& lt, int& li) const return f; } - // from then on, p is coplanar with all the triangulation's points + // From then on, p is coplanar with all the triangulation's points - // check if p is coradial with one existing point + // Check if p is coradial with one existing point Vertices_iterator vi; for(vi=vertices_begin(); vi!=vertices_end(); ++vi) // @todo turns insertion into O(n) { if(are_equal(point(vi), p)) { lt = VERTEX; - li = 1; + f = vi->face(); + li = f->index(vi); // could be simply '1' return f; } } // v1, v2, and v3 are coplanar so this is just checking if they are coplanar with the center of the sphere - Orientation pqr = orientation_on_sphere(point(v1), point(v2), point(v3)); + const Orientation pqr = orientation_on_sphere(point(v1), point(v2), point(v3)); if(pqr == ON_ORIENTED_BOUNDARY) return locate_edge(p, lt, li, true /*on_diametral_plane*/); else @@ -749,7 +827,8 @@ march_locate_2D(Face_handle f, Locate_type& lt, int& li) const { - CGAL_assertion(!f->is_ghost()); + CGAL_triangulation_precondition(dimension() == 2); + CGAL_triangulation_precondition(!is_ghost(f)); boost::rand48 rng; boost::uniform_smallint<> two(0, 1); @@ -760,7 +839,7 @@ march_locate_2D(Face_handle f, for(;;) { - if(f->is_ghost()) + if(is_ghost(f)) { if(orientation(f, t) == ON_POSITIVE_SIDE) // conflict with the corresponding face { @@ -1004,15 +1083,12 @@ locate(const Point& p, return Face_handle(); } - //std::cout << "Locate " << p << ", dimension " << dimension() << std::endl; - switch(dimension()) { case -2: // empty triangulation { lt = OUTSIDE_AFFINE_HULL; - li = 4; // li should not be used in this case - + li = 4; return Face_handle(); } case -1: // 1 vertex @@ -1026,7 +1102,6 @@ locate(const Point& p, lt = VERTEX; Face_handle f = vi->face(); li = f->index(vi); - return f; } } @@ -1039,22 +1114,21 @@ locate(const Point& p, return f; } - case 1: + case 1: // 3+ coplanar vertices { return march_locate_1D(p, lt, li); } } // below is dimension() == 2 - if(start == Face_handle()) start = all_faces_begin(); - if(start->is_ghost()) + if(is_ghost(start)) { - for(All_faces_iterator it = this->_tds.face_iterator_base_begin(); it !=all_faces_end(); it++) + for(All_faces_iterator it=all_faces_begin(); it !=all_faces_end(); ++it) { - if(!it->is_ghost()) + if(!is_ghost(it)) { start = it; break; @@ -1122,8 +1196,7 @@ show_all() const if(dimension() == 1) { std::cerr << " all edges dim 1 " << std::endl; - All_edges_iterator aeit; - for(aeit=all_edges_begin(); aeit!=all_edges_end(); ++aeit) + for(All_edges_iterator aeit=all_edges_begin(); aeit!=all_edges_end(); ++aeit) { show_face(aeit->first); std::cerr << " ------------ " << std::endl; @@ -1133,8 +1206,7 @@ show_all() const } std::cerr << " faces " << std::endl; - All_faces_iterator fi; - for(fi = all_faces_begin(); fi !=all_faces_end(); ++fi) + for(All_faces_iterator fi = all_faces_begin(); fi !=all_faces_end(); ++fi) { show_face(fi); std::cerr << " ------------ " << std::endl; @@ -1143,8 +1215,7 @@ show_all() const if(number_of_vertices() > 1) { std::cerr << "print triangulation vertices:" << std::endl; - Vertices_iterator vi; - for(vi=vertices_begin(); vi!=vertices_end(); ++vi) + for(Vertices_iterator vi=vertices_begin(); vi!=vertices_end(); ++vi) { show_vertex(vi); std::cerr << " / associated face: " << (void*)(&(*(vi->face()))) << std::endl;;