diff --git a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3DelaunayTriangulationTraits_3.h b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3DelaunayTriangulationTraits_3.h index 64016500753..881ae37a8c7 100644 --- a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3DelaunayTriangulationTraits_3.h +++ b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3DelaunayTriangulationTraits_3.h @@ -34,12 +34,7 @@ public: /// @{ /*! -A predicate object that must provide the function operators - -`Oriented_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Point_3 t)`, - -which determines on which side of the oriented sphere circumscribing -`p, q, r, s` the point `t` lies and +A predicate object that must provide the function operator `Oriented_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Point_3 t, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s, Periodic_3_offset_3 o_t)`, @@ -51,12 +46,7 @@ which determines on which side of the oriented sphere circumscribing typedef unspecified_type Side_of_oriented_sphere_3; /*! -A predicate object that must provide the function operators - -`Comparison_result operator()(Point_3 p, Point_3 q, Point_3 r)`, - -which compares the distance between `p` and `q` to the distance -between `p` and `r` and +A predicate object that must provide the function operator `Comparison_result operator()(Point_3 p, Point_3 q, Point_3 r, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r)`, @@ -73,13 +63,7 @@ typedef unspecified_type Compare_distance_3; /// @{ /*! -A predicate object that must provide the function operators - -`Orientation operator()(Point_3 p, Point_3 q, Point_3 r)`, - -which returns `COLLINEAR`, if the points are collinear; otherwise -it must return a consistent orientation for any three points chosen in -a same plane and +A predicate object that must provide the function operator `Orientation operator()(Point_3 p, Point_3 q, Point_3 r Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r)`, @@ -91,12 +75,7 @@ three point-offset pairs chosen in a same plane. typedef unspecified_type Coplanar_orientation_3; /*! -A predicate object that must provide the function operators - -`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`, - -which determines the bounded side of the circle defined by `p, q`, -and `r` on which the point `s` lies and +A predicate object that must provide the function operator `Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`, @@ -115,42 +94,19 @@ typedef unspecified_type Coplanar_side_of_bounded_circle_3; /// @{ /*! -A predicate object that must provide the function operators - -`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 t)`, - -which returns the position of the point `t` relative to the sphere -that has `pq` as its diameter, +A predicate object that must provide the function operator `Bounded_side operator()(Point_3 p, Point_3 q, Point_3 t, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_t)`, which returns the position of the point-offset pair `(t,o_t)` relative to the sphere that has `(p,o_p)(q,o_q)` as its diameter, -`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 t)`, - -which returns the position of the point `t` relative to the sphere -passing through `p, q`, and `r` and whose center is in the -plane defined by these three points, - `Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 t, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_q)`, which returns the position of the point-offset pair `(t,o_t)` relative to the sphere passing through `(p,o_p), (q,o_q)`, and `(r,o_r)` and whose center is in the plane defined by these three -point-offset pairs, - -`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Point_3 t)`, - -which returns the relative position of point `t` to the sphere -defined by `p, q, r`, and `s`; the order of the points `p, q, r`, and `s` does not matter, and - -`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Point_3 t, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s, Periodic_3_offset_3 o_q)`, - -which returns the relative position of the point-offset pair -`(t,o_t)` to the sphere defined by `(p,o_p), (q,o_q), (r,o_r)`, and `(s,o_s)`; the order of the point-offset pairs -`(p,o_p), (q,o_q), (r,o_r)`, and `(s,o_s)` does not matter. -\pre `p, q, r`, and `s` are not coplanar, `(p,o_p), (q,o_q), (r,o_r)`, and `(s,o_s)` are not coplanar, `p`, `q`, `r`, `s`, `t` lie inside the domain. +point-offset pairs. */ typedef unspecified_type Side_of_bounded_sphere_3; @@ -162,11 +118,7 @@ typedef unspecified_type Side_of_bounded_sphere_3; /// @{ /*! -A constructor object that must provide the function operators - -`Point_3 operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`, - -which constructs the circumcenter of four points and +A constructor object that must provide the function operator `Point_3 operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`, diff --git a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3RegularTriangulationTraits_3.h b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3RegularTriangulationTraits_3.h index 0c2e6e0d919..9f7ffb809eb 100644 --- a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3RegularTriangulationTraits_3.h +++ b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3RegularTriangulationTraits_3.h @@ -34,58 +34,41 @@ public: /// @{ /*! -A predicate object that must provide the function operators: - -`Oriented_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s, Weighted_point_3 t)`, - -which determines the position of `t` with respect to the power sphere of `p, q, r, s`. +A predicate object that must provide the function operators `Oriented_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s, Weighted_point_3 t, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s, Periodic_3_offset_3 o_t)`, -which is the same for the point-offset pair `(t,o_t)` with respect to the power sphere of the point-offset pairs -`(p,o_p), (q,o_q), (r,o_r), (s,o_s)`. +which determines the position of the point-offset pair `(t,o_t)` with respect +to the power sphere of the point-offset pairs `(p,o_p), (q,o_q), (r,o_r), (s,o_s)`. \pre `p`, `q`, `r`, `s`, `t` lie inside the domain and `p, q, r, s` are not coplanar.
-When vertex removal is used, the predicate must in addition provide the following function operators: +When vertex removal is used, the predicate must in addition provide the function operators -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 t)`, +`Oriented_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 t, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_t)`, which has a definition similar to the previous method, for coplanar points, with the power circle of `p,q,r`. +\pre `p`, `q`, `r`, `t` lie inside the domain, `p, q, r` are not collinear, +and `(p,o_p), (q,o_q), (r,o_r), (t,o_t)` are coplanar. -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 t, -Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_t)`, +`Oriented_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 t, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_t)`, -which is the same for point-offset pairs. +which is the same for collinear points, and the power segment of `(p,o_p)` and `(q,o_q)`, -\pre `p`, `q`, `r`, `t` lie inside the domain, `p, q, r` are not collinear, and `p, q, r, t` are coplanar. +\pre `p`, `q`, `t` lie inside the domain, `p` and `q` have different Bare_points, and +`(p,o_p), (q,o_q), (t,o_t)` are collinear. +`Oriented_side operator()(Weighted_point_3 p, Weighted_point_3 q, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`, -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 t)`, - -which is the same for collinear points, and the power segment of `p` and `q`, - -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 t, -Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_t)`, - -which is the same for point-offset pairs. - -\pre `p`, `q`, `t` lie inside the domain, `p` and `q` have different Bare_points, and `p, q, t` are collinear. - - -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q)`, - -which is the same for equal points, that is when `p` and `q` +which is the same for equal points, that is when `(p,o_p)` and `(q,o_q)` have equal coordinates, then it returns the comparison of the weights. -`Oriented_side operator()( Weighted_point_3 p, Weighted_point_3 q, -Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`, - -which is the same for point-offset pairs. - \pre `p` and `q` lie inside the domain and have equal Bare_points. */ @@ -118,15 +101,11 @@ typedef unspecified_type Compare_weighted_squared_radius_3; A predicate object, model of `ComparePowerDistance_3`, that must provide the function operator -`Comparison_result operator()(Point_3 p, Weighted_point_3 q, Weighted_point_3 r)`, - -which compares the power distance between `p` and `q` to the power distance -between `p` and `r` and - `Comparison_result operator()(Point_3 p, Weighted_point_3 q, Weighted_point_3 r, -Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r)`, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r)`, -which is the same for point-offset pairs. +which compares the power distance between `(p,o_p)` and `(q,o_q)` to the power distance +between `(p,o_p)` and `(r,o_r)`. \note This predicate is required if a call to `nearest_power_vertex()` or `nearest_power_vertex_in_cell()` is issued.*/ @@ -157,23 +136,51 @@ typedef unspecified_type Coplanar_orientation_3; /// @} +/// \name +/// When `is_Gabriel` functions are used, the traits class must +/// in addition provide the following predicate object: +/// @{ + +/*! +A predicate object that must provide the function operator + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 t, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_t)`, + +which returns the sign of the power test of `(t,o_t)` with respect to the smallest +sphere orthogonal to `(p,o_p)` (which is the sphere with center `(p,o_p)` and squared +radius `-w_p` with `w_p` the weight of `p`), + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 t, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_t)`, + +which returns the sign of the power test of `(t,o_t)` with respect to the smallest +sphere orthogonal to `(p,o_p)` and `(q,o_q)`, + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 t, + Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_q)`, + +which returns the sign of the power test of `(t,o_t)` with respect to the smallest +sphere orthogonal to `(p,o_p)`, `(q,o_q)`, and `(r,o_r)`. +*/ +typedef unspecified_type Power_side_of_bounded_power_sphere_3; + +/// @} + /// \name /// When the dual operations are used, the traits /// class must in addition provide the following constructor object: /// @{ /*! -A constructor object that must provide the function operators: - -`Weighted_point_3 operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s)`, - -which constructs the weighted circumcenter of four points and +A constructor object that must provide the function operator `Weighted_point_3 operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`, which constructs the weighted circumcenter of four point-offset pairs. -\pre `p`, `q`, `r` and `s` as well as `(p,o_p)`, `(q,o_q)`, `(r,o_r)` and `(s,o_s)` must be non coplanar. `p`, `q`, `r`, `s` lie inside the domain. +\pre `p`, `q`, `r`, `s` lie inside the domain. `p`, `q`, `r` and `s`, +as well as `(p,o_p)`, `(q,o_q)`, `(r,o_r)` and `(s,o_s)` must be non coplanar. */ typedef unspecified_type Construct_weighted_circumcenter_3; diff --git a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3TriangulationTraits_3.h b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3TriangulationTraits_3.h index d18e42be86d..ff7448b2e32 100644 --- a/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3TriangulationTraits_3.h +++ b/Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3TriangulationTraits_3.h @@ -80,11 +80,7 @@ of `Kernel::Tetrahedron_3`. typedef unspecified_type Tetrahedron_3; /*! -A predicate object that must provide the function operators - -`Comparison_result operator()(Point_3 p, Point_3 q)`, - -which returns `EQUAL` if the two points are equal and +A predicate object that must provide the function operator `Comparison_result operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`, @@ -96,14 +92,7 @@ in a same line. typedef unspecified_type Compare_xyz_3; /*! -A predicate object that must provide the function operators - -`Orientation operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`, - -which returns `POSITIVE`, if `s` lies on the positive side of -the oriented plane `h` defined by `p`, `q`, and `r`, -returns `NEGATIVE` if `s` lies on the negative side of -`h`, and returns `COPLANAR` if `s` lies on `h` and +A predicate object that must provide the function operator `Orientation operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`, @@ -134,11 +123,7 @@ which constructs a point from a point-offset pair. typedef unspecified_type Construct_point_3; /*! -A constructor object that must provide the function operators - -`Segment_3 operator()(Point_3 p, Point_3 q)`, - -which constructs a segment from two points and +A constructor object that must provide the function operator `Segment_3 operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`, @@ -148,11 +133,7 @@ which constructs a segment from two point-offset pairs. typedef unspecified_type Construct_segment_3; /*! -A constructor object that must provide the function operators - -`Triangle_3 operator()(Point_3 p, Point_3 q, Point_3 r )`, - -which constructs a triangle from three points and +A constructor object that must provide the function operator `Triangle_3 operator()(Point_3 p, Point_3 q, Point_3 r, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r)`, @@ -162,11 +143,7 @@ which constructs a triangle from three point-offset pairs. typedef unspecified_type Construct_triangle_3; /*! -A constructor object that must provide the function operators - -`Tetrahedron_3 operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`, - -which constructs a tetrahedron from four points and +A constructor object that must provide the function operator `Tetrahedron_3 operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`, diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h index 17198cb4c86..8c0edf178d0 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h +++ b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h @@ -733,37 +733,96 @@ public: Vertex_handle nearest_power_vertex(const Bare_point& p, Cell_handle start) const { + CGAL_triangulation_precondition(p.x() < domain().xmax()); + CGAL_triangulation_precondition(p.y() < domain().ymax()); + CGAL_triangulation_precondition(p.z() < domain().zmax()); + CGAL_triangulation_precondition(p.x() >= domain().xmin()); + CGAL_triangulation_precondition(p.y() >= domain().ymin()); + CGAL_triangulation_precondition(p.z() >= domain().zmin()); + if(number_of_vertices() == 0) return Vertex_handle(); + typename Gt::Construct_weighted_point_3 p2wp = + geom_traits().construct_weighted_point_3_object(); + Locate_type lt; int li, lj; + Offset query_offset; + Cell_handle c = locate(p2wp(p), query_offset, lt, li, lj, start); - typename Gt::Construct_weighted_point_3 p2wp = - geom_traits().construct_weighted_point_3_object(); - Cell_handle c = locate(p2wp(p), lt, li, lj, start); - if(lt == Tr_Base::VERTEX) - return c->vertex(li); - const Conflict_tester tester(p2wp(p), this); - Offset o = combine_offsets(Offset(), get_location_offset(tester, c)); +#ifdef CGAL_PERIODIC_DEBUG_NEAREST_POWER_VERTEX + std::cout << "nearest power vertex: " << p << std::endl; + std::cout << "vertices: " << number_of_vertices() << std::endl; + std::cout << "stored vertices: " << this->number_of_stored_vertices() << std::endl; + std::cout << "Locate: " << p << std::endl; + std::cout << "Cell: " << &*c << std::endl; + std::cout << this->point(c, 0) << std::endl; + std::cout << this->point(c, 1) << std::endl; + std::cout << this->point(c, 2) << std::endl; + std::cout << this->point(c, 3) << std::endl; + std::cout << "offset query: " << query_offset << std::endl; + std::cout << "bounded side : " << geom_traits().bounded_side_3_object()( + Tetrahedron(this->point(c, 0).point(), this->point(c, 1).point(), + this->point(c, 2).point(), this->point(c, 3).point()), p) << std::endl; + std::cout << "power side: " << geom_traits().power_side_of_bounded_power_sphere_3_object()( + this->point(c, 0), this->point(c, 1), + this->point(c, 2), this->point(c, 3), p2wp(p)) << std::endl; + std::cout << "power distance: " << geom_traits().compute_power_distance_to_power_sphere_3_object()( + this->point(c, 0), this->point(c, 1), + this->point(c, 2), this->point(c, 3), p2wp(p)) << std::endl; +#endif // - 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(c, p, o); + + // Take the opposite because periodic_locate() returns the offset such that + // cell + offset contains 'p' but here we need to move 'p' + query_offset = this->combine_offsets(Offset(), -query_offset); + + Vertex_handle nearest = nearest_vertex_in_cell(c, p, query_offset); + Offset offset_of_nearest = get_min_dist_offset(p, query_offset, nearest); + +#ifdef CGAL_PERIODIC_DEBUG_NEAREST_POWER_VERTEX + std::cout << "nearest vertex in cell : " << &*nearest << " : " << nearest->point() << std::endl; + std::cout << "offset_of_nearest: " << offset_of_nearest << std::endl; +#endif + std::vector vs; vs.reserve(32); while(true) { Vertex_handle tmp = nearest; - Offset tmp_off = get_min_dist_offset(p, o, tmp); +#ifdef CGAL_PERIODIC_DEBUG_NEAREST_POWER_VERTEX + std::cout << "tmp set to : " << &*nearest << " : " << nearest->point() + << " || offset: " << nearest->offset() << std::endl; +#endif + adjacent_vertices(nearest, std::back_inserter(vs)); - for(typename std::vector::const_iterator vsit = vs.begin(); vsit != vs.end(); ++vsit) - tmp = (compare_distance(p, tmp->point(), (*vsit)->point(), - o, tmp_off, get_min_dist_offset(p, o, *vsit)) - == SMALLER) ? tmp : *vsit; + for(typename std::vector::const_iterator vsit = vs.begin(); + vsit != vs.end(); ++vsit) + { + // Can happen in 27-sheeted triangulations composed of few points + if((*vsit)->point() == nearest->point()) + continue; + + const Offset min_dist_offset = get_min_dist_offset(p, query_offset, *vsit); + if(compare_distance(p, (*vsit)->point(), tmp->point(), + query_offset, min_dist_offset, offset_of_nearest) == SMALLER) + { + tmp = *vsit; + offset_of_nearest = min_dist_offset; +#ifdef CGAL_PERIODIC_DEBUG_NEAREST_POWER_VERTEX + std::cout << " Closer adjacent vertex: " << &*tmp << " : " << tmp->point() + << " || offset " << offset_of_nearest << std::endl; +#endif + } + } + if(tmp == nearest) break; + vs.clear(); nearest = tmp; } @@ -834,7 +893,7 @@ public: return true; } - bool is_Gabriel(const Facet& f)const + bool is_Gabriel(const Facet& f) const { return is_Gabriel(f.first, f.second); } @@ -846,8 +905,20 @@ public: bool is_Gabriel(Vertex_handle v) const { - return nearest_power_vertex( - geom_traits().construct_point_3_object()(v->point()), v->cell()) == v; + 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(); + + const Bare_point& bp = geom_traits().construct_point_3_object()(v->point()); + + Vertex_handle nearest_v = nearest_power_vertex(bp, v->cell()); + + // Need to find the offset such that power distance v->point() + // to nearest_v->point() is minimum + Offset nearest_v_off = get_min_dist_offset_general(bp, nearest_v); + + return (side_of_bounded_orthogonal_sphere( + v->point(), nearest_v->point(), Offset(), nearest_v_off) != CGAL::ON_BOUNDED_SIDE); } Offset get_min_dist_offset(const Bare_point& p, const Offset & o, @@ -878,6 +949,36 @@ public: return combine_offsets(mdo,min_off); } + // In this version, `p` must be in the domain, and we check all possible + // offsets to find the minimum + Offset get_min_dist_offset_general(const Bare_point& p, const Vertex_handle vh) const + { + CGAL_triangulation_precondition(p.x() < domain().xmax()); + CGAL_triangulation_precondition(p.y() < domain().ymax()); + CGAL_triangulation_precondition(p.z() < domain().zmax()); + CGAL_triangulation_precondition(p.x() >= domain().xmin()); + CGAL_triangulation_precondition(p.y() >= domain().ymin()); + CGAL_triangulation_precondition(p.z() >= domain().zmin()); + + Offset min_off = Offset(0,0,0); + + for(int i=-1; i<=1; ++i) { + for(int j=-1; j<=1; ++j) { + for(int k=-1; k<=1; ++k) + { + if(i==0 && j==0 && k==0) + continue; + + Offset loc_off(i, j, k); + if(compare_distance(p, vh->point(), vh->point(), Offset(), min_off, loc_off) == LARGER) + min_off = loc_off; + } + } + } + + return min_off; + } + Vertex_handle nearest_vertex_in_cell(const Cell_handle& c, const Bare_point& p, const Offset & o) const { diff --git a/Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_regular_triangulation_3.cpp b/Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_regular_triangulation_3.cpp index d7555237f02..95600fca004 100644 --- a/Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_regular_triangulation_3.cpp +++ b/Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_regular_triangulation_3.cpp @@ -42,6 +42,7 @@ public: typedef typename P3RT3::Facet Facet; typedef typename P3RT3::Cell Cell; typedef typename P3RT3::Vertex_iterator Vertex_iterator; + typedef typename P3RT3::Unique_vertex_iterator Unique_vertex_iterator; typedef typename P3RT3::Cell_iterator Cell_iterator; typedef typename P3RT3::Segment Segment; typedef typename P3RT3::Triangle Triangle; @@ -77,6 +78,100 @@ public: assert(p3rt3.is_valid()); } + static void test_is_gabriel() + { + std::cout << "--- test is_gabriel" << std::endl; + + P3RT3 p3rt3; + + typename P3RT3::Geom_traits::Power_side_of_bounded_power_sphere_3 + side_of_bounded_power_sphere = + p3rt3.geom_traits().power_side_of_bounded_power_sphere_3_object(); + + Weighted_point_3 t(Point_3(0.5,0.5,0.45), 0.01); + Weighted_point_3 s(Point_3(0.5,0.5,0.49), 0.006); + Weighted_point_3 q(Point_3(0.5,0.5,0.5), 0.015); + Weighted_point_3 p(Point_3(0.95,0.95,0.96), 0.001); + Weighted_point_3 r(Point_3(0.01,0.008,0.01), 0.002); + p3rt3.insert(p); + p3rt3.insert(q); + p3rt3.insert(r); + p3rt3.insert(s); + p3rt3.insert(t); + + assert(p3rt3.is_valid()); + + std::cout << "p3rt3.number_of_vertices() " << p3rt3.number_of_vertices() << std::endl; + assert(p3rt3.number_of_vertices() == 4); + assert(std::distance(p3rt3.unique_vertices_begin(), p3rt3.unique_vertices_end()) == 4); + assert(p3rt3.number_of_stored_vertices() == 108); + + for(Unique_vertex_iterator iter = p3rt3.unique_vertices_begin(), end_iter = p3rt3.unique_vertices_end(); iter != end_iter; ++iter) + { + Vertex_handle vh((Vertex_iterator(iter))); + std::cout << p3rt3.is_Gabriel(vh) << std::endl; + if(p3rt3.is_Gabriel(vh)) + { + for(Unique_vertex_iterator iter2 = p3rt3.unique_vertices_begin(), end_iter2 = p3rt3.unique_vertices_end(); iter2 != end_iter2; ++iter2) + { + Vertex_handle vh2((Vertex_iterator(iter2))); + + if(vh2->point() == vh->point()) + { + assert(p3rt3.is_Gabriel(vh2)); // consistency check + } + else + { + // Check that w/e the offset, the power distance is positive + for(int i = -1; i < 2; ++i) { + for(int j = -1; j < 2; ++j) { + for(int k = -1; k < 2; ++k) + { + const Offset off(i, j, k); +// std::cout << "power distance: " << p3rt3.geom_traits().compute_power_product_3_object()( +// Weighted_point_3(vh->point().point(), vh->point().weight()), +// p3rt3.geom_traits().construct_weighted_point_3_object()(vh2->point(), off)) << std::endl; + if(!(side_of_bounded_power_sphere(vh->point(), vh2->point(), + Offset(), off) != CGAL::ON_BOUNDED_SIDE)) + { + assert(false); + } + } + } + } + } + } + } + else + { + bool found = false; + for(Vertex_iterator iter2 = p3rt3.vertices_begin(), end_iter2 = p3rt3.vertices_end(); iter2 != end_iter2; ++iter2) + { + Vertex_handle vh2 = iter2; + if(vh2->point() == vh->point()) + { + assert(!p3rt3.is_Gabriel(vh2)); // consistency check + } + else + { + for(int i = -1; i < 2; ++i) { + for(int j = -1; j < 2; ++j) { + for(int k = -1; k < 2; ++k) + { + const Offset off = vh->offset() + Offset(i, j, k); + if(!(side_of_bounded_power_sphere(vh->point(), vh2->point(), + vh->offset(), off) != CGAL::ON_BOUNDED_SIDE)) + found = true; + } + } + } + } + } // iter2 + assert(found); // must have found a point that is in the smallest orthogonal power sphere + } // is_gabriel(vh) + } + } + static void test_insert_1 () { std::cout << "--- test_insert_1" << std::endl; @@ -721,6 +816,7 @@ public: static void test () { + test_is_gabriel(); test_find_conflicts(); test_insert_range(800, 7); test_construction_and_insert_range(800, 7); @@ -736,7 +832,6 @@ public: test_insert_two_points_with_the_same_position(); test_remove(); test_27_to_1_sheeted_covering(); - ////// Iso_cuboid unitaire -> 0 <= weight < 0.015625 test_insert_rnd_as_delaunay(100, 0.); test_insert_rnd_as_delaunay(100, 0.01); } diff --git a/Triangulation_3/doc/Triangulation_3/Concepts/DelaunayTriangulationTraits_3.h b/Triangulation_3/doc/Triangulation_3/Concepts/DelaunayTriangulationTraits_3.h index a92e75955bc..f7812c32d9f 100644 --- a/Triangulation_3/doc/Triangulation_3/Concepts/DelaunayTriangulationTraits_3.h +++ b/Triangulation_3/doc/Triangulation_3/Concepts/DelaunayTriangulationTraits_3.h @@ -81,6 +81,29 @@ It is only needed when using the `Fast_location` policy or the typedef unspecified_type Compare_distance_3; +/// @} + +/// \name +/// When `is_Gabriel` functions are used, the traits class must +/// in addition provide the following predicate object: +/// @{ + +/*! +A predicate object that must provide the function operators + +`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 t)`, + +which returns the position of the point `t` relative to the sphere +that has `pq` as its diameter, + +`Bounded_side operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 t)`, + +which returns the position of the point `t` relative to the sphere +passing through `p, q`, and `r` and whose center is in the +plane defined by these three points. +*/ +typedef unspecified_type Side_of_bounded_sphere_3; + /// @} /*! \name @@ -166,8 +189,6 @@ When using the `Fast_location` policy or the `CGAL::Delaunay_triangulation_3::ne */ Compare_distance_3 compare_distance_3_object(); - - /// @} /*! \name diff --git a/Triangulation_3/doc/Triangulation_3/Concepts/RegularTriangulationTraits_3.h b/Triangulation_3/doc/Triangulation_3/Concepts/RegularTriangulationTraits_3.h index bf16e8b37ef..dd6c787c8a5 100644 --- a/Triangulation_3/doc/Triangulation_3/Concepts/RegularTriangulationTraits_3.h +++ b/Triangulation_3/doc/Triangulation_3/Concepts/RegularTriangulationTraits_3.h @@ -210,6 +210,34 @@ typedef unspecified_type Construct_ray_3; /// @} +/// \name +/// When `is_Gabriel` functions are used, the traits class must +/// in addition provide the following predicate object: +/// @{ + +/*! +A predicate object that must provide the function operators + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 t)`, + +which returns the sign of the power test of `t` with respect to the smallest +sphere orthogonal to `p` (which is the sphere with center `p` and squared +radius `-w_p` with `w_p` the weight of `p`), + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 t)`, + +which returns the sign of the power test of `t` with respect to the smallest +sphere orthogonal to `p` and `q`, + +`Bounded_side operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 t)`, + +which returns the sign of the power test of `t` with respect to the smallest +sphere orthogonal to `p`, `q`, and `r`. +*/ +typedef unspecified_type Power_side_of_bounded_power_sphere_3; + +/// @} + /// \name Operations /// @{ diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index d2a06431086..7287944d586 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -2180,8 +2180,16 @@ namespace CGAL { Regular_triangulation_3:: is_Gabriel(Vertex_handle v) const { - return nearest_power_vertex( - geom_traits().construct_point_3_object()(v->point()), v->cell()) == v; + 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(); + + Vertex_handle nearest_v = + nearest_power_vertex(geom_traits().construct_point_3_object()(v->point()), + v->cell()); + + return (side_of_bounded_orthogonal_sphere(v->point(), nearest_v->point()) + != CGAL::ON_BOUNDED_SIDE); } // Returns