From 27a81e581d32b882159bb8fea63c25c00d5fc0c4 Mon Sep 17 00:00:00 2001 From: Mikhail Bogdanov Date: Fri, 8 Mar 2013 19:40:30 +0100 Subject: [PATCH] modified and moved the predicate Is_hyperbolic to the traits --- .../Delaunay_hyperbolic_triangulation_2.h | 136 +-- .../CGAL/Triangulation_hyperbolic_traits_2.h | 810 ++++++++++-------- 2 files changed, 450 insertions(+), 496 deletions(-) diff --git a/Hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h b/Hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h index f79ad2c9156..52e6d4588ec 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h @@ -424,111 +424,33 @@ private: Mark_face test(*this); mark_face(f, test); } - - class Is_hyperbolic - { - public: - Is_hyperbolic() - { - } - bool operator() (const Face_handle& f) const - { - typedef typename Gt::Vector_3 Vector_3; - - Point p0 = f->vertex(0)->point(); - Point p1 = f->vertex(1)->point(); - Point p2 = f->vertex(2)->point(); - - Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), - p1.x()*p1.x() + p1.y()*p1.y(), - p2.x()*p2.x() + p2.y()*p2.y()); - - Vector_3 v1 = Vector_3(p0.x(), p1.x(), p2.x()); - Vector_3 v2 = Vector_3(p0.y(), p1.y(), p2.y()); - Vector_3 v3 = Vector_3(FT(1), FT(1), FT(1)); - - FT dt0 = CGAL::determinant(v0, v1, v3); - FT dt1 = CGAL::determinant(v0, v2, v3); - FT dt2 = CGAL::determinant(v0 - v3, v1, v2); - - return dt0*dt0 + dt1*dt1 - dt2*dt2 < 0; - } - - // assume the incident faces are "non-hyperbolic" - bool operator() (const Edge& e) const - { - typedef typename Gt::Vector_2 Vector_2; - typedef typename Gt::Vector_3 Vector_3; - - // endpoints of the edge - Point p0 = e.first->vertex(cw(e.second))->point(); - Point p1 = e.first->vertex(ccw(e.second))->point(); - - // vertices opposite to p0 and p1 - Point q0 = e.first->vertex(e.second)->point(); - - Face_handle f = e.first->neighbor(e.second); - int ind = f->index(e.first); - Point q1 = f->vertex(ind); - - Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), - p1.x()*p1.x() + p1.y()*p1.y(), - q0.x()*q0.x() + q0.y()*q0.y()); - - Vector_3 v1 = Vector_3(2*p0.y(), 2*p1.y(), 2*q0.y()); - Vector_3 v2 = Vector_3(2*p0.x(), 2*p1.x(), 2*q0.x()); - Vector_3 v3 = Vector_3(FT(-1), FT(-1), FT(-1)); - - FT dt0 = CGAL::determinant(v0, v1, v3); - FT dt1 = CGAL::determinant(v2, v0, v3); - - FT m11 = p0.x() * (p1.x()*p1.x() + p1.y()*p1.y() - 2) - p1.x() * (p0.x()*p0.x() + p0.y()*p0.y() - 2); - FT m12 = p0.y() * (p1.x()*p1.x() + p1.y()*p1.y() - 2) - p1.y() * (p0.x()*p0.x() + p0.y()*p0.y() - 2); - - Vector_3 v4 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), - p1.x()*p1.x() + p1.y()*p1.y(), - q1.x()*q1.x() + q1.y()*q1.y()); - - Vector_3 v5 = Vector_3(2*p0.y(), 2*p1.y(), 2*q1.y()); - Vector_3 v6 = Vector_3(2*p0.x(), 2*p1.x(), 2*q1.x()); - Vector_3 v7 = Vector_3(FT(-1), FT(-1), FT(-1)); - - FT dt3 = CGAL::determinant(v4, v5, v7); - FT dt4 = CGAL::determinant(v6, v4, v7); - - FT big_dt0 = CGAL::determinant(Vector_2(dt0, m11), Vector_2(dt1, m12)); - FT big_dt1 = CGAL::determinant(Vector_2(dt3, m11), Vector_2(dt4, m12)); - - assert(big_dt0 == 0 || big_dt1 == 0); - return (big_dt0 > 0 && big_dt1 > 0) || (big_dt0 < 0 && big_dt1 < 0); - } - - private: - Is_hyperbolic(const Is_hyperbolic&); - Is_hyperbolic& operator= (const Is_hyperbolic&); - }; - class Mark_face { public: - typedef typename Gt::Circle_2 Circle_2; - Mark_face(const Self& tr) : _tr(tr) {} Face_info operator ()(const Face_handle& f) const { + typedef typename Gt::Is_hyperbolic Is_hyperbolic; + Face_info info; if(_tr.has_infinite_vertex(f)) { return info; } - if(Is_hyperbolic()(f) == false) { + Point p0 = f->vertex(0)->point(); + Point p1 = f->vertex(1)->point(); + Point p2 = f->vertex(2)->point(); + int ind = 0; + + Is_hyperbolic is_hyperbolic = _tr.geom_traits().Is_hyperbolic_object(); + if(is_hyperbolic(p0, p1, p2, ind) == false) { info.set_finite_invisible(true); - info.set_invisible_edge(find_invisible_edge(f)); + info.set_invisible_edge(ind); return info; } @@ -537,46 +459,10 @@ private: } private: - - // assume that the circumscribing circle of a given face intersects - // the circle at infinity - unsigned char find_invisible_edge(const Face_handle& f) const - { - typedef Euclidean_geom_traits Egt; - typedef typename Egt::Construct_circumcenter_2 Construct_circumcenter_2; - typedef typename Egt::Direction_2 Direction_2; - - assert(!_tr.has_infinite_vertex(f)); - - Point p0 = f->vertex(0)->point(); - Point p1 = f->vertex(1)->point(); - Point p2 = f->vertex(2)->point(); - - _circle = Circle_2(p0, p1, p2); - Point c = _circle.center(); - - Direction_2 d0(p0.x()-c.x(), p0.y()-c.y()); - Direction_2 d1(p1.x()-c.x(), p1.y()-c.y()); - Direction_2 d2(p2.x()-c.x(), p2.y()-c.y()); - - Direction_2 d(c.x(), c.y()); - - if(d.counterclockwise_in_between(d0, d1)) { - return 2; - } - - if(d.counterclockwise_in_between(d1, d2)) { - return 0; - } - - return 1; - } - + Mark_face(const Mark_face&); Mark_face& operator= (const Mark_face&); - mutable Circle_2 _circle; - const Self& _tr; }; diff --git a/Hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h b/Hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h index 79051b75124..34d324ea838 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h @@ -35,407 +35,475 @@ namespace CGAL { - - - template < class R > - class Triangulation_hyperbolic_traits_2 { +template < class R > +class Triangulation_hyperbolic_traits_2 { +public: + typedef Triangulation_hyperbolic_traits_2 Self; + + typedef R Triangulation_euclidean_traits_2; + + typedef R Rep; + typedef typename R::RT RT; + typedef typename R::Point_2 Point_2; + typedef typename R::Vector_2 Vector_2; + typedef typename R::Triangle_2 Triangle_2; + typedef typename R::Line_2 Line_2; + typedef typename R::Ray_2 Ray_2; + + typedef typename R::Vector_3 Vector_3; + typedef typename R::Point_3 Point_3; + + typedef typename R::Less_x_2 Less_x_2; + typedef typename R::Less_y_2 Less_y_2; + typedef typename R::Compare_x_2 Compare_x_2; + typedef typename R::Compare_y_2 Compare_y_2; + typedef typename R::Orientation_2 Orientation_2; + typedef typename R::Side_of_oriented_circle_2 Side_of_oriented_circle_2; + typedef typename R::Construct_bisector_2 Construct_bisector_2; + typedef typename R::Compare_distance_2 Compare_distance_2; + typedef typename R::Construct_triangle_2 Construct_triangle_2; + typedef typename R::Construct_direction_2 Construct_direction_2; + + typedef typename R::Angle_2 Angle_2; + typedef typename R::Construct_midpoint_2 Construct_midpoint_2; + typedef typename R::Compute_squared_distance_2 Compute_squared_distance_2; + + typedef typename R::Iso_rectangle_2 Iso_rectangle_2; + typedef typename R::Circle_2 Circle_2; + + typedef boost::tuple Arc_2; + typedef typename R::Segment_2 Line_segment_2; + typedef boost::variant Segment_2; + + typedef typename R::Line_2 Euclidean_line_2; + +private: + // Poincaré disk + const Circle_2 _unit_circle; + +public: + const Circle_2& unit_circle() const + { + return _unit_circle; + } + + Angle_2 + angle_2_object() const + { return Angle_2(); } + + Compute_squared_distance_2 + compute_squared_distance_2_object() const + { return Compute_squared_distance_2(); } + + class Construct_segment_2 + { + typedef typename CGAL::Regular_triangulation_filtered_traits_2 Regular_geometric_traits_2; + typedef typename Regular_geometric_traits_2::Construct_weighted_circumcenter_2 Construct_weighted_circumcenter_2; + typedef typename Regular_geometric_traits_2::Weighted_point_2 Weighted_point_2; + typedef typename Regular_geometric_traits_2::Bare_point Bare_point; + public: - typedef Triangulation_hyperbolic_traits_2 Self; + Construct_segment_2(const Circle_2& c) : _unit_circle(c) + { + } - typedef R Triangulation_euclidean_traits_2; - - typedef R Rep; - typedef typename R::RT RT; - typedef typename R::Point_2 Point_2; - typedef typename R::Vector_2 Vector_2; - typedef typename R::Triangle_2 Triangle_2; - typedef typename R::Line_2 Line_2; - typedef typename R::Ray_2 Ray_2; - - typedef typename R::Vector_3 Vector_3; - typedef typename R::Point_3 Point_3; - - typedef typename R::Less_x_2 Less_x_2; - typedef typename R::Less_y_2 Less_y_2; - typedef typename R::Compare_x_2 Compare_x_2; - typedef typename R::Compare_y_2 Compare_y_2; - typedef typename R::Orientation_2 Orientation_2; - typedef typename R::Side_of_oriented_circle_2 Side_of_oriented_circle_2; - typedef typename R::Construct_bisector_2 Construct_bisector_2; - typedef typename R::Compare_distance_2 Compare_distance_2; - typedef typename R::Construct_triangle_2 Construct_triangle_2; - typedef typename R::Construct_direction_2 Construct_direction_2; - - typedef typename R::Angle_2 Angle_2; - typedef typename R::Construct_midpoint_2 Construct_midpoint_2; - typedef typename R::Compute_squared_distance_2 Compute_squared_distance_2; - - typedef typename R::Iso_rectangle_2 Iso_rectangle_2; - typedef typename R::Circle_2 Circle_2; - - typedef boost::tuple Arc_2; - typedef typename R::Segment_2 Line_segment_2; - typedef boost::variant Segment_2; - - typedef typename R::Line_2 Euclidean_line_2; + Segment_2 operator()(const Point_2& p, const Point_2& q) const + { + typedef typename R::Collinear_2 Collinear_2; + if(Collinear_2()(p, q, _unit_circle.center())){ + return Line_segment_2(p, q); + } + + Weighted_point_2 wp(p); + Weighted_point_2 wq(q); + Weighted_point_2 wo(_unit_circle.center(), _unit_circle.squared_radius()); + + Bare_point center = Construct_weighted_circumcenter_2()(wp, wo, wq); + FT radius = Compute_squared_distance_2()(p, center); + + Circle_2 circle( center, radius); + // uncomment!!! + assert(circle.has_on_boundary(p) && circle.has_on_boundary(q)); + + if(Orientation_2()(p, q, center) == LEFT_TURN) { + return Arc_2(circle, p, q); + } + return Arc_2(circle, q, p); + } private: - // Poincaré disk - const Circle_2 _unit_circle; + const Circle_2& _unit_circle; + }; + Construct_segment_2 + construct_segment_2_object() const + { + return Construct_segment_2(_unit_circle); + } + + class Construct_circumcenter_2 + { public: - const Circle_2& unit_circle() const - { - return _unit_circle; - } - - Angle_2 - angle_2_object() const - { return Angle_2(); } - - Compute_squared_distance_2 - compute_squared_distance_2_object() const - { return Compute_squared_distance_2(); } - - class Construct_segment_2 - { - typedef typename CGAL::Regular_triangulation_filtered_traits_2 Regular_geometric_traits_2; - typedef typename Regular_geometric_traits_2::Construct_weighted_circumcenter_2 Construct_weighted_circumcenter_2; - typedef typename Regular_geometric_traits_2::Weighted_point_2 Weighted_point_2; - typedef typename Regular_geometric_traits_2::Bare_point Bare_point; - - public: - Construct_segment_2(const Circle_2& c) : _unit_circle(c) - { - } - - Segment_2 operator()(const Point_2& p, const Point_2& q) const - { - typedef typename R::Collinear_2 Collinear_2; - if(Collinear_2()(p, q, _unit_circle.center())){ - return Line_segment_2(p, q); - } - - Weighted_point_2 wp(p); - Weighted_point_2 wq(q); - Weighted_point_2 wo(_unit_circle.center(), _unit_circle.squared_radius()); - - Bare_point center = Construct_weighted_circumcenter_2()(wp, wo, wq); - FT radius = Compute_squared_distance_2()(p, center); - - Circle_2 circle( center, radius); - // uncomment!!! - //assert(circle.has_on_boundary(p) && circle.has_on_boundary(q)); - - if(Orientation_2()(p, q, center) == LEFT_TURN) { - return Arc_2(circle, p, q); - } - return Arc_2(circle, q, p); - } - - private: - const Circle_2& _unit_circle; - }; - - Construct_segment_2 - construct_segment_2_object() const - { - return Construct_segment_2(_unit_circle); - } - - class Construct_circumcenter_2 - { - public: - Construct_circumcenter_2(const Circle_2& c) : _unit_circle(c) - {} - - // TODO: improve this function - Point_2 operator()(Point_2 p, Point_2 q, Point_2 r) - { - assert(_unit_circle.bounded_side(p) == ON_BOUNDED_SIDE); - assert(_unit_circle.bounded_side(q) == ON_BOUNDED_SIDE); - assert(_unit_circle.bounded_side(r) == ON_BOUNDED_SIDE); - - Circle_2 circle(p, q, r); - // circle must be inside the unit one - assert(CGAL::do_intersect(_unit_circle, circle) == false); - - if(circle.center() <= _unit_circle.center() && circle.center() >= _unit_circle.center()){ - return _unit_circle.center(); - } - - FT x0 = circle.center().x(), y0 = circle.center().y(); - // a*alphaˆ2 + b*alpha + c = 0; - FT a = x0*x0 + y0*y0; - FT b = a - circle.squared_radius() + _unit_circle.squared_radius(); - FT c = _unit_circle.squared_radius(); - FT D = b*b - 4*a*c; - - FT alpha = (b - CGAL::sqrt(to_double(D)))/(2*a); - - Point_2 center(x0*alpha, y0*alpha); - if(!circle.has_on_bounded_side(center)) - { std::cout << "Center does not belong to the pencil of spheres!!!" << std::endl;} ; - return center; - } - - private: - const Circle_2 _unit_circle; - }; - - Construct_circumcenter_2 - construct_circumcenter_2_object() - { - Construct_circumcenter_2(_unit_circle); - } - - Construct_midpoint_2 - construct_midpoint_2_object() const - { return Construct_midpoint_2(); } - - //for natural_neighbor_coordinates_2 - typedef typename R::FT FT; - typedef typename R::Equal_x_2 Equal_x_2; - typedef typename R::Compute_area_2 Compute_area_2; - Compute_area_2 compute_area_2_object () const - { - return Compute_area_2(); - } - - // for compatibility with previous versions - typedef Point_2 Point; - typedef Segment_2 Segment; - typedef Triangle_2 Triangle; - typedef Ray_2 Ray; - //typedef Line_2 Line; - - Triangulation_hyperbolic_traits_2() : - _unit_circle(Point_2(0, 0), 1*1) + Construct_circumcenter_2(const Circle_2& c) : _unit_circle(c) {} - Triangulation_hyperbolic_traits_2(FT r) : - _unit_circle(Point_2(0, 0), r*r) - {} - - Triangulation_hyperbolic_traits_2(const Triangulation_hyperbolic_traits_2 & other) : - _unit_circle(other._unit_circle) - {} - - Triangulation_hyperbolic_traits_2 &operator= - (const Triangulation_hyperbolic_traits_2 &) - { - return *this; - } - - Less_x_2 - less_x_2_object() const - { return Less_x_2();} - - Less_y_2 - less_y_2_object() const - { return Less_y_2();} - - Compare_x_2 - compare_x_2_object() const - { return Compare_x_2();} - - Compare_y_2 - compare_y_2_object() const - { return Compare_y_2();} - - Orientation_2 - orientation_2_object() const - { return Orientation_2();} - - Side_of_oriented_circle_2 - side_of_oriented_circle_2_object() const - {return Side_of_oriented_circle_2();} - - Construct_circumcenter_2 - construct_circumcenter_2_object() const - { - return Construct_circumcenter_2(_unit_circle); - } - - class Construct_hyperbolic_bisector_2 - { - public: - Construct_hyperbolic_bisector_2(const Circle_2& unit_circle) : - _unit_circle(unit_circle) {} + // TODO: improve this function + Point_2 operator()(Point_2 p, Point_2 q, Point_2 r) + { + assert(_unit_circle.bounded_side(p) == ON_BOUNDED_SIDE); + assert(_unit_circle.bounded_side(q) == ON_BOUNDED_SIDE); + assert(_unit_circle.bounded_side(r) == ON_BOUNDED_SIDE); - Segment_2 operator()(Point_2 p, Point_2 q) const - { - // If two points are almost of the same distance to the origin, then - // the bisector is supported by the circle of huge radius etc. - // This circle is computed inexactly. - // At present time, in this case the bisector is supported by the line. + Circle_2 circle(p, q, r); + // circle must be inside the unit one + assert(CGAL::do_intersect(_unit_circle, circle) == false); + + if(circle.center() <= _unit_circle.center() && circle.center() >= _unit_circle.center()){ + return _unit_circle.center(); + } + + FT x0 = circle.center().x(), y0 = circle.center().y(); + // a*alphaˆ2 + b*alpha + c = 0; + FT a = x0*x0 + y0*y0; + FT b = a - circle.squared_radius() + _unit_circle.squared_radius(); + FT c = _unit_circle.squared_radius(); + FT D = b*b - 4*a*c; + + FT alpha = (b - CGAL::sqrt(to_double(D)))/(2*a); + + Point_2 center(x0*alpha, y0*alpha); + if(!circle.has_on_bounded_side(center)) + { std::cout << "Center does not belong to the pencil of spheres!!!" << std::endl;} ; + return center; + } + + private: + const Circle_2 _unit_circle; + }; + + Construct_circumcenter_2 + construct_circumcenter_2_object() + { + Construct_circumcenter_2(_unit_circle); + } + + Construct_midpoint_2 + construct_midpoint_2_object() const + { return Construct_midpoint_2(); } + + //for natural_neighbor_coordinates_2 + typedef typename R::FT FT; + typedef typename R::Equal_x_2 Equal_x_2; + typedef typename R::Compute_area_2 Compute_area_2; + Compute_area_2 compute_area_2_object () const + { + return Compute_area_2(); + } + + // for compatibility with previous versions + typedef Point_2 Point; + typedef Segment_2 Segment; + typedef Triangle_2 Triangle; + typedef Ray_2 Ray; + //typedef Line_2 Line; + + Triangulation_hyperbolic_traits_2() : + _unit_circle(Point_2(0, 0), 1*1) + {} + + Triangulation_hyperbolic_traits_2(FT r) : + _unit_circle(Point_2(0, 0), r*r) + {} + + Triangulation_hyperbolic_traits_2(const Triangulation_hyperbolic_traits_2 & other) : + _unit_circle(other._unit_circle) + {} + + Triangulation_hyperbolic_traits_2 &operator= + (const Triangulation_hyperbolic_traits_2 &) + { + return *this; + } + + Less_x_2 + less_x_2_object() const + { return Less_x_2();} + + Less_y_2 + less_y_2_object() const + { return Less_y_2();} + + Compare_x_2 + compare_x_2_object() const + { return Compare_x_2();} + + Compare_y_2 + compare_y_2_object() const + { return Compare_y_2();} + + Orientation_2 + orientation_2_object() const + { return Orientation_2();} + + Side_of_oriented_circle_2 + side_of_oriented_circle_2_object() const + {return Side_of_oriented_circle_2();} + + Construct_circumcenter_2 + construct_circumcenter_2_object() const + { + return Construct_circumcenter_2(_unit_circle); + } + + class Construct_hyperbolic_bisector_2 + { + public: + Construct_hyperbolic_bisector_2(const Circle_2& unit_circle) : + _unit_circle(unit_circle) {} + + Segment_2 operator()(Point_2 p, Point_2 q) const + { + // If two points are almost of the same distance to the origin, then + // the bisector is supported by the circle of huge radius etc. + // This circle is computed inexactly. + // At present time, in this case the bisector is supported by the line. + + Compute_squared_distance_2 dist = Compute_squared_distance_2(); + Point origin = _unit_circle.center(); + FT dif = dist(origin, p) - dist(origin, q); + FT eps = 0.0000000001; + + // Bisector is straight in euclidean sense + if(dif > -eps && dif < eps){ - Compute_squared_distance_2 dist = Compute_squared_distance_2(); - Point origin = _unit_circle.center(); - FT dif = dist(origin, p) - dist(origin, q); - FT eps = 0.0000000001; - - // Bisector is straight in euclidean sense - if(dif > -eps && dif < eps){ - // ideally //if(Compare_distance_2()(_unit_circle.center(), p, q) == EQUAL){ - // TODO: calling R::Construct_bisector - Euclidean_line_2 l = Construct_bisector_2()(p, q); - // compute the ending points - std::pair points = find_intersection(l); - // TODO: improve - Vector_2 v(points.first, points.second); - if(v*l.to_vector() > 0){ - return Line_segment_2(points.first, points.second); - } - return Line_segment_2(points.second, points.first); - } - - Circle_2 c = construct_supporting_circle(p, q); + // TODO: calling R::Construct_bisector + Euclidean_line_2 l = Construct_bisector_2()(p, q); // compute the ending points - std::pair points = find_intersection(c); - - if(Orientation_2()(points.first, points.second, c.center()) == LEFT_TURN) { - return Arc_2(c, points.first, points.second); + std::pair points = find_intersection(l); + // TODO: improve + Vector_2 v(points.first, points.second); + if(v*l.to_vector() > 0){ + return Line_segment_2(points.first, points.second); } - return Arc_2(c, points.second, points.first); + return Line_segment_2(points.second, points.first); } - private: - // The cirle belongs to the pencil with limit points p and q - Circle_2 construct_supporting_circle(Point_2 p, Point_2 q) const - { - // p, q are zero-circles - // (x, y, xˆ2 + yˆ2 - rˆ2) = alpha*(xp, yp, xpˆ2 + ypˆ2) + (1-alpha)*(xq, yq, xqˆ2 + yqˆ2) - // xˆ2 + yˆ2 - rˆ2 = Rˆ2, where R - is a radius of the given unit circle - FT op = p.x()*p.x() + p.y()*p.y(); - FT oq = q.x()*q.x() + q.y()*q.y(); - FT alpha = (_unit_circle.squared_radius() - oq) / (op - oq); - - FT x = alpha*p.x() + (1-alpha)*q.x(); - FT y = alpha*p.y() + (1-alpha)*q.y(); - FT radius = x*x + y*y - _unit_circle.squared_radius(); - - //improve - typename R::Line_2 l = typename R::Construct_bisector_2()(p, q); - Point_2 middle = Construct_midpoint_2()(p, q); - Point_2 temp = middle + l.to_vector(); - if(Orientation_2()(middle, temp, Point_2(x, y)) == ON_POSITIVE_SIDE){ - return Circle_2(Point_2(x, y), radius, CLOCKWISE); - } - - return Circle_2(Point_2(x, y), radius, COUNTERCLOCKWISE); + Circle_2 c = construct_supporting_circle(p, q); + // compute the ending points + std::pair points = find_intersection(c); + + if(Orientation_2()(points.first, points.second, c.center()) == LEFT_TURN) { + return Arc_2(c, points.first, points.second); } - - // Find intersection of an input circle orthogonal to the Poincaré disk - // and the circle representing this disk - - // TODO: sqrt(to_double()?) - std::pair find_intersection(Circle_2& circle) const - { - FT x = circle.center().x(), y = circle.center().y(); - - // axˆ2 + 2bˆx + c = 0; - FT a = x*x + y*y; - FT b = -_unit_circle.squared_radius() * x; - FT c = _unit_circle.squared_radius()*_unit_circle.squared_radius() - _unit_circle.squared_radius()*y*y; - assert(b*b - a*c > 0); - FT D = CGAL::sqrt(to_double(b*b - a*c)); - - FT x1 = (-b - D)/a; - FT x2 = (-b + D)/a; - FT y1 = (_unit_circle.squared_radius() - x1*x)/y; - FT y2 = (_unit_circle.squared_radius() - x2*x)/y; - - return std::make_pair(Point_2(x1, y1), Point_2(x2, y2)); - } - - // Find intersection of an input line orthogonal to the Poincaré disk - // and the circle representing this disk - - // TODO: sqrt(to_double()?) - std::pair find_intersection(Euclidean_line_2& l) const - { - typedef typename R::Vector_2 Vector_2; - Vector_2 v = l.to_vector(); - - // normalize the vector - FT squared_coeff = _unit_circle.squared_radius()/v.squared_length(); - FT coeff = CGAL::sqrt(to_double(squared_coeff)); - - Point_2 p1(coeff*v.x(), coeff*v.y()); - Point_2 p2(-p1.x(), -p1.y()); - return std::make_pair(p1, p2); - } - - private: - const Circle_2 _unit_circle; - }; + return Arc_2(c, points.second, points.first); + } - Construct_hyperbolic_bisector_2 - construct_hyperbolic_bisector_2_object() const - { return Construct_hyperbolic_bisector_2(_unit_circle);} - - Construct_bisector_2 - construct_bisector_2_object() const - {return Construct_bisector_2();} - - Compare_distance_2 - compare_distance_2_object() const - {return Compare_distance_2();} - - Construct_triangle_2 construct_triangle_2_object() const - {return Construct_triangle_2();} - - Construct_direction_2 construct_direction_2_object() const - {return Construct_direction_2();} - - class Construct_ray_2 + private: + // The cirle belongs to the pencil with limit points p and q + Circle_2 construct_supporting_circle(Point_2 p, Point_2 q) const { - public: - Construct_ray_2(Circle_2 c) : - _unit_circle(c) {} + // p, q are zero-circles + // (x, y, xˆ2 + yˆ2 - rˆ2) = alpha*(xp, yp, xpˆ2 + ypˆ2) + (1-alpha)*(xq, yq, xqˆ2 + yqˆ2) + // xˆ2 + yˆ2 - rˆ2 = Rˆ2, where R - is a radius of the given unit circle + FT op = p.x()*p.x() + p.y()*p.y(); + FT oq = q.x()*q.x() + q.y()*q.y(); + FT alpha = (_unit_circle.squared_radius() - oq) / (op - oq); - Segment_2 operator()(Point_2 p, Segment_2 l) const - { - if(typename R::Segment_2* s = boost::get(&l)){ - return operator()(p, *s); - } - if(Arc_2* arc = boost::get(&l)){ - if(arc->get<0>().orientation() == CLOCKWISE){ - arc->get<1>() = p; - return *arc; - } - arc->get<2>() = p; + FT x = alpha*p.x() + (1-alpha)*q.x(); + FT y = alpha*p.y() + (1-alpha)*q.y(); + FT radius = x*x + y*y - _unit_circle.squared_radius(); + + //improve + typename R::Line_2 l = typename R::Construct_bisector_2()(p, q); + Point_2 middle = Construct_midpoint_2()(p, q); + Point_2 temp = middle + l.to_vector(); + if(Orientation_2()(middle, temp, Point_2(x, y)) == ON_POSITIVE_SIDE){ + return Circle_2(Point_2(x, y), radius, CLOCKWISE); + } + + return Circle_2(Point_2(x, y), radius, COUNTERCLOCKWISE); + } + + // Find intersection of an input circle orthogonal to the Poincaré disk + // and the circle representing this disk + + // TODO: sqrt(to_double()?) + std::pair find_intersection(Circle_2& circle) const + { + FT x = circle.center().x(), y = circle.center().y(); + + // axˆ2 + 2bˆx + c = 0; + FT a = x*x + y*y; + FT b = -_unit_circle.squared_radius() * x; + FT c = _unit_circle.squared_radius()*_unit_circle.squared_radius() - _unit_circle.squared_radius()*y*y; + assert(b*b - a*c > 0); + FT D = CGAL::sqrt(to_double(b*b - a*c)); + + FT x1 = (-b - D)/a; + FT x2 = (-b + D)/a; + FT y1 = (_unit_circle.squared_radius() - x1*x)/y; + FT y2 = (_unit_circle.squared_radius() - x2*x)/y; + + return std::make_pair(Point_2(x1, y1), Point_2(x2, y2)); + } + + // Find intersection of an input line orthogonal to the Poincaré disk + // and the circle representing this disk + + // TODO: sqrt(to_double()?) + std::pair find_intersection(Euclidean_line_2& l) const + { + typedef typename R::Vector_2 Vector_2; + Vector_2 v = l.to_vector(); + + // normalize the vector + FT squared_coeff = _unit_circle.squared_radius()/v.squared_length(); + FT coeff = CGAL::sqrt(to_double(squared_coeff)); + + Point_2 p1(coeff*v.x(), coeff*v.y()); + Point_2 p2(-p1.x(), -p1.y()); + return std::make_pair(p1, p2); + } + + private: + const Circle_2 _unit_circle; + }; + + Construct_hyperbolic_bisector_2 + construct_hyperbolic_bisector_2_object() const + { return Construct_hyperbolic_bisector_2(_unit_circle);} + + Construct_bisector_2 + construct_bisector_2_object() const + {return Construct_bisector_2();} + + Compare_distance_2 + compare_distance_2_object() const + {return Compare_distance_2();} + + Construct_triangle_2 construct_triangle_2_object() const + {return Construct_triangle_2();} + + Construct_direction_2 construct_direction_2_object() const + {return Construct_direction_2();} + + class Construct_ray_2 + { + public: + Construct_ray_2(Circle_2 c) : + _unit_circle(c) {} + + Segment_2 operator()(Point_2 p, Segment_2 l) const + { + if(typename R::Segment_2* s = boost::get(&l)){ + return operator()(p, *s); + } + if(Arc_2* arc = boost::get(&l)){ + if(arc->get<0>().orientation() == CLOCKWISE){ + arc->get<1>() = p; return *arc; } - assert(false); - return Segment_2(); + arc->get<2>() = p; + return *arc; } - - Segment_2 operator()(Point_2 p, typename R::Segment_2 s) const - { - return typename R::Segment_2(p, s.target()); - } - - private: - - const Circle_2 _unit_circle; - }; + assert(false); + return Segment_2(); + } - Construct_ray_2 construct_ray_2_object() const - {return Construct_ray_2(_unit_circle);} + Segment_2 operator()(Point_2 p, typename R::Segment_2 s) const + { + return typename R::Segment_2(p, s.target()); + } + + private: + + const Circle_2 _unit_circle; }; + + Construct_ray_2 construct_ray_2_object() const + {return Construct_ray_2(_unit_circle);} + + // For details see the rapport RR-8146 + class Is_hyperbolic + { + public: + + bool operator() (const Point& p0, const Point& p1, const Point& p2) const + { + Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), + p1.x()*p1.x() + p1.y()*p1.y(), + p2.x()*p2.x() + p2.y()*p2.y()); + + Vector_3 v1 = Vector_3(p0.x(), p1.x(), p2.x()); + Vector_3 v2 = Vector_3(p0.y(), p1.y(), p2.y()); + Vector_3 v3 = Vector_3(FT(1), FT(1), FT(1)); + + FT dt0 = CGAL::determinant(v0, v1, v3); + FT dt1 = CGAL::determinant(v0, v2, v3); + FT dt2 = CGAL::determinant(v0 - v3, v1, v2); + + return dt0*dt0 + dt1*dt1 - dt2*dt2 < 0; + } + + bool operator() (const Point& p0, const Point& p1, const Point& p2, int& ind) const + { + if (this->operator()(p0, p1, p2) == false) { + ind = find_non_hyperbolic_edge(p0, p1, p2); + return false; + } + return true; + } + + private: + + // assume the face (p0, p1, p2) is non-hyperbolic + int find_non_hyperbolic_edge(const Point& p0, const Point& p1, const Point& p2) const + { + typedef typename R::Direction_2 Direction_2; + + Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), + p1.x()*p1.x() + p1.y()*p1.y(), + p2.x()*p2.x() + p2.y()*p2.y()); + + Vector_3 v1 = Vector_3(p0.x(), p1.x(), p2.x()); + Vector_3 v2 = Vector_3(p0.y(), p1.y(), p2.y()); + Vector_3 v3 = Vector_3(FT(1), FT(1), FT(1)); + + FT dt0 = CGAL::determinant(v0, 2*v2, -v3); + FT dt1 = CGAL::determinant(2*v1, v0, -v3); + FT dt2 = CGAL::determinant(2*v1, 2*v2, -v3); + + Direction_2 d0(p0.x()*dt2 - dt0, p0.y()*dt2 - dt1); + Direction_2 d1(p1.x()*dt2 - dt0, p1.y()*dt2 - dt1); + Direction_2 d2(p2.x()*dt2 - dt0, p2.y()*dt2 - dt1); + + Direction_2 d(dt0, dt1); + + if(d.counterclockwise_in_between(d0, d1)) { + return 2; + } + + if(d.counterclockwise_in_between(d1, d2)) { + return 0; + } + + return 1; + } + }; + + Is_hyperbolic Is_hyperbolic_object() const + { return Is_hyperbolic(); } +}; - // Take out the code below to some separate file - + #ifdef CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H template <> struct Triangulation_structural_filtering_traits< Triangulation_hyperbolic_traits_2 > { - typedef Tag_true Use_structural_filtering_tag; + typedef Tag_true Use_structural_filtering_tag; }; #endif // CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H @@ -445,7 +513,7 @@ struct Triangulation_structural_filtering_traits< Triangulation_hyperbolic_trait typedef Tag_true Use_structural_filtering_tag; }; #endif // CGAL_EXACT_PREDICATES_INEXACT_CONSTRUCTIONS_KERNEL_H - + } //namespace CGAL #endif // CGAL_TRIANGULATION_HYPERBOLIC_TRAITS_2_H