From 78db95fcb9d352eece37f2290fb74dd2bd25a8a7 Mon Sep 17 00:00:00 2001 From: Menelaos Karavelas Date: Wed, 23 Jul 2003 03:52:48 +0000 Subject: [PATCH] added the first version for the files in the predicates subdir --- .../Segment_Voronoi_diagram_predicates_C2.h | 1759 +++++++++++++++++ .../Segment_Voronoi_diagram_predicates_ftC2.h | 725 +++++++ .../include/CGAL/predicates/Square_root_1.h | 327 +++ .../include/CGAL/predicates/Square_root_2.h | 379 ++++ 4 files changed, 3190 insertions(+) create mode 100644 Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_C2.h create mode 100644 Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_ftC2.h create mode 100644 Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_1.h create mode 100644 Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_2.h diff --git a/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_C2.h b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_C2.h new file mode 100644 index 00000000000..09c02559def --- /dev/null +++ b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_C2.h @@ -0,0 +1,1759 @@ +#ifndef CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_2_H +#define CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_2_H + + +#include +#include + +CGAL_BEGIN_NAMESPACE + +#define PRED_PRINT 0 + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + +template< class K > +class Svd_are_same_points_2 +{ +public: + typedef typename K::Point_2 Point_2; + + typedef typename K::Compare_x_2 compare_x_2; + typedef typename K::Compare_y_2 compare_y_2; +public: + + inline + bool operator()(const Point_2& p, const Point_2& q) const + { + if ( compare_x_2()(p, q) != EQUAL ) { return false; } + Comparison_result res = compare_y_2()(p, q); + return res == EQUAL; + } +}; + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + +template< class K > +class Svd_is_endpoint_of_segment_2 +{ +public: + typedef typename K::Point_2 Point_2; + typedef typename K::Segment_2 Segment_2; + typedef typename CGAL::Svd_are_same_points_2 Are_same_points_2; + + typedef typename K::Compare_x_2 compare_x_2; + typedef typename K::Compare_y_2 compare_y_2; +public: + + inline + bool operator()(const Point_2& p, const Segment_2& s) const + { + Are_same_points_2 are_same_points; + return ( are_same_points(p, s.source()) || + are_same_points(p, s.target()) ); + } +}; + + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + + +template +class Svd_incircle_2 +{ +private: + typedef typename R::Point_2 Point; + typedef typename R::Segment_2 Segment; + typedef typename R::Site_2 Site; + typedef CGAL::Svd_voronoi_vertex_2 Voronoi_vertex; + + typedef typename R::FT FT; + typedef typename R::RT RT; + +private: + Sign incircle(const Site& p, const Site& q, + const Point& t) const + { + if ( p.is_point() && q.is_point() ) { + Orientation o = orientation(p.point(), q.point(), t); + + if ( o != COLLINEAR ) { +#if PRED_PRINT + std::cout << ((o == LEFT_TURN) ? POSITIVE : NEGATIVE) << std::endl; +#endif + return (o == LEFT_TURN) ? POSITIVE : NEGATIVE; + } + + + RT dtpx = p.point().x() - t.x(); + RT dtpy = p.point().y() - t.y(); + RT dtqx = q.point().x() - t.x(); + RT dtqy = q.point().y() - t.y(); + + Sign s = sign_of_determinant2x2(dtpx, dtpy, -dtqy, dtqx); + + CGAL_assertion( s != ZERO ); +#if PRED_PRINT + std::cout << s << std::endl; +#endif + return s; + } + + CGAL_assertion( p.is_point() || q.is_point() ); + + Orientation o; + if ( p.is_point() && q.is_segment() ) { + Point pq = (p.point() == q.source()) ? q.target() : q.source(); + o = orientation(p.point(), pq, t); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + } else { // p is a segment and q is a point + Point pp = (q.point() == p.source()) ? p.target() : p.source(); + o = orientation(pp, q.point(), t); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + } + return ( o == RIGHT_TURN ) ? NEGATIVE : POSITIVE; + } + + //----------------------------------------------------------------------- + + + Sign incircle(const Point& p, const Point& q, + const Segment& t) const + { + if ( (p == t.source() || p == t.target()) && + (q == t.source() || q == t.target()) ) { + // if t is the segment joining p and q then t must be a vertex + // on the convex hull +#if PRED_PRINT + std::cout << NEGATIVE << std::endl; +#endif + return NEGATIVE; + } else if ( p == t.source() || p == t.target() ) { + // p is an endpoint of t + // in this case the p,q,oo vertex is destroyed only if the + // other endpoint of t is beyond + Point pt = (p == t.source()) ? t.target() : t.source(); + Orientation o = orientation(p, q, pt); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + return (o == RIGHT_TURN) ? NEGATIVE : POSITIVE; + } else if ( q == t.source() || q == t.target() ) { + Point pt = (q == t.source()) ? t.target() : t.source(); + Orientation o = orientation(p, q, pt); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + return (o == RIGHT_TURN) ? NEGATIVE : POSITIVE; + } else { + // maybe here I should immediately return POSITIVE; + // since we insert endpoints of segments first, p and q cannot + // be consecutive points on the convex hull if one of the + // endpoints of t is to the right of the line pq. + Orientation o1 = orientation(p, q, t.source()); + Orientation o2 = orientation(p, q, t.target()); + + if ( o1 == RIGHT_TURN || o2 == RIGHT_TURN ) { +#if PRED_PRINT + std::cout << NEGATIVE << std::endl; +#endif + return NEGATIVE; + } +#if PRED_PRINT + std::cout << POSITIVE << std::endl; +#endif + return POSITIVE; + } + } + + + Sign incircle(const Segment& p, const Point& q, + const Segment& t) const + { + if ( q == t.source() && q == t.target() ) { + Point pp = (q == p.source()) ? p.target() : p.source(); + Point pt = (q == t.source()) ? t.target() : t.source(); + + Orientation o = orientation(pp, q, pt); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + return (o == RIGHT_TURN) ? NEGATIVE : POSITIVE; + } else { +#if PRED_PRINT + std::cout << POSITIVE << std::endl; +#endif + return POSITIVE; + } + } + + + Sign incircle(const Point& p, const Segment& q, + const Segment& t) const + { + if ( p == t.source() || p == t.target() ) { + Point pq = (p == q.source()) ? q.target() : q.source(); + Point pt = (p == t.source()) ? t.target() : t.source(); + + Orientation o = orientation(p, pq, pt); + +#if PRED_PRINT + std::cout << ((o == RIGHT_TURN) ? NEGATIVE : POSITIVE) << std::endl; +#endif + return (o == RIGHT_TURN) ? NEGATIVE : POSITIVE; + } else { + // if p is not an endpoint of t, then either p and q should + // not be on the convex hull or t does not affect the vertex + // of p and q. +#if PRED_PRINT + std::cout << POSITIVE << std::endl; +#endif + return POSITIVE; + } + } + + + Sign incircle(const Site& p, const Site& q, + const Segment& t) const + { + if ( p.is_point() && q.is_point() ) { + return incircle(p.point(), q.point(), t); + } else if ( p.is_point() && q.is_segment() ) { + return incircle(p.point(), q.segment(), t); + } else { // p is a segment and q is a point + return incircle(p.segment(), q.point(), t); + } + } + + +public: + + Sign operator()(const Site& p, const Site& q, + const Site& r, const Site& t) const + { +#if PRED_PRINT + std::cout << "Incircle(p,q,r,t): " << std::endl; + std::cout << p << std::endl; + std::cout << q << std::endl; + std::cout << r << std::endl; + std::cout << t << std::endl; +#endif + + Voronoi_vertex v(p, q, r); + +#if PRED_PRINT + std::cout << v.incircle(t, Method_tag()) << std::endl; +#endif + + return v.incircle(t, Method_tag()); + } + + + + + inline Sign operator()(const Site& p, const Site& q, + const Site& t) const + { +#if PRED_PRINT + std::cout << "Incircle(p,q,t): " << std::endl; + std::cout << p << std::endl; + std::cout << q << std::endl; + std::cout << t << std::endl; +#endif + + CGAL_assertion( !(p.is_segment() && q.is_segment()) ); + + if ( p.is_point() && q.is_segment() ) { + // p must be an endpoint of q + CGAL_assertion( (p.point() == q.source()) || + (p.point() == q.target()) ); + } else if ( p.is_segment() && q.is_point() ) { + // q must be an endpoint of p + CGAL_assertion( (p.source() == q.point()) || + (p.target() == q.point()) ); + } + + if ( t.is_point() ) { + return incircle(p, q, t.point()); + } + + return incircle(p, q, t.segment()); + } + + +}; + + + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + + + +template +class Svd_finite_edge_interior_2 +{ +public: + typedef typename K::Point_2 Point; + typedef typename K::Segment_2 Segment; + typedef typename K::Line_2 Line; + typedef CGAL::Svd_voronoi_vertex_2 Voronoi_vertex; + typedef typename K::Site_2 Site; + typedef typename K::FT FT; + typedef typename K::RT RT; + + typedef typename Voronoi_vertex::vertex_t vertex_t; + typedef typename Voronoi_vertex::PPP_Type PPP_Type; + typedef typename Voronoi_vertex::PPS_Type PPS_Type; + typedef typename Voronoi_vertex::PSS_Type PSS_Type; + typedef typename Voronoi_vertex::SSS_Type SSS_Type; + + typedef typename Voronoi_vertex::Sqrt_1 Sqrt_1; + typedef typename Voronoi_vertex::Sqrt_2 Sqrt_2; + typedef typename Voronoi_vertex::Sqrt_3 Sqrt_3; + +#if 0 + class Line + { + private: + RT a_, b_, c_; + + public: + Line() : a_(1), b_(0), c_(0) {} + Line(const RT& a, const RT& b, const RT& c) + : a_(a), b_(b), c_(c) {} + + RT a() const { return a_; } + RT b() const { return b_; } + RT c() const { return c_; } + + }; +#endif + + + struct Homogeneous_point_2 + { + RT hx_, hy_, hw_; + + Homogeneous_point_2() : hx_(0), hy_(0), hw_(1) {} + Homogeneous_point_2(const RT& hx, const RT& hy, const RT& hw) + : hx_(hx), hy_(hy), hw_(hw) + { + CGAL_precondition( !(CGAL_NTS is_zero(hw_)) ); + } + + Homogeneous_point_2(const Point& p) + : hx_(p.x()), hy_(p.y()), hw_(1) {} + + Homogeneous_point_2(const Homogeneous_point_2& other) + : hx_(other.hx_), hy_(other.hy_), hw_(other.hw_) {} + + RT hx() const { return hx_; } + RT hy() const { return hy_; } + RT hw() const { return hw_; } + + FT x() const { return hx_ / hw_; } + FT y() const { return hy_ / hw_; } + }; + +public: + static + Homogeneous_point_2 + projection_on_line(const Line& l, const Point& p) + { + RT ab = l.a() * l.b(); + + RT hx = CGAL_NTS square(l.b()) * p.x() + - ab * p.y() - l.a() * l.c(); + RT hy = CGAL_NTS square(l.a()) * p.y() + - ab * p.x() - l.b() * l.c(); + RT hw = CGAL_NTS square(l.a()) + CGAL_NTS square(l.b()); + + return Homogeneous_point_2(hx, hy, hw); + } + + + static + Homogeneous_point_2 + midpoint(const Point& p1, const Point& p2) + { + RT hx = p1.x() + p2.x(); + RT hy = p1.y() + p2.y(); + RT hw = RT(2); + + return Homogeneous_point_2(hx, hy, hw); + } + + static + Homogeneous_point_2 + midpoint(const Homogeneous_point_2& p1, + const Homogeneous_point_2& p2) + { + RT hx = p1.hx() * p2.hw() + p2.hx() * p1.hw(); + RT hy = p1.hy() * p2.hw() + p2.hy() * p1.hw(); + RT hw = RT(2) * p1.hw() * p2.hw(); + + return Homogeneous_point_2(hx, hy, hw); + } + + + static + Line compute_supporting_line(const Segment& s) + { +#if 1 + RT a, b, c; + Voronoi_vertex::compute_supporting_line(s, a, b, c); + return Line(a, b, c); +#else + return Line(s); +#endif + } + + static + Comparison_result + compare_squared_distances_to_line(const Line& l, const Point& p, + const Point& q) + { +#if 1 + RT d2_lp = CGAL_NTS square(l.a() * p.x() + l.b() * p.y() + l.c()); + RT d2_lq = CGAL_NTS square(l.a() * q.x() + l.b() * q.y() + l.c()); + + return CGAL_NTS compare(d2_lp, d2_lq); +#else + FT d2_lp = squared_distance_to_line(p, l); + FT d2_lq = squared_distance_to_line(q, l); + + return CGAL_NTS compare(d2_lp, d2_lq); +#endif + } + + static + Comparison_result + compare_squared_distances_to_lines(const Point& p, const Line& l1, + const Line& l2) + { +#if 1 + RT d2_l1 = CGAL_NTS square(l1.a() * p.x() + l1.b() * p.y() + l1.c()); + RT d2_l2 = CGAL_NTS square(l2.a() * p.x() + l2.b() * p.y() + l2.c()); + + RT n1 = CGAL_NTS square(l1.a()) + CGAL_NTS square(l1.b()); + RT n2 = CGAL_NTS square(l2.a()) + CGAL_NTS square(l2.b()); + + return CGAL_NTS compare(d2_l1 * n2, d2_l2 * n1); +#else + + FT d2_p_from_l1 = squared_distance_to_line(p, l1); + FT d2_p_from_l2 = squared_distance_to_line(p, l2); + + return CGAL_NTS compare(d2_p_from_l1, d2_p_from_l2); +#endif + } + + static + Oriented_side + oriented_side_of_line(const Line& l, const Point& p) + { +#if 1 + Sign s = CGAL_NTS sign(l.a() * p.x() + l.b() * p.y() + l.c()); + if ( s == ZERO ) { return ON_ORIENTED_BOUNDARY; } + return ( s == POSITIVE ) ? ON_POSITIVE_SIDE : ON_NEGATIVE_SIDE; +#else + return l.oriented_side(p); +#endif + } + + + static + Oriented_side + oriented_side_of_line(const Line& l, const Homogeneous_point_2& p) + { + Sign s1 = + CGAL_NTS sign(l.a() * p.hx() + l.b() * p.hy() + l.c() * p.hw()); + Sign s_hw = CGAL_NTS sign(p.hw()); + + Sign s = Sign(s1 * s_hw); + + if ( s == ZERO ) { return ON_ORIENTED_BOUNDARY; } + return ( s == POSITIVE ) ? ON_POSITIVE_SIDE : ON_NEGATIVE_SIDE; + } + + + + static + Line compute_perpendicular(const Line& l, const Point& p) + { +#if 1 + RT a, b, c; + a = -l.b(); + b = l.a(); + c = l.b() * p.x() - l.a() * p.y(); + return Line(a, b, c); +#else + return l.perpendicular(p); +#endif + } + + static + Line opposite_line(const Line& l) + { +#if 1 + return Line(-l.a(), -l.b(), -l.c()); +#else + return l.opposite(); +#endif + } + +private: + static + FT squared_distance_to_line(const Point& p, const Line& l) + { + Point proj = l.projection(p); + Segment s(p, proj); + return s.squared_length(); + } + + static + FT + squared_distance_to_line(const Point& p, const Segment& s) + { + return squared_distance_to_line(p, Line(s)); + } + + +private: + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- +#if 1 + bool + is_interior_in_conflict_both(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + if ( p.is_point() && q.is_point() ) { + return + is_interior_in_conflict_both(p.point(), q.point(), + r, s, t, tag); + } else if ( p.is_segment() && q.is_segment() ) { + return + is_interior_in_conflict_both(p.segment(), q.segment(), + r, s, t, tag); + } else if ( p.is_point() && q.is_segment() ) { + return + is_interior_in_conflict_both(p.point(), q.segment(), + r, s, t, tag); + } else { // p is a segment and q is a point + return + is_interior_in_conflict_both(p.segment(), q.point(), + r, s, t, tag); + } + } +#endif + //-------------------------------------------------------------------- + + bool + is_interior_in_conflict_both(const Point& p, const Point& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + if ( t.is_point() ) { return true; } + + Line lt = compute_supporting_line(t.segment()); + + // MK:: THIS IS REALLY IMPORTANT BECAUSE IT DEPENDS ON HOW LINES + // ARE IMPLEMENTED ******************************************* + // change from passing a line to passing a triple of numbers + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + Oriented_side op, oq; + + if ( p == t.segment().source() || + p == t.segment().target() ) { + op = ON_ORIENTED_BOUNDARY; + } else { + op = oriented_side_of_line(lt, p); + } + + if ( q == t.segment().source() || + q == t.segment().target() ) { + oq = ON_ORIENTED_BOUNDARY; + } else { + oq = oriented_side_of_line(lt, q); + } + + + if ((op == ON_POSITIVE_SIDE && oq == ON_NEGATIVE_SIDE) || + (op == ON_NEGATIVE_SIDE && oq == ON_POSITIVE_SIDE) || + (op == ON_ORIENTED_BOUNDARY || oq == ON_ORIENTED_BOUNDARY)) { + return true; + } + + Comparison_result res = + compare_squared_distances_to_line(lt, p, q); + + if ( res == EQUAL ) { return true; } + + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + + + Line lperp; + if ( res == SMALLER ) { + // p is closer to lt than q + lperp = compute_perpendicular(lt, p); + } else { + // q is closer to lt than p + lperp = compute_perpendicular(lt, q); + } + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + return ( opqr == oqps ); + } + + //-------------------------------------------------------------------- + + bool + is_interior_in_conflict_both(const Segment& p, const Segment& q, + const Site& r, const Site& s, + const Site& t, Method_tag) const + { + return true; + } + + //-------------------------------------------------------------------- + + bool + is_interior_in_conflict_both(const Point& p, const Segment& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + if ( p == q.source() || p == q.target() ) { + return false; + } + + if ( t.is_point() ) { + return is_interior_in_conflict_both(p, q, r, s, t.point(), tag); + } + return is_interior_in_conflict_both(p, q, r, s, t.segment(), tag); + } + + bool + is_interior_in_conflict_both(const Point& p, const Segment& q, + const Site& r, const Site& s, + const Point& t, Method_tag tag) const + { + Line lq = compute_supporting_line(q); + + Comparison_result res = + compare_squared_distances_to_line(lq, p, t); + + if ( res != SMALLER ) { return true; } + + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + + Line lperp = compute_perpendicular(lq, p); + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + return (opqr == oqps); + } + + bool + is_interior_in_conflict_both(const Point& p, const Segment& q, + const Site& r, const Site& s, + const Segment& t, Method_tag tag) const + { + Line lt = compute_supporting_line(t); + Line lq = compute_supporting_line(q); + + if ( oriented_side_of_line(lq, p) == ON_NEGATIVE_SIDE ) { + lq = opposite_line(lq); + } + + if ( p == t.source() || p == t.target() ) { + Line lqperp = compute_perpendicular(lq, p); + + + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + + Oriented_side opqr = vpqr.oriented_side(lqperp, tag); + Oriented_side oqps = vqps.oriented_side(lqperp, tag); + + bool on_different_parabola_arcs = + ((opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE) || + (opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE)); + + if ( !on_different_parabola_arcs ) { return true; } + + Point t1; + if ( p == t.source() ) { + t1 = t.target(); + } else { + t1 = t.source(); + } + + Oriented_side o_p = oriented_side_of_line(lq, p); + Oriented_side o_t1 = oriented_side_of_line(lq, t1); + + if ( (o_p == ON_POSITIVE_SIDE && o_t1 == ON_NEGATIVE_SIDE) || + (o_p == ON_NEGATIVE_SIDE && o_t1 == ON_POSITIVE_SIDE) ) { + return true; + } + + + Comparison_result res = + compare_squared_distances_to_line(lq, p, t1); + + return ( res == LARGER ); + } + + if ( oriented_side_of_line(lt, p) == ON_NEGATIVE_SIDE ) { + lt = opposite_line(lt); + } + + + Comparison_result res = + CGAL_NTS compare(lt.a() * lq.b(), lt.b() * lq.a()); + bool are_parallel = (res == EQUAL); + + if ( are_parallel ) { + Sign sgn = CGAL_NTS sign(lt.a() * lq.a() + lt.b() * lq.b()); + bool have_opposite_directions = (sgn == NEGATIVE); + if ( have_opposite_directions ) { lq = opposite_line(lq); } + + if ( oriented_side_of_line(lq, p) == oriented_side_of_line(lt, p) ) { + return true; + } + + if ( have_opposite_directions ) { + lq = opposite_line(lq); + } + } + + Line l = compute_perpendicular(lt, p); + + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + + Oriented_side o_l_pqr = vpqr.oriented_side(l, tag); + Oriented_side o_l_qps = vqps.oriented_side(l, tag); + if ( o_l_pqr == ON_POSITIVE_SIDE && + o_l_qps == ON_NEGATIVE_SIDE ) { return false; } + if ( o_l_pqr == ON_NEGATIVE_SIDE && + o_l_qps == ON_POSITIVE_SIDE ) { return true; } + + //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //>>>>>>>>>> HERE I NEED TO CHECK THE BOUNDARY CASES <<<<<< + //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + Line lqperp = compute_perpendicular(lq, p); + + Oriented_side opqr = vpqr.oriented_side(lqperp, tag); + Oriented_side oqps = vqps.oriented_side(lqperp, tag); + + bool on_different_parabola_arcs = + ((opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE) || + (opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE)); + + if ( !on_different_parabola_arcs ) { return true; } + + + Homogeneous_point_2 pv = projection_on_line(lq, p); + Homogeneous_point_2 hp(p); + + pv = midpoint(pv, hp); + + Oriented_side o_l_pv = oriented_side_of_line(l, pv); + + CGAL_assertion( o_l_pv != ON_ORIENTED_BOUNDARY ); + + CGAL_assertion( o_l_pqr != ON_ORIENTED_BOUNDARY || + o_l_qps != ON_ORIENTED_BOUNDARY ); + + if ( o_l_pqr == ON_ORIENTED_BOUNDARY ) { + return ( o_l_qps == o_l_pv ); + } else { + return ( o_l_pqr == o_l_pv ); + } + + } + + //-------------------------------------------------------------------- + + bool + is_interior_in_conflict_both(const Segment& p, const Point& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + return is_interior_in_conflict_both(q, p, s, r, t, tag); + } + + + //-------------------------------------------------------------------- + +#if 0 + bool is_interior_in_conflict_both(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Ring_tag tag) + const + { + // THIS NEEDS TO BE CORRECTED + return + is_interior_in_conflict_both(p, q, r, s, t, Sqrt_field_tag()); + } + + bool is_interior_in_conflict_both(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Sqrt_field_tag tag) + const + { + // checks if interior of voronoi edge is in conflict if both extrema + // of the voronoi edge are in conclift + // return true if interior is in conflict; false otherwise + + if ( p.is_point() && q.is_point() ) { + // both p and q are points + if ( t.is_point() ) { return true; } + + Line lt = compute_supporting_line(t.segment()); + + Oriented_side op, oq; + + if ( p.point() == t.segment().source() || + p.point() == t.segment().target() ) { + op = ON_ORIENTED_BOUNDARY; + } else { + op = oriented_side_of_line(lt, p.point()); + } + + if ( q.point() == t.segment().source() || + q.point() == t.segment().target() ) { + oq = ON_ORIENTED_BOUNDARY; + } else { + oq = oriented_side_of_line(lt, q.point()); + } + + + if ((op == ON_POSITIVE_SIDE && oq == ON_NEGATIVE_SIDE) || + (op == ON_NEGATIVE_SIDE && oq == ON_POSITIVE_SIDE) || + (op == ON_ORIENTED_BOUNDARY || oq == ON_ORIENTED_BOUNDARY)) { + return true; + } + + Comparison_result res = + compare_squared_distances_to_line(lt, p.point(), q.point()); + // FT d2_p_from_lt = squared_distance_to_line(p.point(), lt); + // FT d2_q_from_lt = squared_distance_to_line(q.point(), lt); + + // res = CGAL_NTS compare(d2_p_from_lt, d2_q_from_lt); + + if ( res == EQUAL ) { return true; } + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + Line lperp; + if ( res == SMALLER ) { + // p is closer to lt than q + lperp = compute_perpendicular(lt, p.point()); + } else { + // q is closer to lt than p + lperp = compute_perpendicular(lt, q.point()); + } + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + return ( opqr == oqps ); + } else if ( p.is_segment() && q.is_segment() ) { + // both p and q are segments + return true; + + } else if ( p.is_point() && q.is_segment() ) { + // p is a point q is a segment + if ( p.point() == q.segment().source() || + p.point() == q.segment().target() ) { + return false; + } + if ( t.is_point() ) { + Line lq = compute_supporting_line(q.segment()); + + // FT d2_p_from_lq = squared_distance_to_line(p.point(), lq); + // FT d2_t_from_lq = squared_distance_to_line(t.point(), lq); + + // Comparison_result res = + // CGAL_NTS compare(d2_p_from_lq, d2_t_from_lq); + + Comparison_result res = + compare_squared_distances_to_line(lq, p.point(), t.point()); + + if ( res != SMALLER ) { return true; } + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + // Line lperp = lq.perpendicular(p.point()); + Line lperp = compute_perpendicular(lq, p.point()); + + // Oriented_side opqr = lperp.oriented_side(vpqr.point()); + // Oriented_side oqps = lperp.oriented_side(vqps.point()); + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + return (opqr == oqps); + } else { + // t is a segment + Line lt = compute_supporting_line(t.segment()); + Line lq = compute_supporting_line(q.segment()); + + // if ( lq.has_on_negative_side(p.point()) ) { + if ( oriented_side_of_line(lq, p.point()) == + ON_NEGATIVE_SIDE ) { + lq = opposite_line(lq); + // lq = lq.opposite(); + } + + if ( p.point() == t.source() || p.point() == t.target() ) { + // Line lqperp = lq.perpendicular(p.point()); + Line lqperp = compute_perpendicular(lq, p.point()); + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + // Oriented_side opqr = lqperp.oriented_side(vpqr.point()); + // Oriented_side oqps = lqperp.oriented_side(vqps.point()); + + Oriented_side opqr = vpqr.oriented_side(lqperp, tag); + Oriented_side oqps = vqps.oriented_side(lqperp, tag); + + bool on_different_parabola_arcs = + ((opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE) || + (opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE)); + + if ( !on_different_parabola_arcs ) { return true; } + + Point t1; + if ( p.point() == t.source() ) { + t1 = t.target(); + } else { + t1 = t.source(); + } + + // Oriented_side o_p = lq.oriented_side(p.point()); + // Oriented_side o_t1 = lq.oriented_side(t1); + + Oriented_side o_p = oriented_side_of_line(lq, p.point()); + Oriented_side o_t1 = oriented_side_of_line(lq, t1); + + if ( (o_p == ON_POSITIVE_SIDE && o_t1 == ON_NEGATIVE_SIDE) || + (o_p == ON_NEGATIVE_SIDE && o_t1 == ON_POSITIVE_SIDE) ) { + return true; + } + + // FT d2_p_from_lq = squared_distance_to_line(p.point(), lq); + // FT d2_t1_from_lq = squared_distance_to_line(t1, lq); + + // Comparison_result res = + // CGAL_NTS compare(d2_p_from_lq, d2_t1_from_lq); + + Comparison_result res = + compare_squared_distances_to_line(lq, p.point(), t1); + + return ( res == LARGER ); + } + + // if ( lt.has_on_negative_side(p.point()) ) { + if ( oriented_side_of_line(lt, p.point()) == + ON_NEGATIVE_SIDE) { + // lt = lt.opposite(); + lt = opposite_line(lt); + } + + + Comparison_result res = + CGAL_NTS compare(lt.a() * lq.b(), lt.b() * lq.a()); + bool are_parallel = (res == EQUAL); + + if ( are_parallel ) { + Sign sgn = CGAL_NTS sign(lt.a() * lq.a() + lt.b() * lq.b()); + bool have_opposite_directions = (sgn == NEGATIVE); + if ( have_opposite_directions ) { + // lq = lq.opposite(); + lq = opposite_line(lq); + } + + // if ( lq.oriented_side(p.point()) == + // lt.oriented_side(p.point()) ) { + if ( oriented_side_of_line(lq, p.point()) == + oriented_side_of_line(lt, p.point()) ) { + return true; + } + if ( have_opposite_directions ) { + lq = opposite_line(lq); + // lq = lq.opposite(); + } + } + + // Line l = lt.perpendicular(p.point()); + Line l = compute_perpendicular(lt, p.point()); + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + // Oriented_side o_l_pqr = l.oriented_side(vpqr.point()); + // Oriented_side o_l_qps = l.oriented_side(vqps.point()); + + Oriented_side o_l_pqr = vpqr.oriented_side(l, tag); + Oriented_side o_l_qps = vqps.oriented_side(l, tag); + + + if ( o_l_pqr == ON_POSITIVE_SIDE && + o_l_qps == ON_NEGATIVE_SIDE ) { return false; } + if ( o_l_pqr == ON_NEGATIVE_SIDE && + o_l_qps == ON_POSITIVE_SIDE ) { return true; } + + //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //>>>>>>>>>> HERE I NEED TO CHECK THE BOUNDARY CASES <<<<<< + //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + // Line lqperp = lq.perpendicular(p.point()); + Line lqperp = compute_perpendicular(lq, p.point()); + + // Oriented_side opqr = lqperp.oriented_side(vpqr.point()); + // Oriented_side oqps = lqperp.oriented_side(vqps.point()); + + Oriented_side opqr = vpqr.oriented_side(lqperp, tag); + Oriented_side oqps = vqps.oriented_side(lqperp, tag); + + bool on_different_parabola_arcs = + ((opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE) || + (opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE)); + + if ( !on_different_parabola_arcs ) { return true; } + + + Point pv = lq.projection(p.point()); + pv = midpoint(pv, p.point()); + + Oriented_side o_l_pv = oriented_side_of_line(l, pv); + + CGAL_assertion( o_l_pv != ON_ORIENTED_BOUNDARY ); + + CGAL_assertion( o_l_pqr != ON_ORIENTED_BOUNDARY || + o_l_qps != ON_ORIENTED_BOUNDARY ); + + if ( o_l_pqr == ON_ORIENTED_BOUNDARY ) { + return ( o_l_qps == o_l_pv ); + } else { + return ( o_l_pqr == o_l_pv ); + } + + } + } else { + // q is a point p is a segment + return is_interior_in_conflict_both(q, p, s, r, t, tag); + } + } + +#endif + + //------------------------------------------------------------------------ + //------------------------------------------------------------------------ + //------------------------------------------------------------------------ + + bool + is_interior_in_conflict_touch(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + // checks if interior of voronoi edge is in conflict if both extrema + // of the voronoi edge touch the corresponding circles. + // return true if interior is in conflict; false otherwise + + if ( t.is_segment() ) { return false; } + + if ( (p.is_point() && q.is_point()) || + (p.is_segment() && q.is_segment()) ) { + return true; + } + + if ( p.is_point() && q.is_segment() ) { + Line lq = compute_supporting_line(q.segment()); + +#if 0 + std::cout << "p: " << p << std::endl; + std::cout << "q: " << q << std::endl; + std::cout << "r: " << r << std::endl; + std::cout << "s: " << s << std::endl; + std::cout << "t: " << t << std::endl; + std::cout << std::endl; +#endif + + Comparison_result res = + compare_squared_distances_to_line(lq, p.point(), t.point()); + + return (res != SMALLER); + } + + return is_interior_in_conflict_touch(q, p, s, r, t, tag); + } + + //------------------------------------------------------------------------ + //------------------------------------------------------------------------ + //------------------------------------------------------------------------ + + +#if 1 + + bool + is_interior_in_conflict_none(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Method_tag tag) const + { + if ( t.is_segment() ) { return false; } + + Point tp = t.point(); + + if ( p.is_point() && q.is_point() ) { + return + is_interior_in_conflict_none(p.point(), q.point(), + r, s, tp, tag); + } else if ( p.is_point() && q.is_segment() ) { + return + is_interior_in_conflict_none(p.point(), q.segment(), + r, s, tp, tag); + } else if ( p.is_segment() && q.is_point() ) { + return + is_interior_in_conflict_none(p.segment(), q.point(), + r, s, tp, tag); + } else { // both p and q are segments + return + is_interior_in_conflict_none(p.segment(), q.segment(), + r, s, tp, tag); + } + } +#endif + + //------------------------------------------------------------------------ + + + bool + is_interior_in_conflict_none(const Point& p, const Point& q, + const Site& r, const Site& s, + const Point& t, Method_tag tag) const + { + return false; + } + + //------------------------------------------------------------------------ + + bool + is_interior_in_conflict_none(const Point& p, const Segment& q, + const Site& r, const Site& s, + const Point& t, Method_tag tag) const + { + if ( p == q.source() || p == q.target() ) { + return false; + } + +#if 0 + std::cout << "p: " << p << std::endl; + std::cout << "q: " << q << std::endl; + std::cout << "r: " << r << std::endl; + std::cout << "s: " << s << std::endl; + std::cout << "t: " << t << std::endl; +#endif + + Line lq = compute_supporting_line(q); + + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + +#if 0 + std::cout << "Vpqr: " << vpqr.point() << std::endl; + std::cout << "Vqps: " << vqps.point() << std::endl; +#endif + + Line lperp = compute_perpendicular(lq, t); + + Oriented_side op = oriented_side_of_line(lq, p); + Oriented_side ot = oriented_side_of_line(lq, t); + +#if 0 + std::cout << "op: " << int(op) << std::endl; + std::cout << "ot: " << int(ot) << std::endl; +#endif + + bool on_same_side = + ((op == ON_POSITIVE_SIDE && ot == ON_POSITIVE_SIDE) || + (op == ON_NEGATIVE_SIDE && ot == ON_NEGATIVE_SIDE)); + + Comparison_result res = + compare_squared_distances_to_line(lq, t, p); + +#if 0 + std::cout << "res: " << res << std::endl; +#endif + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + bool on_different_side = + ((opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE) || + (opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE)); + + +#if 0 + std::cout << "opqr: " << int(opqr) << std::endl; + std::cout << "oqps: " << int(oqps) << std::endl; + + std::cout << std::endl; +#endif + + return ( on_same_side && (res == SMALLER) && on_different_side ); + } + + //------------------------------------------------------------------------ + + bool + is_interior_in_conflict_none(const Segment& p, const Point& q, + const Site& r, const Site& s, + const Point& t, Method_tag tag) const + { + return is_interior_in_conflict_none(q, p, s, r, t, tag); + } + + //------------------------------------------------------------------------ + + bool + is_interior_in_conflict_none(const Segment& p, const Segment& q, + const Site& r, const Site& s, + const Point& t, Method_tag tag) const + { + Site sp(p), sq(q); + Voronoi_vertex vpqr(sp, sq, r); + Voronoi_vertex vqps(sq, sp, s); + + Line lp = compute_supporting_line(p); + Line lq = compute_supporting_line(q); + + // first orient lp according to its Voronoi vertices + if ( ( vpqr.is_degenerate_Voronoi_circle() && + vpqr.point() == p.source() ) || + ( vpqr.is_degenerate_Voronoi_circle() && + vpqr.point() == p.target() ) ) { + // CGAL_assertion + // ( !vqps.is_same_point(p.source(), tag) && + // !vqps.is_same_point(p.target(), tag) ); + if ( vqps.oriented_side(lp, tag) != ON_POSITIVE_SIDE ) { + lp = opposite_line(lp); + } + } else { + if ( vpqr.oriented_side(lp, tag) != ON_POSITIVE_SIDE ) { + lp = opposite_line(lp); + } + } + + // then orient lq according to its Voronoi vertices + if ( ( vpqr.is_degenerate_Voronoi_circle() && + vpqr.point() == q.source() ) || + ( vpqr.is_degenerate_Voronoi_circle() && + vpqr.point() == q.target() ) ) { + // CGAL_assertion + // ( !vqps.is_same_point(q.source(), tag) && + // !vqps.is_same_point(q.target(), tag) ); + if ( vqps.oriented_side(lq, tag) != ON_POSITIVE_SIDE ) { + lq = opposite_line(lq); + } + } else { + if ( vpqr.oriented_side(lq, tag) != ON_POSITIVE_SIDE ) { + lq = opposite_line(lq); + } + } + + // check if t is on the same side as the Voronoi vertices + Oriented_side ot_lp = oriented_side_of_line(lp, t); + Oriented_side ot_lq = oriented_side_of_line(lq, t); + + if ( ot_lp != ON_POSITIVE_SIDE || ot_lq != ON_POSITIVE_SIDE ) { + return false; + } + + Line lperp; + + Comparison_result res = + compare_squared_distances_to_lines(t, lp, lq); + +#if 0 + FT d2_t_from_lp = squared_distance_to_line(t.point(), lp); + FT d2_t_from_lq = squared_distance_to_line(t.point(), lq); + + Comparison_result res = + CGAL_NTS compare(d2_t_from_lp, d2_t_from_lq); +#endif + + if ( res == SMALLER ) { + lperp = compute_perpendicular(lp, t); + } else { + lperp = compute_perpendicular(lq, t); + } + + CGAL_precondition( ot_lp != ON_ORIENTED_BOUNDARY && + ot_lq != ON_ORIENTED_BOUNDARY ); + + // check of lperp separates the two Voronoi vertices + Oriented_side opqr_perp = vpqr.oriented_side(lperp, tag); + Oriented_side oqps_perp = vqps.oriented_side(lperp, tag); + + bool on_different_side = + (opqr_perp == ON_POSITIVE_SIDE && + oqps_perp == ON_NEGATIVE_SIDE) || + (opqr_perp == ON_NEGATIVE_SIDE && + oqps_perp == ON_POSITIVE_SIDE); + + return ( on_different_side ); + } + + //------------------------------------------------------------------------ + + + +#if 0 + bool is_interior_in_conflict_none(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Ring_tag tag) const + { + // THIS NEEDS TO BE CORRECTED + return + is_interior_in_conflict_none(p, q, r, s, t, Sqrt_field_tag()); + } + + + + bool + is_interior_in_conflict_none(const Site& p, const Site& q, + const Site& r, const Site& s, + const Site& t, Sqrt_field_tag tag) const + { + // checks if interior of voronoi edge is in conflict if both extrema + // of the voronoi edge are not in conclift + // return true if interior is in conflict; false otherwise + + if ( t.is_segment() ) { return false; } + + if ( p.is_point() && q.is_point() ) { return false; } + + if ( p.is_point() && q.is_segment() ) { + if ( p.point() == q.source() || p.point() == q.target() ) { + return false; + } + + // Line lq(q.segment()); + + Line lq = compute_supporting_line(q.segment()); + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + // Line lperp = lq.perpendicular(t.point()); + Line lperp = compute_perpendicular(lq, t.point()); + + // Oriented_side op = lq.oriented_side(p.point()); + // Oriented_side ot = lq.oriented_side(t.point()); + + Oriented_side op = oriented_side_of_line(lq, p.point()); + Oriented_side ot = oriented_side_of_line(lq, t.point()); + + bool on_same_side = + ((op == ON_POSITIVE_SIDE && ot == ON_POSITIVE_SIDE) || + (op == ON_NEGATIVE_SIDE && ot == ON_NEGATIVE_SIDE)); + +#if 0 + FT d2_t_from_lq = squared_distance_to_line(t.point(), lq); + FT d2_p_from_lq = squared_distance_to_line(p.point(), lq); + + Comparison_result res = + CGAL_NTS compare(d2_t_from_lq, d2_p_from_lq); +#endif + + Comparison_result res = + compare_squared_distances_to_line(lq, t.point(), p.point()); + + // Oriented_side opqr = lperp.oriented_side(vpqr.point()); + // Oriented_side oqps = lperp.oriented_side(vqps.point()); + + Oriented_side opqr = vpqr.oriented_side(lperp, tag); + Oriented_side oqps = vqps.oriented_side(lperp, tag); + + bool on_different_side = + ((opqr == ON_POSITIVE_SIDE && oqps == ON_NEGATIVE_SIDE) || + (opqr == ON_NEGATIVE_SIDE && oqps == ON_POSITIVE_SIDE)); + + + return ( on_same_side && (res == SMALLER) && + on_different_side ); + } else if ( p.is_segment() && q.is_point() ) { + return is_interior_in_conflict_none(q, p, s, r, t, tag); + } else { + // both p and q are segments + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + Line lp = compute_supporting_line(p.segment()); + Line lq = compute_supporting_line(q.segment()); + + // first orient lp according to its Voronoi vertices + if ( vpqr.is_same_point(p.source(), tag) || + vpqr.is_same_point(p.target(), tag) ) { + CGAL_assertion + ( !vqps.is_same_point(p.source(), tag) && + !vqps.is_same_point(p.target(), tag) ); + // if ( lp.oriented_side(vqps.point()) != ON_POSITIVE_SIDE ) { + if ( vqps.oriented_side(lp, tag) != ON_POSITIVE_SIDE ) { + lp = opposite_line(lp); + } + } else { + if ( vpqr.oriented_side(lp, tag) != ON_POSITIVE_SIDE ) { + lp = opposite_line(lp); + } + } + + // then orient lq according to its Voronoi vertices + if ( vpqr.is_same_point(q.source(), tag) || + vpqr.is_same_point(q.target(), tag) ) { + CGAL_assertion + ( !vqps.is_same_point(q.source(), tag) && + !vqps.is_same_point(q.target(), tag) ); + if ( vqps.oriented_side(lq, tag) != ON_POSITIVE_SIDE ) { + lq = opposite_line(lq); + } + } else { + if ( vpqr.oriented_side(lq, tag) != ON_POSITIVE_SIDE ) { + lq = opposite_line(lq); + } + } + + // check if t is on the same side as the Voronoi vertices + Oriented_side ot_lp = oriented_side_of_line(lp, t.point()); + Oriented_side ot_lq = oriented_side_of_line(lq, t.point()); + + if ( ot_lp != ON_POSITIVE_SIDE || + ot_lq != ON_POSITIVE_SIDE ) { + return false; + } + + Line lperp; + +#if 0 + FT d2_t_from_lp = squared_distance_to_line(t.point(), lp); + FT d2_t_from_lq = squared_distance_to_line(t.point(), lq); + + Comparison_result res = + CGAL_NTS compare(d2_t_from_lp, d2_t_from_lq); +#endif + + Comparison_result res = + compare_squared_distances_to_lines(t.point(), lp, lq); + + if ( res == SMALLER ) { + lperp = compute_perpendicular(lp, t.point()); + } else { + lperp = compute_perpendicular(lq, t.point()); + } + + CGAL_precondition( ot_lp != ON_ORIENTED_BOUNDARY && + ot_lq != ON_ORIENTED_BOUNDARY ); + + // check of lperp separates the two Voronoi vertices + // Oriented_side opqr_perp = lperp.oriented_side(vpqr.point()); + // Oriented_side oqps_perp = lperp.oriented_side(vqps.point()); + + Oriented_side opqr_perp = vpqr.oriented_side(lperp, tag); + Oriented_side oqps_perp = vqps.oriented_side(lperp, tag); + + + + bool on_different_side = + (opqr_perp == ON_POSITIVE_SIDE && + oqps_perp == ON_NEGATIVE_SIDE) || + (opqr_perp == ON_NEGATIVE_SIDE && + oqps_perp == ON_POSITIVE_SIDE); + + return ( on_different_side ); + } + } +#endif + + //------------------------------------------------------------------------ + +public: + inline + bool operator()(const Site& p, const Site& q, const Site& r, + const Site& s, const Site& t, Sign sgn) const + { +#if PRED_PRINT + std::cout << "finite edge interior (p,q,r,s,t): " << std::endl; + std::cout << p << std::endl; + std::cout << q << std::endl; + std::cout << r << std::endl; + std::cout << s << std::endl; + std::cout << t << std::endl; + std::cout << sgn << std::endl; + + Voronoi_vertex vpqr(p, q, r); + Voronoi_vertex vqps(q, p, s); + + std::cout << "---" << std::endl; + std::cout << "Voronoi vertex pqr: " << vpqr.point() << std::endl; + std::cout << "Voronoi vertex qps: " << vqps.point() << std::endl; + std::cout << "---" << std::endl; + +#endif + + bool res; + if ( sgn == POSITIVE ) { + res = is_interior_in_conflict_none(p, q, r, s, t, Method_tag()); + } else if ( sgn == NEGATIVE ) { + res = is_interior_in_conflict_both(p, q, r, s, t, Method_tag()); + } else { + res = is_interior_in_conflict_touch(p, q, r, s, t, Method_tag()); + } +#if PRED_PRINT + std::cout << "result: " << res << std::endl; +#endif + + if ( sgn == POSITIVE ) { + return is_interior_in_conflict_none(p, q, r, s, t, Method_tag()); + } else if ( sgn == NEGATIVE ) { + return is_interior_in_conflict_both(p, q, r, s, t, Method_tag()); + } else { + return is_interior_in_conflict_touch(p, q, r, s, t, Method_tag()); + } + } + + + inline + bool operator()(const Site& p, const Site& q, const Site& r, + const Site& t, Sign sgn) const + { +#if PRED_PRINT + std::cout << "finite edge interior (p,q,r,t): " << std::endl; + std::cout << p << std::endl; + std::cout << q << std::endl; + std::cout << r << std::endl; + std::cout << t << std::endl; + std::cout << sgn << std::endl; +#endif + + + + if ( t.is_point() ) { +#if PRED_PRINT + std::cout << false << std::endl; +#endif + return ( sgn == NEGATIVE ); + } + + if ( sgn != NEGATIVE ) { +#if PRED_PRINT + std::cout << false << std::endl; +#endif + return false; + } + + if ( p.is_segment() || q.is_segment() ) { +#if PRED_PRINT + std::cout << false << std::endl; +#endif + return false; + } + + bool p_is_endpoint = (p.point() == t.source() || p.point() == t.target()); + bool q_is_endpoint = (q.point() == t.source() || q.point() == t.target()); + +#if PRED_PRINT + std::cout << ( p_is_endpoint && q_is_endpoint ) << std::endl; +#endif + + return ( p_is_endpoint && q_is_endpoint ); + } + + inline + bool operator()(const Site& p, const Site& q, const Site& t, + Sign sgn) const + { +#if PRED_PRINT + std::cout << "finite edge interior (p,q,t): " << std::endl; + std::cout << p << std::endl; + std::cout << q << std::endl; + std::cout << t << std::endl; + std::cout << sgn << std::endl; +#endif + + + if ( sgn != ZERO ) { +#if PRED_PRINT + std::cout << false << std::endl; +#endif + return false; + } + + if ( p.is_segment() || q.is_segment()) { +#if PRED_PRINT + std::cout << false << std::endl; +#endif + return false; + } + + // both p and q are points + + if ( t.is_point() ) { + RT dtpx = p.point().x() - t.point().x(); + RT dtpy = p.point().y() - t.point().y(); + RT dtqx = q.point().x() - t.point().x(); + RT dtqy = q.point().y() - t.point().y(); + + Sign s1 = sign_of_determinant2x2(dtpx, -dtpy, dtqy, dtqx); + + // Sign s1 = sign_of_determinant2x2(tp.x(), -tp.y(), + // tq.y(), tq.x()); + + CGAL_assertion( s1 != ZERO ); +#if PRED_PRINT + std::cout << ( s1 == NEGATIVE ) << std::endl; +#endif + return ( s1 == NEGATIVE ); + } + + bool bp = ( (p.point() == t.source()) || (p.point() == t.target()) ); + bool bq = ( (q.point() == t.source()) || (q.point() == t.target()) ); + +#if PRED_PRINT + std::cout << ( bp && bq ) << std::endl; +#endif + + return ( bp && bq ); + } + +}; + + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + + +template +class Svd_infinite_edge_interior_2 +{ +public: + typedef typename R::Point_2 Point; + typedef typename R::Segment_2 Segment; + typedef typename R::Line_2 Line; + typedef typename R::Vector_2 Vector; + typedef typename CGAL::Svd_voronoi_vertex_2 Voronoi_vertex; + typedef typename R::Site_2 Site; + typedef typename R::FT FT; + +public: + inline + bool operator()(const Site& q, const Site& s, const Site& r, + const Site& t, Sign sgn) + { + +#if PRED_PRINT + std::cout << "infinite edge interior(q, s, r): " << std::endl; + std::cout << q << std::endl; + std::cout << s << std::endl; + std::cout << r << std::endl; + std::cout << t << std::endl; + std::cout << sgn << std::endl; +#endif + + + if ( t.is_segment() ) { +#if PRED_PRINT + std::cout << false << std::endl; + return false; +#endif + } + +#if 0 + if ( q.is_segment() ) { + // in this case r and s must be endpoints of q + return ( sgn == NEGATIVE ); + } +#endif + +#if PRED_PRINT + std::cout << ( sgn == NEGATIVE ) << std::endl; +#endif + return ( sgn == NEGATIVE ); + } + +}; + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + + +template +class Svd_is_degenerate_edge_test_2 +{ +public: + typedef typename R::Point_2 Point; + typedef typename R::Segment_2 Segment; + typedef typename R::Line_2 Line; + typedef typename R::Vector_2 Vector; + typedef typename CGAL::Svd_voronoi_vertex_2 Voronoi_vertex; + typedef typename R::Site_2 Site; + typedef typename R::FT FT; + +public: + inline + bool operator()(const Site& p, const Site& q, const Site& r, + const Site& s) + { + return false; + } + +}; + +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- + + +CGAL_END_NAMESPACE + + +#endif // CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_2_H diff --git a/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_ftC2.h b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_ftC2.h new file mode 100644 index 00000000000..03fc49ff6f6 --- /dev/null +++ b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Segment_Voronoi_diagram_predicates_ftC2.h @@ -0,0 +1,725 @@ +#ifndef CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H +#define CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H + + +#include +#include + + + +CGAL_BEGIN_NAMESPACE + +template +Comparison_result +svd_compare_distanceC2(const FT& qx, const FT& qy, + const FT& sx, const FT& sy, + const FT& tx, const FT& ty, + const FT& px, const FT& py, Method_tag) +{ + // first check if (qx,qy) is inside, the boundary or at the exterior + // of the band of the segment s + + must_be_filtered(qx); + + FT dx = sx - tx; + FT dy = sy - ty; + + + FT d1x = qx - sx; + FT d1y = qy - sy; + FT d2x = qx - tx; + FT d2y = qy - ty; + + if ( px == sx && py == sy ) { + FT o = dx * d1x + dy * d1y; + if ( o >= 0 ) { + return LARGER; + } + } + + if ( px == tx && py == ty ) { + FT o = dx * d2x + dy * d2y; + if ( o <= 0 ) { + return LARGER; + } + } + + + FT d2_from_p = CGAL_NTS square (qx-px) + CGAL_NTS square(qy-py); + + FT dot1 = dx * d1x + dy * d1y; + if ( dot1 >= 0 ) { + // q is outside (or the boundary of) the band on the side of s. + FT d2_from_s = CGAL_NTS square(d1x) + CGAL_NTS square(d1y); + return CGAL_NTS compare(d2_from_s, d2_from_p); + } + + FT dot2 = dx * d2x + dy * d2y; + if ( dot2 <= 0 ) { + // q is outside (or the boundary of) the band on the side of t. + FT d2_from_t = CGAL_NTS square(d2x) + CGAL_NTS square(d2y); + return CGAL_NTS compare(d2_from_t, d2_from_p); + } + + // q is strictly in the interior of the band, so I have to compare + // its distance from the supporting line. + FT c = sx * ty - sy * tx; + FT n = CGAL_NTS square(dx) + CGAL_NTS square(dy); + FT d2_from_l = CGAL_NTS square(qx * dy - qy * dx + c); + return CGAL_NTS compare(d2_from_l, d2_from_p * n); +} + +template +Comparison_result +svd_compare_distanceC2(const FT& qx, const FT& qy, + const FT& s1x, const FT& s1y, + const FT& t1x, const FT& t1y, + const FT& s2x, const FT& s2y, + const FT& t2x, const FT& t2y, Method_tag) +{ + // first check if (qx,qy) is inside, the boundary or at the exterior + // of the band of the segments s1, s2 + + must_be_filtered(qx); + + FT d1x = s1x - t1x; + FT d1y = s1y - t1y; + + FT d2x = s2x - t2x; + FT d2y = s2y - t2y; + + FT dqs1x = qx - s1x; + FT dqs1y = qy - s1y; + FT dqt1x = qx - t1x; + FT dqt1y = qy - t1y; + + FT dqs2x = qx - s2x; + FT dqs2y = qy - s2y; + FT dqt2x = qx - t2x; + FT dqt2y = qy - t2y; + + FT dot_s1 = d1x * dqs1x + d1y * dqs1y; + int idx1; // 0 for s1, 1 for interior of 1, 2 for t1; + + if ( qx == s1x && qy == s1y ) { + idx1 = 0; + } else if ( qx == t1x && qy == t1y ) { + idx1 = 2; + } else if ( dot_s1 >= 0 ) { + // q is outside (or the boundary of) the band of 1 on the side of s1. + idx1 = 0; + } else { + FT dot_t1 = d1x * dqt1x + d1y * dqt1y; + if ( dot_t1 <= 0 ) { + // q is outside (or the boundary of) the band of 1 on the side of t1. + idx1 = 2; + } else { + idx1 = 1; + } + } + + FT dot_s2 = d2x * dqs2x + d2y * dqs2y; + int idx2; // 0 for s2, 1 for interior of 2, 2 for t2; + + if ( qx == s2x && qy == s2y ) { + idx2 = 0; + } else if ( qx == t2x && qy == t2y ) { + idx2 = 2; + } else if ( dot_s2 >= 0 ) { + // q is outside (or the boundary of) the band of 2 on the side of s2. + idx2 = 0; + } else { + FT dot_t2 = d2x * dqt2x + d2y * dqt2y; + if ( dot_t2 <= 0 ) { + // q is outside (or the boundary of) the band of 2 on the side of t2. + idx2 = 2; + } else { + idx2 = 1; + } + } + + if ( idx1 == 0 ) { + FT d2_from_s1 = CGAL_NTS square(dqs1x) + CGAL_NTS square(dqs1y); + // if ( qx == s1x && qy == s1y ) { d2_from_s1 = FT(0); } + if ( idx2 == 0 ) { + FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y); + // if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); } + + if ( s1x == s2x && s1y == s2y ) { return EQUAL; } + + // std::cout << "checkpoint 2.1" << std::endl; + + + return CGAL_NTS compare(d2_from_s1, d2_from_s2); + } else if ( idx2 == 2 ) { + + // std::cout << "checkpoint 2.2" << std::endl; + + + + FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y); + // if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); } + + if ( s1x == t2x && s1y == t2y ) { return EQUAL; } + + return CGAL_NTS compare(d2_from_s1, d2_from_t2); + } else { // idx2 == 1 + FT c2 = s2x * t2y - s2y * t2x; + FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y); + FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2); + + // std::cout << "checkpoint 2.3" << std::endl; + + return CGAL_NTS compare(d2_from_s1 * n2, d2_from_l2); + } + } else if ( idx1 == 2 ) { + FT d2_from_t1 = CGAL_NTS square(dqt1x) + CGAL_NTS square(dqt1y); + // if ( qx == t1x && qy == t1y ) { d2_from_t1 = FT(0); } + + if ( idx2 == 0 ) { + FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y); + // if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); } + + /* + std::cout << "checkpoint 3.1" << std::endl; + + std::pair ivl; + + ivl = to_interval(d2_from_t1); + std::cout << "interval of d2_from_t1: " + << ivl.first << " " << ivl.second << std::endl; + + ivl = to_interval(d2_from_s2); + std::cout << "interval of d2_from_s2: " + << ivl.first << " " << ivl.second << std::endl; + */ + if ( t1x == s2x && t1y == s2y ) { return EQUAL; } + + return CGAL_NTS compare(d2_from_t1, d2_from_s2); + } else if ( idx2 == 2 ) { + FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y); + // if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); } + + // std::cout << "checkpoint 3.2" << std::endl; + + if ( t1x == t2x && t1y == t2y ) { return EQUAL; } + + return CGAL_NTS compare(d2_from_t1, d2_from_t2); + } else { // idx2 == 1 + FT c2 = s2x * t2y - s2y * t2x; + FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y); + FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2); + + + // std::cout << "checkpoint 3.3" << std::endl; + + + return CGAL_NTS compare(d2_from_t1 * n2, d2_from_l2); + } + } else { // idx1 == 1 + FT c1 = s1x * t1y - s1y * t1x; + FT n1 = CGAL_NTS square(d1x) + CGAL_NTS square(d1y); + FT d2_from_l1 = CGAL_NTS square(qx * d1y - qy * d1x + c1); + if ( idx2 == 0 ) { + FT d2_from_s2 = CGAL_NTS square(dqs2x) + CGAL_NTS square(dqs2y); + // if ( qx == s2x && qy == s2y ) { d2_from_s2 = FT(0); } + + // std::cout << "checkpoint 4.1" << std::endl; + + + return CGAL_NTS compare(d2_from_l1, d2_from_s2 * n1); + } else if ( idx2 == 2 ) { + FT d2_from_t2 = CGAL_NTS square(dqt2x) + CGAL_NTS square(dqt2y); + // if ( qx == t2x && qy == t2y ) { d2_from_t2 = FT(0); } + + // std::cout << "checkpoint 4.2" << std::endl; + + + return CGAL_NTS compare(d2_from_l1, d2_from_t2 * n1); + } else { // idx2 == 1 + FT c2 = s2x * t2y - s2y * t2x; + FT n2 = CGAL_NTS square(d2x) + CGAL_NTS square(d2y); + FT d2_from_l2 = CGAL_NTS square(qx * d2y - qy * d2x + c2); + + // std::cout << "checkpoint 4.3" << std::endl; + + + return CGAL_NTS compare(d2_from_l1 * n2, d2_from_l2 * n1); + } + } +} + + + +//-------------------------------------------------------------------------- + +#if 0 +template +inline +Edge_conflict_test_result +svd_edge_predicate_ftC2(const std::vector& v, + char site_types[], int num_sites, + bool is_degenerate, Method_tag) +{ + FT a; + must_be_filtered(a); + + typedef Simple_cartesian Rep; + typedef Point_2 Point; + typedef Segment_2 Segment; + typedef Segment_Voronoi_diagram_site_2 Site; + + Site t[num_sites]; + + for (int i = 0, k = 0; i < num_sites; i++) { + if ( site_types[i] == 'p' ) { + Point p(v[k], v[k+1]); + t[i].set_point( p ); + } else if ( site_types[i] == 's' ) { + Point p1(v[k], v[k+1]), p2(v[k+2], v[k+3]); + Segment s(p1, p2); + t[i].set_segment( s ); + } else { + CGAL_assertion( site_types[i] == 'p' || + site_types[i] == 's' ); + } + k += ( (site_types[i] == 'p') ? 2 : 4 ); + } + + CGAL_precondition( num_sites >= 3 && num_sites <= 5 ); + if ( num_sites == 3 ) { + return svd_edge_test_degenerated_2(t[0], t[1], t[2]); + } + + if ( num_sites == 4 ) { + if ( is_degenerate ) { + return svd_edge_test_degenerated_2(t[0], t[1], t[2], t[3]); + } else { + return svd_infinite_edge_test_2(t[0], t[1], t[2], t[3]); + } + } + + return svd_edge_test_2(t[0], t[1], t[2], t[3], t[4]); +} +#endif + +//-------------------------------------------------------------------------- + +template +inline +Sign +svd_vertex_conflict_ftC2(const typename K::Site_2 t[], + int num_sites, Method_tag mtag) +{ + typedef typename K::FT FT; + char site_types[num_sites]; + + std::vector v; + + for (int i = 0; i < num_sites; i++) { + if ( t[i].is_point() ) { + v.push_back( t[i].point().x() ); + v.push_back( t[i].point().y() ); + site_types[i] = 'p'; + } else { + v.push_back( t[i].source().x() ); + v.push_back( t[i].source().y() ); + v.push_back( t[i].target().x() ); + v.push_back( t[i].target().y() ); + site_types[i] = 's'; + } + } + + return svd_vertex_conflict_ftC2(v, site_types, num_sites, mtag); +} + + +template +inline +Sign +svd_vertex_conflict_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& t, + Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, t}; + return + svd_vertex_conflict_ftC2(site_vec, 3, mtag); +} + + +template +inline +Sign +svd_vertex_conflict_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& r, + const typename K::Site_2& t, + Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, r, t}; + return + svd_vertex_conflict_ftC2(site_vec, 4, mtag); +} + + +template +inline +Sign +svd_vertex_conflict_ftC2(const std::vector& v, + char site_types[], int num_sites, + Method_tag) +{ + CGAL_precondition( num_sites == 3 || num_sites == 4 ); + + must_be_filtered(FT()); + + typedef Simple_cartesian Rep; + typedef CGAL::Segment_Voronoi_diagram_kernel_wrapper_2 Kernel; + + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Segment_2 Segment_2; + typedef typename Kernel::Site_2 Site_2; + typedef Svd_incircle_2 Incircle; + + + Site_2 t[num_sites]; + + for (int i = 0, k = 0; i < num_sites; i++) { + if ( site_types[i] == 'p' ) { + Point_2 p(v[k], v[k+1]); + t[i].set_point( p ); + } else if ( site_types[i] == 's' ) { + Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]); + Segment_2 s(p1, p2); + t[i].set_segment( s ); + } else { + CGAL_assertion( site_types[i] == 'p' || + site_types[i] == 's' ); + } + k += ( (site_types[i] == 'p') ? 2 : 4 ); + } + + if ( num_sites == 3 ) { + return Incircle()(t[0], t[1], t[2]); + } + + return Incircle()(t[0], t[1], t[2], t[3]); +} + + +//-------------------------------------------------------------------------- + +template +inline +bool +svd_finite_edge_conflict_ftC2(const typename K::Site_2 t[], + Sign sgn, int num_sites, Method_tag mtag) +{ + typedef typename K::FT FT; + char site_types[num_sites]; + + std::vector v; + + for (int i = 0; i < num_sites; i++) { + if ( t[i].is_point() ) { + v.push_back( t[i].point().x() ); + v.push_back( t[i].point().y() ); + site_types[i] = 'p'; + } else { + v.push_back( t[i].source().x() ); + v.push_back( t[i].source().y() ); + v.push_back( t[i].target().x() ); + v.push_back( t[i].target().y() ); + site_types[i] = 's'; + } + } + + return svd_finite_edge_conflict_ftC2(v, sgn, site_types, + num_sites, mtag); +} + + +template +inline +bool +svd_finite_edge_conflict_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& t, + Sign sgn, Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, t}; + return + svd_finite_edge_conflict_ftC2(site_vec, sgn, 3, mtag); +} + + +template +inline +bool +svd_finite_edge_conflict_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& r, + const typename K::Site_2& t, + Sign sgn, Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, r, t}; + return + svd_finite_edge_conflict_ftC2(site_vec, sgn, 4, mtag); +} + + + + + +template +inline +bool +svd_finite_edge_conflict_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& r, + const typename K::Site_2& s, + const typename K::Site_2& t, + Sign sgn, Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, r, s, t}; + return + svd_finite_edge_conflict_ftC2(site_vec, sgn, 5, mtag); +} + +template +inline +bool +svd_finite_edge_conflict_ftC2(const std::vector& v, Sign sgn, + char site_types[], int num_sites, + Method_tag) +{ + CGAL_precondition( num_sites >= 3 || num_sites <= 5 ); + + must_be_filtered(FT()); + + typedef Simple_cartesian Rep; + typedef Segment_Voronoi_diagram_kernel_wrapper_2 Kernel; + + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Segment_2 Segment_2; + typedef typename Kernel::Site_2 Site_2; + typedef Svd_finite_edge_interior_2 Edge_interior_test; + + + Site_2 t[num_sites]; + + for (int i = 0, k = 0; i < num_sites; i++) { + if ( site_types[i] == 'p' ) { + Point_2 p(v[k], v[k+1]); + t[i].set_point( p ); + } else if ( site_types[i] == 's' ) { + Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]); + Segment_2 s(p1, p2); + t[i].set_segment( s ); + } else { + CGAL_assertion( site_types[i] == 'p' || + site_types[i] == 's' ); + } + k += ( (site_types[i] == 'p') ? 2 : 4 ); + } + + if ( num_sites == 3 ) { + return Edge_interior_test()(t[0], t[1], t[2], sgn); + } else if ( num_sites == 4 ) { + return Edge_interior_test()(t[0], t[1], t[2], t[3], sgn); + } + + return Edge_interior_test()(t[0], t[1], t[2], t[3], t[4], sgn); +} + +//-------------------------------------------------------------------------- + +template +inline +bool +svd_infinite_edge_conflict_ftC2(const typename K::Site_2 t[], + Sign sgn, int num_sites, Method_tag mtag) +{ + typedef typename K::FT FT; + char site_types[num_sites]; + + std::vector v; + + for (int i = 0; i < num_sites; i++) { + if ( t[i].is_point() ) { + v.push_back( t[i].point().x() ); + v.push_back( t[i].point().y() ); + site_types[i] = 'p'; + } else { + v.push_back( t[i].source().x() ); + v.push_back( t[i].source().y() ); + v.push_back( t[i].target().x() ); + v.push_back( t[i].target().y() ); + site_types[i] = 's'; + } + } + + return svd_infinite_edge_conflict_ftC2(v, sgn, site_types, + num_sites, mtag); +} + + +template +inline +bool +svd_infinite_edge_conflict_ftC2(const typename K::Site_2& q, + const typename K::Site_2& r, + const typename K::Site_2& s, + const typename K::Site_2& t, + Sign sgn, Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {q, r, s, t}; + return + svd_infinite_edge_conflict_ftC2(site_vec, sgn, 4, mtag); +} + + +template +inline +bool +svd_infinite_edge_conflict_ftC2(const std::vector& v, Sign sgn, + char site_types[], int num_sites, + Method_tag) +{ + CGAL_precondition( num_sites == 4 ); + + must_be_filtered(FT()); + + typedef Simple_cartesian Rep; + typedef Segment_Voronoi_diagram_kernel_wrapper_2 Kernel; + + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Segment_2 Segment_2; + typedef typename Kernel::Site_2 Site_2; + typedef + Svd_infinite_edge_interior_2 Edge_interior_test; + + + Site_2 t[num_sites]; + + for (int i = 0, k = 0; i < num_sites; i++) { + if ( site_types[i] == 'p' ) { + Point_2 p(v[k], v[k+1]); + t[i].set_point( p ); + } else if ( site_types[i] == 's' ) { + Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]); + Segment_2 s(p1, p2); + t[i].set_segment( s ); + } else { + CGAL_assertion( site_types[i] == 'p' || + site_types[i] == 's' ); + } + k += ( (site_types[i] == 'p') ? 2 : 4 ); + } + + + return Edge_interior_test()(t[0], t[1], t[2], t[3], sgn); +} + +//-------------------------------------------------------------------------- + +template +inline +bool +svd_is_degenerate_edge_ftC2(const typename K::Site_2 t[], + int num_sites, Method_tag mtag) +{ + typedef typename K::FT FT; + char site_types[num_sites]; + + std::vector v; + + for (int i = 0; i < num_sites; i++) { + if ( t[i].is_point() ) { + v.push_back( t[i].point().x() ); + v.push_back( t[i].point().y() ); + site_types[i] = 'p'; + } else { + v.push_back( t[i].source().x() ); + v.push_back( t[i].source().y() ); + v.push_back( t[i].target().x() ); + v.push_back( t[i].target().y() ); + site_types[i] = 's'; + } + } + + return svd_is_degenerate_edge_ftC2(v, site_types, + num_sites, mtag); +} + + +template +inline +bool +svd_is_degenerate_edge_ftC2(const typename K::Site_2& p, + const typename K::Site_2& q, + const typename K::Site_2& r, + const typename K::Site_2& t, + Method_tag mtag) +{ + typename K::Site_2 site_vec[] = {p, q, r, t}; + return svd_is_degenerate_edge_ftC2(site_vec, 4, mtag); +} + + +template +inline +bool +svd_is_degenerate_edge_ftC2(const std::vector& v, + char site_types[], int num_sites, + Method_tag) +{ + CGAL_precondition( num_sites == 4 ); + + must_be_filtered(FT()); + + typedef Simple_cartesian Rep; + typedef Segment_Voronoi_diagram_kernel_wrapper_2 Kernel; + + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Segment_2 Segment_2; + typedef typename Kernel::Site_2 Site_2; + typedef + Svd_is_degenerate_edge_test_2 Is_degenerate_edge; + + + Site_2 t[num_sites]; + + for (int i = 0, k = 0; i < num_sites; i++) { + if ( site_types[i] == 'p' ) { + Point_2 p(v[k], v[k+1]); + t[i].set_point( p ); + } else if ( site_types[i] == 's' ) { + Point_2 p1(v[k], v[k+1]), p2(v[k+2], v[k+3]); + Segment_2 s(p1, p2); + t[i].set_segment( s ); + } else { + CGAL_assertion( site_types[i] == 'p' || + site_types[i] == 's' ); + } + k += ( (site_types[i] == 'p') ? 2 : 4 ); + } + + + return Is_degenerate_edge()(t[0], t[1], t[2], t[3]); +} + +//-------------------------------------------------------------------------- + +CGAL_END_NAMESPACE + + +#ifdef CGAL_ARITHMETIC_FILTER_H +#ifndef CGAL_ARITHMETIC_FILTER_SVD_PREDICATES_FTC2_H +#include +#endif // CGAL_ARITHMETIC_FILTER_SVD_PREDICATES_FTC2_H +#endif + +#endif // CGAL_SEGMENT_VORONOI_DIAGRAM_PREDICATES_FTC2_H + diff --git a/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_1.h b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_1.h new file mode 100644 index 00000000000..172db813f8a --- /dev/null +++ b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_1.h @@ -0,0 +1,327 @@ +#ifndef CGAL_SQUARE_ROOT_1_H +#define CGAL_SQUARE_ROOT_1_H + +#include +#include + +#include + +#define CHECK_CGAL_PRECONDITIONS 0 + +CGAL_BEGIN_NAMESPACE + + +template +class Square_root_1 +{ +private: + NT x, y, r; + +private: + typedef Square_root_1 Self; + +public: + typedef NT FT; + typedef NT RT; + + Square_root_1() : x(0), y(0), r(0) {} + Square_root_1(int i) : x(i), y(0), r(0) {} + Square_root_1(const NT& a) : x(a), y(0), r(0) {} + Square_root_1(const NT& a, const NT& b, const NT& c) + : x(a), y(b), r(c) + { + // CGAL_assertion( !(CGAL_NTS is_negative(r)) ); + } + + Square_root_1(const Square_root_1& other) + : x(other.x), y(other.y), r(other.r) {} + + + NT a() const { return x; } + NT b() const { return y; } + NT c() const { return r; } + + Self square() const + { + NT xy = x * y; + + return Self(CGAL_NTS square(x) + CGAL_NTS square(y) * r, + xy + xy, + r); + } + + Sign sign() const + { + Sign sx = CGAL_NTS sign(x); + + if ( CGAL_NTS sign(r) == ZERO ) { return sx; } + + Sign sy = CGAL_NTS sign(y); + + if ( sx == sy ) { return sx; } + if ( sx == ZERO ) { return sy; } + + return Sign( sx * CGAL_NTS compare( CGAL_NTS square(x), + r * CGAL_NTS square(y) ) + ); + } + + Self operator+() const + { + return (*this); + } + + Self operator-() const + { + return Self(-x, -y, r); + } + + double to_double() const + { + // THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS + double xd = CGAL_NTS to_double(x); + double yd = CGAL_NTS to_double(y); + double rd = CGAL_NTS to_double(r); + + return (xd + yd * CGAL_NTS sqrt(rd)); + } + + std::pair to_interval() const + { + // THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS + std::pair x_ivl = CGAL::to_interval(x); + std::pair y_ivl = CGAL::to_interval(y); + std::pair r_ivl = CGAL::to_interval(r); + + std::pair sqrt_r_ivl(CGAL_NTS sqrt(r_ivl.first), + CGAL_NTS sqrt(r_ivl.second)); + + std::pair + ivl(x_ivl.first + y_ivl.first * sqrt_r_ivl.first, + x_ivl.second + y_ivl.second * sqrt_r_ivl.second); + + return ivl; + } +}; + + +// operator * +template +inline +Square_root_1 +operator*(const Square_root_1& x, const NT& n) +{ + return Square_root_1(x.a() * n, x.b() * n, x.c()); +} + + +template +inline +Square_root_1 +operator*(const NT& n, const Square_root_1& x) +{ + return (x * n); +} + +template +inline +Square_root_1 +operator*(const Square_root_1& x, const Square_root_1& y) +{ +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL ); +#endif + + NT a = x.a() * y.a() + x.b() * y.b() * x.c(); + NT b = x.a() * y.b() + x.b() * y.a(); + + return Square_root_1(a, b, x.c()); +} + + +// operator + +template +inline +Square_root_1 +operator+(const Square_root_1& x, const NT& n) +{ + return Square_root_1(x.a() + n, x.b(), x.c()); +} + + +template +inline +Square_root_1 +operator+(const NT& n, const Square_root_1& x) +{ + return (x + n); +} + +template +inline +Square_root_1 +operator+(const Square_root_1& x, const Square_root_1& y) +{ +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL ); +#endif + + return Square_root_1(x.a() + y.a(), x.b() + y.b(), x.c()); +} + + + +// operator - +template +inline +Square_root_1 +operator-(const Square_root_1& x, const NT& n) +{ + return x + (-n); +} + + +template +inline +Square_root_1 +operator-(const NT& n, const Square_root_1& x) +{ + return -(x - n); +} + +template +inline +Square_root_1 +operator-(const Square_root_1& x, const Square_root_1& y) +{ + return (x + (-y)); +} + + + + + +template +inline +std::pair +to_interval(const Square_root_1& x) +{ + return x.to_interval(); +} + + +namespace NTS { + + + template + inline + bool + is_positive(const Square_root_1& x) + { + return sign(x) == POSITIVE; + } + + template + inline + bool + is_negative(const Square_root_1& x) + { + return sign(x) == NEGATIVE; + } + + template + inline + bool + is_zero(const Square_root_1& x) + { + return sign(x) == ZERO; + } + + + template + inline + Sign + sign(const Square_root_1& x) + { + return x.sign(); + } + + template + inline + Square_root_1 + square(const Square_root_1& x) + { + return x.square(); + } + + template + inline + Comparison_result + compare(const Square_root_1& x, const Square_root_1& y) + { +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(x.c(), y.c()) == EQUAL ); +#endif + + Sign s = CGAL_NTS sign(x - y); + + if ( s == ZERO ) { return EQUAL; } + return (s == POSITIVE) ? LARGER : SMALLER; + } + + template + inline + double + to_double(const Square_root_1& x) + { + return x.to_double(); + } + +} // namespace NTS + + + +// operator << +template +inline +Stream& +operator<<(Stream& os, const Square_root_1& x) +{ + os << "(" << x.a() << ")+(" << x.b() << ") sqrt{" << x.c() << "}"; + return os; + +#if 0 + if ( CGAL_NTS is_zero(x) ) { + os << "0"; + return os; + } + + // Square_root_1 One(NT(1), NT(0), x.r()); + + if ( CGAL_NTS sign(x.a()) != ZERO ) { + os << x.a(); + if ( CGAL_NTS is_positive(x.b()) ) { + os << "+"; + } + } + + if ( CGAL_NTS sign(x.b()) != ZERO && + CGAL_NTS sign(x.c()) != ZERO ) { + // if ( CGAL_NTS sign(x.b() - One) != ZERO ) { + os << x.b() << " "; + // } + + os << "sqrt{" << x.c() << "}"; + } + + return os; +#endif +} + + + +CGAL_END_NAMESPACE + + + +#endif // CGAL_SQUARE_ROOT_1_H diff --git a/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_2.h b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_2.h new file mode 100644 index 00000000000..48b20692258 --- /dev/null +++ b/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Square_root_2.h @@ -0,0 +1,379 @@ +#ifndef CGAL_SQUARE_ROOT_2_H +#define CGAL_SQUARE_ROOT_2_H + +#include + + + + +CGAL_BEGIN_NAMESPACE + + +template +class Square_root_2 +{ +private: + typedef Square_root_2 Self; + typedef Square_root_1 Sqrt_1; + + NT a0_, a1_, a2_, a3_; + NT A_, B_; + +public: + typedef NT FT; + typedef NT RT; + +public: + Square_root_2() + : a0_(0), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {} + Square_root_2(int i) + : a0_(i), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {} + Square_root_2(const NT& a) + : a0_(a), a1_(0), a2_(0), a3_(0), A_(0), B_(0) {} + Square_root_2(const NT& a0, const NT& a1, const NT& a2, + const NT& a3, const NT& A, const NT& B) + : a0_(a0), a1_(a1), a2_(a2), a3_(a3), A_(A), B_(B) + { + // CGAL_assertion( !(CGAL_NTS is_negative(A_)) ); + // CGAL_assertion( !(CGAL_NTS is_negative(B_)) ); + } + + Square_root_2(const Square_root_2& other) + : a0_(other.a0_), a1_(other.a1_), a2_(other.a2_), + a3_(other.a3_), A_(other.A_), B_(other.B_) {} + + + NT a() const { return a0_; } + NT b() const { return a1_; } + NT c() const { return a2_; } + NT d() const { return a3_; } + NT e() const { return A_; } + NT f() const { return B_; } + + + Self operator*(const Self& b) const + { +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(A_, b.A_) == EQUAL ); + CGAL_precondition( CGAL_NTS compare(B_, b.B_) == EQUAL ); +#endif + + NT a0 = a0_ * b.a0_ + a1_ * b.a1_ * A_ + a2_ * b.a2_ * B_ + + a3_ * b.a3_ * A_ * B_; + NT a1 = a0_ * b.a1_ + a1_ * b.a0_ + (a2_ * b.a3_ + a3_ * b.a2_) * B_; + NT a2 = a0_ * b.a2_ + a2_ * b.a0_ + (a1_ * b.a3_ + a3_ * b.a1_) * A_; + NT a3 = a0_ * b.a3_ + a3_ * b.a0_ + a1_ * b.a2_ + a2_ * b.a1_; + + return Self(a0, a1, a2, a3, A_, B_); + } + + Self operator-() const + { + return Self(-a0_, -a1_, -a2_, -a3_, A_, B_); + } + + Self operator+() const + { + return (*this); + } + + Self square() const + { + NT a0 = CGAL_NTS square(a0_) + CGAL_NTS square(a1_) * A_ + + CGAL_NTS square(a2_) * B_ + CGAL_NTS square(a3_) * A_ * B_; + NT a1_half = a0_ * a1_ + a2_ * a3_ * B_; + NT a2_half = a0_ * a2_ + a1_ * a3_ * A_; + NT a3_half = a0_ * a3_ + a1_ * a2_; + + NT a1 = a1_half + a1_half; + NT a2 = a2_half + a2_half; + NT a3 = a3_half + a3_half; + + return Self(a0, a1, a2, a3, A_, B_); + } + + Sign sign() const + { + // std::cout << "cp sgn 1" << std::flush; + + Sqrt_1 x(a0_, a1_, A_); + Sqrt_1 y(a2_, a3_, A_); + +#if 0 + std::cout << "a0: " << a0_ << std::endl; + std::cout << "a1: " << a1_ << std::endl; + std::cout << "A : " << A_ << std::endl; + std::cout << "x: " << x << std::endl; + + std::cout << " 2" << std::flush; +#endif + + Sign s_x = CGAL_NTS sign(x); + + // std::cout << " 3" << std::flush; + + Sign s_y = CGAL_NTS sign(y); + + // std::cout << " 4" << std::flush; + + Sign s_B = CGAL_NTS sign(B_); + + // std::cout << " 5" << std::flush; + // std::cout << std::endl; + + if ( s_B == ZERO ) { + return s_x; + } else if ( s_x == s_y ) { + return s_x; + } else if ( s_x == ZERO ) { + return s_y; + } else if ( s_y == ZERO ) { + return s_x; + } else { + Sqrt_1 Q = CGAL_NTS square(x) - CGAL_NTS square(y) * B_; + return Sign(s_x * CGAL_NTS sign(Q)); + } + } + + double to_double() const + { + // THIS MUST BE CHECK WITH SYLVAIN FOR CORRECTNESS + double a0d = CGAL_NTS to_double(a0_); + double a1d = CGAL_NTS to_double(a1_); + double a2d = CGAL_NTS to_double(a2_); + double a3d = CGAL_NTS to_double(a3_); + double Ad = CGAL_NTS to_double(A_); + double Bd = CGAL_NTS to_double(B_); + + return (a0d + a1d * CGAL_NTS sqrt(Ad) + a2d * CGAL_NTS sqrt(Bd) + + a3d * CGAL_NTS sqrt(Ad * Bd)); + } +}; + + + + + +// operator * +template +inline +Square_root_2 +operator*(const Square_root_2& x, const NT& n) +{ + return Square_root_2(x.a() * n, x.b() * n, x.c() * n, + x.d() * n, x.e(), x.f()); +} + + +template +inline +Square_root_2 +operator*(const NT& n, const Square_root_2& x) +{ + return (x * n); +} + + + +// operator + +template +inline +Square_root_2 +operator+(const Square_root_2& x, const NT& n) +{ + return Square_root_2(x.a() + n, x.b(), x.c(), x.d(), + x.e(), x.f()); +} + + +template +inline +Square_root_2 +operator+(const NT& n, const Square_root_2& x) +{ + return (x + n); +} + +template +inline +Square_root_2 +operator+(const Square_root_2& x, const Square_root_2& y) +{ +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(x.e(), y.e()) == EQUAL ); + CGAL_precondition( CGAL_NTS compare(x.f(), y.f()) == EQUAL ); +#endif + + return Square_root_2(x.a() + y.a(), x.b() + y.b(), + x.c() + y.c(), x.d() + y.d(), + x.e(), x.f()); +} + + + +// operator - +template +inline +Square_root_2 +operator-(const Square_root_2& x, const NT& n) +{ + return x + (-n); +} + + +template +inline +Square_root_2 +operator-(const NT& n, const Square_root_2& x) +{ + return -(x - n); +} + +template +inline +Square_root_2 +operator-(const Square_root_2& x, const Square_root_2& y) +{ + return (x + (-y)); +} + + + + + +namespace NTS { + + + template + inline + bool + is_positive(const Square_root_2& x) + { + return sign(x) == POSITIVE; + } + + template + inline + bool + is_negative(const Square_root_2& x) + { + return sign(x) == NEGATIVE; + } + + template + inline + bool + is_zero(const Square_root_2& x) + { + return sign(x) == ZERO; + } + + + template + inline + Sign + sign(const Square_root_2& x) + { + return x.sign(); + } + + template + inline + Square_root_2 + square(const Square_root_2& x) + { + return x.square(); + } + + template + inline + Comparison_result + compare(const Square_root_2& x, + const Square_root_2& y) + { +#if CHECK_CGAL_PRECONDITIONS + CGAL_precondition( CGAL_NTS compare(x.e(), y.e()) == EQUAL ); + CGAL_precondition( CGAL_NTS compare(x.f(), y.f()) == EQUAL ); +#endif + + Sign s = CGAL_NTS sign(x - y); + + if ( s == ZERO ) { return EQUAL; } + return (s == POSITIVE) ? LARGER : SMALLER; + } + +} // namespace NTS + + + +// operator << +template +inline +Stream& +operator<<(Stream& os, const Square_root_2& x) +{ + os << "(" << x.a() << ")+(" << x.b() << ") sqrt{" << x.e() << "}"; + os << "+(" << x.c() << ") sqrt{" << x.f() << "}"; + os << "+(" << x.d() << ") sqrt{(" << x.e() << ") (" << x.f() + << ")}"; + return os; + +#if 0 + if ( CGAL_NTS is_zero(x) ) { + os << "0"; + return os; + } + + Square_root_2 One(NT(1), NT(0), NT(0), NT(0), x.e(), x.f()); + + if ( CGAL_NTS sign(x.a()) != ZERO ) { + os << x.a(); + if ( CGAL_NTS is_positive(x.b()) ) { + os << "+"; + } + } + + if ( CGAL_NTS sign(x.b()) != ZERO && + CGAL_NTS sign(x.e()) != ZERO ) { + if ( CGAL_NTS sign(x.b() - One) != ZERO ) { + os << x.b() << " "; + } + os << "sqrt{" << x.e() << "}"; + if ( CGAL_NTS is_positive(x.c()) ) { + os << "+"; + } + } + + if ( CGAL_NTS sign(x.c()) != ZERO && + CGAL_NTS sign(x.f()) != ZERO ) { + if ( CGAL_NTS sign(x.c() - One) != ZERO ) { + os << x.c() << " "; + } + os << "sqrt{" << x.f() << "}"; + if ( CGAL_NTS is_positive(x.d()) ) { + os << "+"; + } + } + + if ( CGAL_NTS sign(x.d()) != ZERO && + CGAL_NTS sign(x.e()) != ZERO && + CGAL_NTS sign(x.f()) != ZERO ) { + if ( CGAL_NTS sign(x.d() - One) != ZERO ) { + os << x.d() << " "; + } + os << "sqrt{" << x.e() << " " << x.f() << "}"; + } + + return os; +#endif +} + + + + + + +CGAL_END_NAMESPACE + + +#endif // CGAL_SQUARE_ROOT_2_H