From bd9b60f7cc27245592f4706f3c8b94316cb7bc1a Mon Sep 17 00:00:00 2001 From: Monique Teillaud Date: Thu, 11 Aug 2016 17:15:18 +0200 Subject: [PATCH] The traits class now uses circular kernel. - dual functions in Delaunay class fixed (but code should still be improved) - examples and demos modified accordingly - broken: Hyperbolic_random_points_in_disc_2 and benchmarks (problems with number types) --- Hyperbolic_triangulation_2/TODO | 54 ++- .../test_hyperbolic_triangulation.cpp | 4 +- .../test_hyperbolic_triangulation_dynamic.cpp | 4 +- ...perbolic_Delaunay_triangulation_2_demo.cpp | 17 +- .../Hyperbolic_Dirichlet_region_2_demo.cpp | 15 +- .../hyperbolic_color.cpp | 4 +- ...erbolic_delaunay_triangulation_example.cpp | 45 +- .../Hyperbolic_Delaunay_triangulation_2.h | 63 +-- ...perbolic_Delaunay_triangulation_traits_2.h | 451 ++++++++++-------- .../CGAL/Qt/HyperbolicPainterOstream.h | 90 ++-- ...Item.h => HyperbolicVoronoiGraphicsItem.h} | 16 +- 11 files changed, 420 insertions(+), 343 deletions(-) rename Hyperbolic_triangulation_2/include/CGAL/Qt/{VoronoiGraphicsItem.h => HyperbolicVoronoiGraphicsItem.h} (90%) diff --git a/Hyperbolic_triangulation_2/TODO b/Hyperbolic_triangulation_2/TODO index b73ae36f60f..65e01ddc03a 100644 --- a/Hyperbolic_triangulation_2/TODO +++ b/Hyperbolic_triangulation_2/TODO @@ -1,28 +1,55 @@ ==== new ==== -add operator == for triangulations +add operator == for triangulations (?) -** clean approximations ---> use CK +** traits bisectors -the construction of the (hyperbolic) segment should not be approximated in the traits -the approximation should be done for the demo only +class Construct_hyperbolic_bisector_2 +Hyperbolic_segment_2 operator()(Point_2 p, Point_2 q, Point_2 r) +at the end +the following lines are giving the wrong arc (wrong orientation), aren't they? + CGAL_triangulation_assertion(assign(pair,inters[1])); + if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE ) + return Circular_arc_2( *c_pq, pair.first, a); + return Circular_arc_2( *c_pq, a, pair.first); +but this is what works on the demo... +understand why -similar: clean find_intersection -(using the circular kernel will do) +remove variant for supporting circle or line of bisector +call it only when we know that it is a circle +it will simplyfy the code of Construct_hyperbolic_bisector_2 at least in some cases -fix construct_circumcenter. no approx allowed -(again, CK will do) +test bisectors dual functions in special cases of euclidean line segments -CK, or at least Cartesian with sqrt (?) +** Hyperbolic_random_points_in_disc + +repair +number types have changed (using CK) + +** benchmarks + +idem, repair number types ** doc -** tests +** test suite ** only keep relevant files for the submission ========== demo +--- input from file does not do anything... + +--- clean files +do we need demo/Hyperbolic_triangulation_2/include/* for the basic demo? +todo diff with similar .h in Graphics view (after updating the branch) + +--- clean buttons +remove translations, add circumcircle, etc +and put appropriate window title + +--- fix aretes delaunay presque rectilineaires pas au bon endroit + --- create an adapted Graphicsitem that does not require finite_* stuff. Replace by hyperbolic_* => understand apply_to_range the graphicsitem should also not require Segment_2 or Line_segment_2 (does it?) @@ -60,13 +87,6 @@ CGAL::https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/Unique_sqrt_ex Done. std::sqrt is replaced by CGAL::sqrt everywhere, because it's used only for computation of the Voronoi diagram. ---- use the Circular_kernel? -because Do_intersect is in the Circular_kernel - -Done. CGAL::do_intersect is a solution. -No need for the Circular_kernel. ---> (MT) wrong for constrcutions. Approximations are used in several places. - ========== future --- make Constrained_delaunay_... work with Hyperbolic_triangulation_2 diff --git a/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation.cpp b/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation.cpp index 16f118cf101..6f54dd16b23 100644 --- a/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation.cpp +++ b/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation.cpp @@ -3,12 +3,12 @@ #include // CGAL headers -#include +#include #include #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_circular_kernel_2 K; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 Gt; typedef K::Point_2 Point_2; diff --git a/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation_dynamic.cpp b/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation_dynamic.cpp index 3ed49197f08..4a72c51d76b 100644 --- a/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation_dynamic.cpp +++ b/Hyperbolic_triangulation_2/benchmark/test_hyperbolic_triangulation_dynamic.cpp @@ -3,12 +3,12 @@ #include // CGAL headers -#include +#include #include #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_circular_kernel_2 K; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 Gt; typedef K::Point_2 Point_2; diff --git a/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Delaunay_triangulation_2_demo.cpp b/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Delaunay_triangulation_2_demo.cpp index c942a63ae35..7c0bbc6dc01 100644 --- a/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Delaunay_triangulation_2_demo.cpp +++ b/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Delaunay_triangulation_2_demo.cpp @@ -1,14 +1,12 @@ #include // CGAL headers -#include -#include -#include -#include -// to be deleted +#include +#include +#include + #include -// #include @@ -22,13 +20,13 @@ // GraphicsView items and event filters (input classes) #include "TriangulationCircumcircle.h" - #include "TriangulationMovingPoint.h" #include "TriangulationConflictZone.h" #include "TriangulationRemoveVertex.h" #include "TriangulationPointInputAndConflictZone.h" #include -#include + +#include // for viewportsBbox #include @@ -37,7 +35,8 @@ #include "ui_Hyperbolic_translations_2.h" #include -typedef CGAL::Exact_predicates_exact_constructions_kernel R; +//typedef CGAL::Exact_predicates_exact_constructions_kernel R; +typedef CGAL::Exact_circular_kernel_2 R; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 K; typedef K::Point_2 Point_2; diff --git a/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Dirichlet_region_2_demo.cpp b/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Dirichlet_region_2_demo.cpp index 05248631fe1..939a82e8f5b 100644 --- a/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Dirichlet_region_2_demo.cpp +++ b/Hyperbolic_triangulation_2/demo/Hyperbolic_triangulation_2/Hyperbolic_Dirichlet_region_2_demo.cpp @@ -1,15 +1,12 @@ #include // CGAL headers -#include -#include -#include -#include +#include #include +#include +#include -// to be deleted #include -// #include @@ -31,7 +28,8 @@ #include "TriangulationRemoveVertex.h" #include "TriangulationPointInputAndConflictZone.h" //#include -#include + +#include // store color #include @@ -55,8 +53,7 @@ #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel InR; -typedef CGAL::Exact_predicates_exact_constructions_kernel R; +typedef CGAL::Exact_circular_kernel_2 R; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 K; typedef K::Point_2 Point_2; diff --git a/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_color.cpp b/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_color.cpp index aed81c5131c..c1e75a4e04c 100644 --- a/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_color.cpp +++ b/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_color.cpp @@ -3,7 +3,7 @@ // CGAL headers #include -#include +#include #include #include @@ -11,7 +11,7 @@ #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_circular_kernel_2 K; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2< K > Gt; typedef Gt::Point_2 Point_2; diff --git a/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_delaunay_triangulation_example.cpp b/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_delaunay_triangulation_example.cpp index 6d0c488e273..8d8d689dc9a 100644 --- a/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_delaunay_triangulation_example.cpp +++ b/Hyperbolic_triangulation_2/examples/Hyperbolic_triangulation_2/hyperbolic_delaunay_triangulation_example.cpp @@ -3,8 +3,7 @@ // CGAL headers #include -#include -#include +#include #include #include @@ -13,7 +12,7 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_circular_kernel_2 K; typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 Gt; typedef Gt::Point_2 Point_2; @@ -25,31 +24,31 @@ typedef CGAL::Hyperbolic_Delaunay_triangulation_2 Dt; int main() { CGAL::Timer timer; - // typedef CGAL::Creator_uniform_2 Creator; + typedef CGAL::Creator_uniform_2 Creator; - // CGAL::Random_points_in_disc_2 in_disc(FT(1)); + CGAL::Random_points_in_disc_2 in_disc(FT(1)); - // int n = 10; - // std::cout << "Number of points: " << n << std::endl; + int n = 100; + std::cout << "Number of points: " << n << std::endl; - // std::vector pts(n); - // std::vector::iterator ip; - - // // Generating n random points - // for (int i=0 ; i < n ; i++) { - // pts.at(i) = *in_disc; - // in_disc++; - // } - - std::vector pts; + std::vector pts(n); std::vector::iterator ip; - Point_2 p; - - std::ifstream ifs("input-file"); - while(ifs >> p) { - pts.push_back(p); + + // Generating n random points + for (int i=0 ; i < n ; i++) { + pts.at(i) = *in_disc; + in_disc++; } - std::cout << "number of points " << std::distance(pts.begin(),pts.end()) << std::endl << std::endl; + + // std::vector pts; + // std::vector::iterator ip; + // Point_2 p; + + // std::ifstream ifs("input-file"); + // while(ifs >> p) { + // pts.push_back(p); + // } + // std::cout << "number of points " << std::distance(pts.begin(),pts.end()) << std::endl << std::endl; std::cout << "check for hyperbolic faces during insertion" << std::endl; diff --git a/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h b/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h index e7426ad4cfd..01dba8b3a9a 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h @@ -64,8 +64,8 @@ public: typedef Gt Geom_traits; typedef typename Geom_traits::FT FT; typedef typename Geom_traits::Point_2 Point; - typedef typename Geom_traits::Hyperbolic_segment_2 Segment; - + typedef typename Geom_traits::Circular_arc_point_2 Circular_arc_point; + typedef typename Geom_traits::Hyperbolic_segment_2 Hyperbolic_segment; Hyperbolic_Delaunay_triangulation_2(const Gt& gt = Gt()) : Delaunay_triangulation_2(gt) {} @@ -477,7 +477,7 @@ public: Finite_edges_iterator finite_edges_begin() const { return hyperbolic_edges_begin(); } Finite_edges_iterator finite_edges_end() const { return hyperbolic_edges_end(); } - Point + Circular_arc_point dual(Face_handle f) const { CGAL_triangulation_precondition (!this->is_non_hyperbolic(f)); @@ -486,65 +486,70 @@ public: ( f->vertex(0)->point(), f->vertex(1)->point(), f->vertex(2)->point()); } - Object + Hyperbolic_segment dual(const Edge& e) const { return dual(e.first, e.second); } - Object + Hyperbolic_segment dual(Face_handle f, int i) const { CGAL_triangulation_precondition (!this->is_non_hyperbolic(f,i)); if(this->dimension() == 1) { - const Point& p = f->vertex(cw(i))->point(); - const Point& q = f->vertex(ccw(i))->point(); + Point p = f->vertex(cw(i))->point(); + Point q = f->vertex(ccw(i))->point(); // hyperbolic line - Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - return make_object(line); + Hyperbolic_segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); + return line; } Face_handle n = f->neighbor(i); int in = n->index(f); - + //TODO MT store values of bools to avoid recomputing is-hyperbolic several times + // boths faces are non_hyperbolic, but the incident edge is hyperbolic - if(is_non_hyperbolic(f) && is_non_hyperbolic(n)){ - const Point& p = f->vertex(cw(i))->point(); - const Point& q = f->vertex(ccw(i))->point(); + if( is_non_hyperbolic(f) && is_non_hyperbolic(n) ){ + const Point& p = f->vertex(ccw(i))->point(); + const Point& q = f->vertex(cw(i))->point(); // hyperbolic line - Segment line = + Hyperbolic_segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - return make_object(line); + return line; } - // both faces are finite - if(!is_non_hyperbolic(f) && !is_non_hyperbolic(n)) { - - Segment s = - this->geom_traits().construct_segment_2_object()(dual(f),dual(n)); - - return make_object(s); + // both faces are hyperbolic + if( !is_non_hyperbolic(f) && !is_non_hyperbolic(n) ) { + const Point& p = f->vertex(ccw(i))->point(); + const Point& q = f->vertex(cw(i))->point(); + + Hyperbolic_segment s = + this->geom_traits().construct_hyperbolic_bisector_2_object() + (p,q,f->vertex(i)->point(),n->vertex(in)->point()); + //TODO MT cut edge at dual points !!!! + return s; } // one of the incident faces is non_hyperbolic - Face_handle finite_face = f; + Face_handle hyp_face = f; if(is_non_hyperbolic(f)) { - finite_face = n; + hyp_face = n; i = in; } - const Point& p = finite_face->vertex(cw(i))->point(); - const Point& q = finite_face->vertex(ccw(i))->point(); + const Point& p = hyp_face->vertex(ccw(i))->point(); + const Point& q = hyp_face->vertex(cw(i))->point(); // ToDo: Line or Segment? // hyperbolic line and ray - Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - Segment ray = this->geom_traits().construct_ray_2_object()(dual(finite_face), line); - return make_object(ray); + Hyperbolic_segment ray = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q,hyp_face->vertex(i)->point()); + // TODO MT cut edge at dual point !!! + // Segment ray = this->geom_traits().construct_ray_2_object()(dual(finite_face), line); + return ray; } }; diff --git a/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h b/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h index b15221d91b8..e08ffe7c15e 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h @@ -22,7 +22,7 @@ #ifndef CGAL_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H #define CGAL_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H -#include +#include #include "boost/tuple/tuple.hpp" #include "boost/variant.hpp" @@ -31,19 +31,25 @@ namespace CGAL { template < class R > class Hyperbolic_Delaunay_triangulation_traits_2 : public R + // R is supposed to be a model of CircularKernel2 { public: typedef typename R::FT FT; typedef typename R::Point_2 Point_2; typedef typename R::Circle_2 Circle_2; + typedef typename R::Line_2 Euclidean_line_2; + typedef boost::variant Euclidean_circle_or_line_2; - typedef boost::tuple Arc_2; + typedef typename R::Circular_arc_2 Circular_arc_2; + typedef typename R::Line_arc_2 Line_arc_2; + typedef typename R::Circular_arc_point_2 Circular_arc_point_2; typedef typename R::Segment_2 Euclidean_segment_2; //only used internally here - typedef boost::variant Hyperbolic_segment_2; + typedef boost::variant Hyperbolic_segment_2; typedef typename R::Compare_x_2 Compare_x_2; typedef typename R::Compare_y_2 Compare_y_2; + typedef typename R::Compare_distance_2 Compare_distance_2; typedef typename R::Orientation_2 Orientation_2; typedef typename R::Side_of_oriented_circle_2 Side_of_oriented_circle_2; @@ -57,12 +63,7 @@ public: typedef typename R::Construct_bisector_2 Construct_Euclidean_bisector_2; typedef typename R::Construct_midpoint_2 Construct_Euclidean_midpoint_2; typedef typename R::Compute_squared_distance_2 Compute_squared_Euclidean_distance_2; - typedef typename R::Line_2 Euclidean_line_2; - typedef typename R::Vector_2 Vector_2; - // used by Is_hyperbolic - typedef typename R::Vector_3 Vector_3; - typedef typename R::Point_3 Point_3; - + typedef typename R::Has_on_bounded_side_2 Has_on_bounded_side_2; // MT useless? // typedef Hyperbolic_Delaunay_triangulation_traits_2 Self; // typedef typename R::RT RT; @@ -74,9 +75,8 @@ public: // typedef typename R::Iso_rectangle_2 Iso_rectangle_2; // typedef typename R::Angle_2 Angle_2; - // typedef typename R::Less_x_2 Less_x_2; - // typedef typename R::Less_y_2 Less_y_2; - // typedef typename R::Compare_distance_2 Compare_distance_2; + typedef typename R::Less_x_2 Less_x_2; + typedef typename R::Less_y_2 Less_y_2; // typedef typename R::Construct_triangle_2 Construct_triangle_2; // typedef typename R::Construct_direction_2 Construct_direction_2; @@ -84,15 +84,14 @@ public: class Construct_hyperbolic_segment_2 { - typedef typename CGAL::Regular_triangulation_filtered_traits_2 Regular_geometric_traits_2; + typedef typename CGAL::Regular_triangulation_euclidean_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_hyperbolic_segment_2() - { - } + Construct_hyperbolic_segment_2() + {} Hyperbolic_segment_2 operator()(const Point_2& p, const Point_2& q) const { @@ -106,19 +105,19 @@ public: Weighted_point_2 wo(Point_2(o), FT(1)); // Poincaré circle Bare_point center = Construct_weighted_circumcenter_2()(wp, wo, wq); - FT radius = Compute_squared_Euclidean_distance_2()(p, center); + FT sq_radius = Compute_squared_Euclidean_distance_2()(p, center); - Circle_2 circle( center, radius); + Circle_2 circle(center, sq_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 Circular_arc_2(circle, p, q); } - return Arc_2(circle, q, p); + return Circular_arc_2(circle, q, p); } - }; + }; // end Construct_hyperbolic_segment_2 Construct_hyperbolic_segment_2 construct_hyperbolic_segment_2_object() const @@ -134,38 +133,71 @@ public: { public: - // TODO: improve this function - Point_2 operator()(Point_2 p, Point_2 q, Point_2 r) - { - CGAL_triangulation_assertion_code(Origin oo; Point_2 poo(oo); Circle_2 co(poo,FT(1))); - CGAL_triangulation_assertion(co.bounded_side(p) == ON_BOUNDED_SIDE); - CGAL_triangulation_assertion(co.bounded_side(q) == ON_BOUNDED_SIDE); - CGAL_triangulation_assertion(co.bounded_side(r) == ON_BOUNDED_SIDE); - - Circle_2 circle(p, q, r); - // circle must be inside the unit one - CGAL_triangulation_assertion(do_intersect(co, circle) == false); - + Circular_arc_point_2 operator()(Point_2 p, Point_2 q, Point_2 r) + { Origin o; Point_2 po = Point_2(o); - if (circle.center() == po) - { return po; } + Circle_2 l_inf(po, FT(1)); + + Euclidean_circle_or_line_2 bis_pq = Construct_circle_or_line_supporting_bisector()(p,q); + Euclidean_circle_or_line_2 bis_qr = Construct_circle_or_line_supporting_bisector()(q,r); - 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() + FT(1); // Poincare disc has radius 1 - FT D = b*b - 4*a; - - 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; + if ( Compare_distance_2()(po,p,q) == EQUAL && + Compare_distance_2()(po,p,r) == EQUAL ) + return po; + // now supporting objects cannot both be Euclidean lines + + std::pair pair; + + Euclidean_line_2* l; + Circle_2* c; + + if ( Circle_2* c_pq = boost::get(&bis_pq) ) + { + if ( Circle_2* c_qr = boost::get(&bis_qr) ) + { + typedef typename CK2_Intersection_traits::type Intersection_result; + std::vector< Intersection_result > inters; + intersection(*c_pq, *c_qr, std::back_inserter(inters)); + + CGAL_triangulation_assertion(assign(pair,inters[0])); + if ( pair.second == 1 ) + { + if ( Has_on_bounded_side_2()( l_inf, pair.first ) ) + return pair.first; + + CGAL_triangulation_assertion(assign(pair,inters[1])); + return pair.first; + } + return pair.first; + } + + // here bis_qr is a line + l = boost::get(&bis_qr); + c = c_pq; + } + + // here bis_pq is a line + l = boost::get(&bis_pq); + c = boost::get(&bis_qr); + + typedef typename CK2_Intersection_traits::type Intersection_result; + std::vector< Intersection_result > inters; + intersection(*l, *c, std::back_inserter(inters)); + + CGAL_triangulation_assertion(assign(pair,inters[0])); + if ( pair.second == 1 ) + { + if ( Has_on_bounded_side_2()( l_inf, pair.first ) ) + return pair.first; + + CGAL_triangulation_assertion(assign(pair,inters[1])); + return pair.first; + } + return pair.first; } - - }; + + }; // end Construct_hyperbolic_circumcenter_2 Hyperbolic_Delaunay_triangulation_traits_2() {} @@ -202,120 +234,141 @@ public: class Construct_hyperbolic_bisector_2 { public: - Construct_hyperbolic_bisector_2() - {} + Construct_hyperbolic_bisector_2() + {} + // constructs a hyperbolic line Hyperbolic_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_Euclidean_distance_2 dist = Compute_squared_Euclidean_distance_2(); Origin o; Point_2 po = Point_2(o); - FT dif = dist(po, p) - dist(po, q); - FT eps = 0.0000000001; + Circle_2 l_inf = Circle_2(po,FT(1)); - // Bisector is straight in euclidean sense - if(dif > -eps && dif < eps){ - - // ideally - //if(Compare_distance_2()(origin, p, q) == EQUAL){ - - // TODO: calling R::Construct_bisector - Euclidean_line_2 l = Construct_Euclidean_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 Euclidean_segment_2(points.first, points.second); - } - return Euclidean_segment_2(points.second, points.first); + if ( Compare_distance_2()(po,p,q) == EQUAL ){ + Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p,q); + + if ( Less_y_2()(p,q) ) + { return Line_arc_2( l, l_inf, false, l_inf, true ); } + return Line_arc_2( l, l_inf, true, l_inf, false ); } - Circle_2 c = construct_supporting_circle_of_bisector(p, q); - // compute the ending points - std::pair points = find_intersection(c); + Euclidean_circle_or_line_2 + bis_pq = Construct_circle_or_line_supporting_bisector()(p, q); + Circle_2* c = boost::get(&bis_pq); - if(Orientation_2()(points.first, points.second, c.center()) == LEFT_TURN) { - return Arc_2(c, points.first, points.second); + if ( Less_y_2()(po,c->center()) ) + { return Circular_arc_2( *c, l_inf, true, l_inf, false ); } + else if ( Less_y_2()(c->center(),po) ) + { return Circular_arc_2( *c, l_inf, false, l_inf, true ); } + // the center of the circle is on the x-axis + if ( Less_x_2()(po,c->center()) ) + { return Circular_arc_2( *c, l_inf, true, l_inf, false ); } + return Circular_arc_2( *c, l_inf, false, l_inf, true); + } + + // constructs the hyperbolic bisector of segment [p,q] limited by + // circumcenter(p,q,r) on one side + // and circumcenter(p,s,q) on the other side + Hyperbolic_segment_2 + operator()(Point_2 p, Point_2 q, Point_2 r, Point_2 s) + { + CGAL_triangulation_precondition + ( (Orientation_2()(p,q,r) == ON_POSITIVE_SIDE) + && (Orientation_2()(p,s,q) == ON_POSITIVE_SIDE) ); + CGAL_triangulation_precondition + ( (Side_of_oriented_circle_2()(p,q,r,s) == ON_NEGATIVE_SIDE) + && (Side_of_oriented_circle_2()(p,s,q,r) == ON_NEGATIVE_SIDE) ); + + Origin o; + Point_2 po = Point_2(o); + + // TODO MT this is non-optimal... + // the bisector is already computed here + // and it will be recomputed below + Circular_arc_point_2 a = Construct_hyperbolic_circumcenter_2()(p,q,r); + Circular_arc_point_2 b = Construct_hyperbolic_circumcenter_2()(p,s,q); + + if ( Compare_distance_2()(po, p, q) == EQUAL ){ + Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p,q); + return Line_arc_2(l,a,b); } - return Arc_2(c, points.second, points.first); - } - - private: - // The cirle belongs to the pencil with limit points p and q - Circle_2 construct_supporting_circle_of_bisector(Point_2 p, Point_2 q) const + + Euclidean_circle_or_line_2 + bis_pq = Construct_circle_or_line_supporting_bisector()(p, q); + Circle_2* c_pq = boost::get(&bis_pq); + + if ( Compare_distance_2()(po,p,q) == POSITIVE ) + // then p is inside the supporting circle + { return Circular_arc_2(*c_pq,b,a);} + return Circular_arc_2(*c_pq,a,b); + } + + // constructs the hyperbolic bisector of segment [p,q] + // limited by hyperbolic circumcenter(p,q,r) on one side + // and going to the infinite line on the other side + Hyperbolic_segment_2 operator()(Point_2 p, Point_2 q, Point_2 r) { - // 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 = (FT(1) - oq) / (op - oq); // Poincare disc has radius 1 - - FT x = alpha*p.x() + (1-alpha)*q.x(); - FT y = alpha*p.y() + (1-alpha)*q.y(); - FT radius = x*x + y*y - FT(1); - - //improve - Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p, q); - Point_2 middle = Construct_Euclidean_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); + CGAL_triangulation_precondition + ( Orientation_2()(p,q,r) == POSITIVE ); + + Origin o; + Point_2 po = Point_2(o); + Circle_2 l_inf(po, FT(1)); + + // TODO MT this is non-optimal... + // the bisector is computed (at least) twice + Circular_arc_point_2 a = Construct_hyperbolic_circumcenter_2()(p,q,r); + + if ( Compare_distance_2()(po, p, q) == EQUAL ){ + Euclidean_line_2 bis_pq = Construct_Euclidean_bisector_2()(p,q); + typedef typename + CK2_Intersection_traits::type + Intersection_result; + std::vector< Intersection_result > inters; + intersection(bis_pq, l_inf, std::back_inserter(inters)); + std::pair pair; + + CGAL_triangulation_assertion(assign(pair,inters[0])); + CGAL_triangulation_assertion(pair.second == 1); + if ( Less_y_2()(p,q) ) + return Line_arc_2(bis_pq,a,pair.first); + + CGAL_triangulation_assertion(assign(pair,inters[1])); + CGAL_triangulation_assertion(pair.second == 1); + return Line_arc_2(bis_pq,a,pair.first); } + + Euclidean_circle_or_line_2 + bis_pq = Construct_circle_or_line_supporting_bisector()(p, q); + Circle_2* c_pq = boost::get(&bis_pq); - return Circle_2(Point_2(x, y), radius, COUNTERCLOCKWISE); + Point_2 approx_a(to_double(a.x()),to_double(a.y())); + + typedef typename + CK2_Intersection_traits::type + Intersection_result; + std::vector< Intersection_result > inters; + intersection(*c_pq, l_inf, std::back_inserter(inters)); + std::pair pair; + + CGAL_triangulation_assertion(assign(pair,inters[0])); + CGAL_triangulation_assertion(pair.second == 1); + + Point_2 approx_pinf(to_double(pair.first.x()), to_double(pair.first.y())); + Point_2 approx_c(to_double(c_pq->center().x()), + to_double(c_pq->center().y())); + if ( Orientation_2()(p,q,approx_pinf) == NEGATIVE ) { + if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE ) + return Circular_arc_2( *c_pq, a, pair.first ); + return Circular_arc_2( *c_pq, pair.first, a); + } + + CGAL_triangulation_assertion(assign(pair,inters[1])); + if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE ) + return Circular_arc_2( *c_pq, pair.first, a); + return Circular_arc_2( *c_pq, a, pair.first); } - - // Find intersection of an input circle orthogonal to the Poincare 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; */ - FT b = -x; - FT c = 1-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 = (FT(1) - x1*x)/y; - FT y2 = (FT(1) - x2*x)/y; - - return std::make_pair(Point_2(x1, y1), Point_2(x2, y2)); - } - - // Find intersection of an input line orthogonal to the Poincare 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 = FT(1)/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); - } - - }; + }; // end Construct_hyperbolic_bisector_2 Construct_hyperbolic_bisector_2 construct_hyperbolic_bisector_2_object() const @@ -325,45 +378,13 @@ public: construct_Euclidean_bisector_2_object() const { return Construct_Euclidean_bisector_2(); } - class Construct_ray_2 - { - public: - Construct_ray_2() - {} - - Hyperbolic_segment_2 operator()(Point_2 p, Hyperbolic_segment_2 l) const - { - if(Euclidean_segment_2* s = boost::get(&l)){ - return operator()(p, *s); - } - if(Arc_2* arc = boost::get(&l)){ - if(get<0>(*arc).orientation() == CLOCKWISE){ - get<1>(*arc) = p; - return *arc; - } - get<2>(*arc) = p; - return *arc; - } - assert(false); - return Hyperbolic_segment_2(); - } - - Hyperbolic_segment_2 operator()(Point_2 p, Euclidean_segment_2 s) const - { - return Euclidean_segment_2(p, s.target()); - } - - }; - - Construct_ray_2 - construct_ray_2_object() const - { return Construct_ray_2(); } - // For details see the JoCG paper (5:56-85, 2014) class Is_hyperbolic { public: - + typedef typename R::Vector_3 Vector_3; + typedef typename R::Point_3 Point_3; + bool operator() (const Point_2& p0, const Point_2& p1, const Point_2& p2) const { Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), @@ -425,13 +446,71 @@ public: return 1; } - }; - + }; // end Is_hyperbolic + Is_hyperbolic Is_hyperbolic_object() const { return Is_hyperbolic(); } + + // do not document + // constructs the Euclidean circle or line supporting the hyperbolic + // bisector of two points + class Construct_circle_or_line_supporting_bisector + { + public: + Construct_circle_or_line_supporting_bisector() + {} + + Euclidean_circle_or_line_2 + operator()(Point_2 p, Point_2 q) const + { + Origin o; + Point_2 po = Point_2(o); + typedef typename R::Point_3 Point_3; + + if ( Compare_distance_2()(po,p,q) == EQUAL ) + { return Construct_Euclidean_bisector_2()(p,q); } + + FT dop2 = p.x()*p.x() + p.y()*p.y(); + FT doq2 = q.x()*q.x() + q.y()*q.y(); + Point_3 p3( p.x(), p.y(), dop2 ); + Point_3 q3( q.x(), q.y(), doq2 ); + + // TODO MT improve + + // The cirle belongs to the pencil with limit points p and q + // 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 = 1 (= radius of the Poincare disc) + FT op = p.x()*p.x() + p.y()*p.y(); + FT oq = q.x()*q.x() + q.y()*q.y(); + FT alpha = (FT(1) - oq) / (op - oq); + + FT x = alpha*p.x() + (1-alpha)*q.x(); + FT y = alpha*p.y() + (1-alpha)*q.y(); + FT sq_radius = x*x + y*y - FT(1); + + // TODO MT improve + // ?? orientation should depend on + // Compare_distance(O,p,q) + // so that p always on positive side + // ??? + // CK does not care about orientation, circular arcs are + // considered in CCW order in any case + + Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p, q); + Point_2 middle = Construct_Euclidean_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), sq_radius, CLOCKWISE); } + return Circle_2(Point_2(x, y), sq_radius, COUNTERCLOCKWISE); + } + }; // end Construct_supporting_circle_of_bisector + }; + // Take out the code below to some separate file #ifdef CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H diff --git a/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicPainterOstream.h b/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicPainterOstream.h index 4b3d5f521cb..444b9163e18 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicPainterOstream.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicPainterOstream.h @@ -16,79 +16,65 @@ // $Id: // // -// Author(s) : Mikhail Bogdanov +// Author(s) : Monique Teillaud +// Mikhail Bogdanov #ifndef CGAL_HYPERBOLIC_PAINTER_OSTREAM_H #define CGAL_HYPERBOLIC_PAINTER_OSTREAM_H #include -#include - namespace CGAL{ namespace Qt { template - class PainterOstream > : public PainterOstream { + class PainterOstream > + : public PainterOstream + { public: typedef PainterOstream Base; - typedef PainterOstream > Self; typedef Hyperbolic_Delaunay_triangulation_traits_2 Gt; + typedef PainterOstream Self; + + typedef typename Gt::Hyperbolic_segment_2 Hyperbolic_segment_2; + + typedef typename Gt::Point_2 Point_2; + typedef Line_arc_2 Line_arc_2; + typedef Circular_arc_2 Circular_arc_2; + typedef Circular_arc_point_2 Circular_arc_point_2; - private: - typedef typename Gt::Segment_2 Segment_2; - typedef typename Gt::Line_segment_2 Line_segment_2; - typedef typename Gt::Arc_2 Arc_2; - //typedef typename Gt::Line_2 Line_2; - - typedef typename K::Point_2 Point_2; - typedef typename K::Circle_2 Circle_2; - - public: PainterOstream(QPainter* p, QRectF rect = QRectF()) - : Base(p, rect), qp(p), convert(rect) {} + : Base(p, rect), qp(p), convert(rect) + {} using Base::operator <<; - PainterOstream& operator << (const Segment_2& s) - { - if(const Line_segment_2* line_segment = boost::get(&s)){ - operator << (*line_segment); - return *this; + PainterOstream& operator << (Hyperbolic_segment_2 s) + { + if (const Line_arc_2* seg = boost::get(&s)) { + operator << (*seg); + return *this; + } + + Circular_arc_2* arc = boost::get(&s); + + if (arc->squared_radius() > 10 ) + // due to rounding, the arc drawn does not look like it + // passes through the endpoints + // so we replace the arc by a line segment + { + Point_2 source(to_double(arc->source().x()),to_double(arc->source().y())); + Point_2 target(to_double(arc->target().x()),to_double(arc->target().y())); + const Line_arc_2 seg(source,target); + operator << (seg); + return *this; + } + + operator << (*arc); + return *this; } - if(const Arc_2* arc = boost::get(&s)){ - - const Circle_2& circle = get<0>(*arc); - const Point_2& center = circle.center(); - const Point_2& source = get<1>(*arc); - const Point_2& target = get<2>(*arc); - - if (circle.squared_radius() > 10) { - const Line_segment_2 seg(source,target); - operator << (seg); - return *this; - } - - double asource = std::atan2( -to_double(source.y() - center.y()), - to_double(source.x() - center.x())); - double atarget = std::atan2( -to_double(target.y() - center.y()), - to_double(target.x() - center.x())); - - std::swap(asource, atarget); - double aspan = atarget - asource; - - if(aspan < 0.) - aspan += 2 * CGAL_PI; - - const double coeff = 180*16/CGAL_PI; - this->qp->drawArc(this->convert(circle.bbox()), - (int)(asource * coeff), - (int)(aspan * coeff)); - } - return *this; - } private: // ToDo: These objects must be deleted diff --git a/Hyperbolic_triangulation_2/include/CGAL/Qt/VoronoiGraphicsItem.h b/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicVoronoiGraphicsItem.h similarity index 90% rename from Hyperbolic_triangulation_2/include/CGAL/Qt/VoronoiGraphicsItem.h rename to Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicVoronoiGraphicsItem.h index 07b199fca17..e863ab96ce3 100644 --- a/Hyperbolic_triangulation_2/include/CGAL/Qt/VoronoiGraphicsItem.h +++ b/Hyperbolic_triangulation_2/include/CGAL/Qt/HyperbolicVoronoiGraphicsItem.h @@ -108,20 +108,12 @@ VoronoiGraphicsItem
::paint(QPainter *painter, const QStyleOptionGraphicsItem for(typename DT::Finite_edges_iterator eit = dt->finite_edges_begin(); eit != dt->finite_edges_end(); - eit++){ - CGAL::Object o = dt->dual(*eit); - typename DT::Segment s; - typename DT::Geom_traits::Ray_2 r; - typename DT::Geom_traits::Line_2 l; - if(CGAL::assign(s,o)){ + eit++) + { + typename DT::Hyperbolic_segment s = dt->dual(*eit); pos << s; - } else if(CGAL::assign(r,o)) { - pos << r; - }else if(CGAL::assign(l,o)) { - pos << l; } - } - + // delete painter->setPen(old); }