diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h index 5f7818d43e6..fc93ecc8c44 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h @@ -540,38 +540,28 @@ public: template Bounded_side operator()(const Point_2 p, const Face_handle fh, const Offset o) const { - cout << "o = " << o << ", off[0] = " << fh->offset(0) << ", off[1] = " << fh->offset(1) << ", off[2] = " << fh->offset(2) << endl; + + //cout << "Checking face " << fh->get_number() << " with offset " << o << endl; + Point_2 p1 = o.append(fh->offset(0)).apply(fh->vertex(0)->point()); - Point_2 p2 = o.append(fh->offset(2)).apply(fh->vertex(1)->point()); - Point_2 p3 = o.append(fh->offset(1)).apply(fh->vertex(2)->point()); + Point_2 p2 = o.append(fh->offset(1)).apply(fh->vertex(1)->point()); + Point_2 p3 = o.append(fh->offset(2)).apply(fh->vertex(2)->point()); - Construct_segment_2 bld; - Segment_2 s1, s2, s3; + // Construct_segment_2 bld; + // Segment_2 s1, s2, s3; - s1 = bld(p2, p3); - s2 = bld(p3, p1); - s3 = bld(p1, p2); + // s1 = bld(p2, p3); + // s2 = bld(p3, p1); + // s3 = bld(p1, p2); - Bounded_side cs1 = side_of_segment_2(p1, s1); - Bounded_side cp1 = side_of_segment_2(p, s1); + Bounded_side cs1 = side_of_segment_2(p1, p2, p3); + Bounded_side cp1 = side_of_segment_2(p, p2, p3); - cout << (cs1 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << (cp1 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << " "; + Bounded_side cs2 = side_of_segment_2(p2, p3, p1); + Bounded_side cp2 = side_of_segment_2(p, p3, p1); - Bounded_side cs2 = side_of_segment_2(p2, s2); - Bounded_side cp2 = side_of_segment_2(p, s2); - - cout << (cs2 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << (cp2 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << " "; - - Bounded_side cs3 = side_of_segment_2(p3, s3); - Bounded_side cp3 = side_of_segment_2(p, s3); - - cout << (cs3 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << (cp3 == ON_BOUNDED_SIDE ? "+" : "-"); - cout << endl; + Bounded_side cs3 = side_of_segment_2(p3, p1, p2); + Bounded_side cp3 = side_of_segment_2(p, p1, p2); if (cs1 != cp1 || cs2 != cp2 || cs3 != cp3) { return ON_UNBOUNDED_SIDE; @@ -593,21 +583,21 @@ public: Point_2 p2 = fh->offset(2).apply(fh->vertex(1)->point()); Point_2 p3 = fh->offset(1).apply(fh->vertex(2)->point()); - Construct_segment_2 bld; - Segment_2 s1, s2, s3; + // Construct_segment_2 bld; + // Segment_2 s1, s2, s3; - s1 = bld(p2, p3); - s2 = bld(p3, p1); - s3 = bld(p1, p2); + // s1 = bld(p2, p3); + // s2 = bld(p3, p1); + // s3 = bld(p1, p2); - Bounded_side cs1 = side_of_segment_2(p1, s1); - Bounded_side cp1 = side_of_segment_2(p, s1); + Bounded_side cs1 = side_of_segment_2(p1, p2, p3); + Bounded_side cp1 = side_of_segment_2(p, p2, p3); - Bounded_side cs2 = side_of_segment_2(p2, s2); - Bounded_side cp2 = side_of_segment_2(p, s2); + Bounded_side cs2 = side_of_segment_2(p2, p3, p1); + Bounded_side cp2 = side_of_segment_2(p, p3, p1); - Bounded_side cs3 = side_of_segment_2(p3, s3); - Bounded_side cp3 = side_of_segment_2(p, s3); + Bounded_side cs3 = side_of_segment_2(p3, p1, p2); + Bounded_side cp3 = side_of_segment_2(p, p1, p2); if (cs1 != cp1 || cs2 != cp2 || cs3 != cp3) { return ON_UNBOUNDED_SIDE; @@ -622,23 +612,37 @@ public: } -// if ( Circle_2* c_pq = boost::get(&bis_pq) ) - private: - Bounded_side side_of_segment_2(const Point_2 p, const Segment_2 s) const { - if (const Circular_arc_2* c = boost::get(&s)) { - return c->supporting_circle().bounded_side(p); + Bounded_side side_of_segment_2(const Point_2 query, const Point_2 p, const Point_2 q) const { + + Point_2 o(0, 0); + + // Invert p or q through the unit circle. + // The inversion depends on the distance from the origin, so to increase + // numerical stability we choose to invert the point further from (0,0). + Point_2 inv; + FT dp = squared_distance(o, p), dq = squared_distance(o, q); + if (dq < dp) { + inv = Point_2( p.x()/dp, p.y()/dp ); } else { - const Line_arc_2* la = boost::get(&s); - Euclidean_line_2 l = la->supporting_line(); - if (l.has_on(p)) { + inv = Point_2( q.x()/dq, q.y()/dq ); + } + + // If iq is on the line defined by p and q, we need to work with a line instead of a circle + if (orientation(p, q, inv) == COLLINEAR) { + Orientation oquery = orientation(p, q, query); + if (oquery == COLLINEAR) { return ON_BOUNDARY; - } else if (l.has_on_positive_side(p)) { - return ON_BOUNDED_SIDE; + } else if (oquery == LEFT_TURN) { + return ON_BOUNDED_SIDE; // this is just a convention } else { return ON_UNBOUNDED_SIDE; } + } else { // this means that we work in the circle + Circle_2 c(p, q, inv); + return c.bounded_side(query); } + } }; diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h index b1c1195be6b..7ce5d0bdf60 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h @@ -1033,27 +1033,35 @@ hyperbolic_locate(const Point& p, Offset& lo, Locate_type& lt, int& li, int& lj, typedef typename GT::Side_of_hyperbolic_face_2 Side_of_hyperbolic_face_2; Side_of_hyperbolic_face_2 sf; - cout << "Euclidean location says face " << lf->get_number() << " with offset " << lo << endl; + cout << "Euclidean locate says face " << lf->get_number() << " with offset " << lo; Bounded_side bs = sf(p, lf, lo); if (bs == ON_BOUNDED_SIDE) { - cout << "Location was correct! The point is inside the face!" << endl; + // Nothing to do here, the face is really in lf. } else if (bs == ON_BOUNDARY) { - cout << "Not clear if correct! The point is on the boundary of the face!" << endl; + // here the point is on the boundary of the face, so it is on an edge or on a vertex. Still, it's acceptable to say it is in lf. } else { - Bounded_side bs1 = sf(p, lf->neighbor(0), lo); - Bounded_side bs2 = sf(p, lf->neighbor(1), lo); - Bounded_side bs3 = sf(p, lf->neighbor(2), lo); - cout << "The location was not correct! The point is found inside "; - if (bs1 == ON_BOUNDED_SIDE) - cout << " neighbor 0! "; - if (bs2 == ON_BOUNDED_SIDE) - cout << " neighbor 1! "; - if (bs3 == ON_BOUNDED_SIDE) - cout << " neighbor 2! "; - cout << endl; + // Here we have to find the face containing the point, it's one of the neighbors of lf. + Bounded_side bs1 = sf(p, lf->neighbor(0), lo.append(lf->neighbor_offset(0))); + if (bs1 != ON_UNBOUNDED_SIDE) { + lo = lo.append(lf->neighbor_offset(0)); + lf = lf->neighbor(0); + } else { + Bounded_side bs2 = sf(p, lf->neighbor(1), lo.append(lf->neighbor_offset(1))); + if (bs2 != ON_UNBOUNDED_SIDE) { + lo = lo.append(lf->neighbor_offset(1)); + lf = lf->neighbor(1); + } else { + Bounded_side bs3 = sf(p, lf->neighbor(2), lo.append(lf->neighbor_offset(2))); + CGAL_triangulation_assertion(bs3 != ON_UNBOUNDED_SIDE); + lo = lo.append(lf->neighbor_offset(2)); + lf = lf->neighbor(2); + } + } } + cout << "; hyperbolic locate says face " << lf->get_number() << " with offset " << lo << endl; + // Orientation o1 = orientation(f->vertex(0)->point(), f->vertex(1)->point(), p, // f->offset(0), f->offset(1), Offset()); //bu_encode_offset());