diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h index 50345f5c747..aef95e5d509 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h @@ -31,6 +31,16 @@ #include #include +#include +#include + +#include +#include +#include +#include +#include + + namespace CGAL { @@ -113,6 +123,28 @@ namespace CGAL { typedef Triple< Vertex_handle, Vertex_handle, Vertex_handle > Vertex_triple; + + class Dummy_point { + private: + Point _pt; + bool _is_inserted; + Vertex_handle _vh; + + public: + + Dummy_point(FT x, FT y): _pt(x,y), _is_inserted(true) {} + Dummy_point(Point p): _pt(p), _is_inserted(true) {} + + Point operator()() { return _pt; } + bool is_inserted() { return _is_inserted; } + Vertex_handle vertex() { return _vh; } + void set_inserted(bool val) { _is_inserted = val; } + void set(Point val) { _pt = val; } + void set_vertex(Vertex_handle v) { _vh = v; } + }; + + std::vector dummy_points; + public: typedef Periodic_4_hyperbolic_triangulation_triangle_iterator_2 Periodic_triangle_iterator; @@ -127,7 +159,7 @@ namespace CGAL { protected: mutable std::vector v_offsets; - int f_cnt, v_cnt; + int f_cnt, v_cnt, n_dpt; private: @@ -135,12 +167,12 @@ namespace CGAL { public: Periodic_4_hyperbolic_Delaunay_triangulation_2(Geometric_traits gt) : - Periodic_4_hyperbolic_triangulation_2(gt) { } + Periodic_4_hyperbolic_triangulation_2(gt) { n_dpt = 14; } Periodic_4_hyperbolic_Delaunay_triangulation_2( const Circle_2 domain = Circle_2(Point_2(FT(0),FT(0)), FT(1*1)), const Geometric_traits > = Geometric_traits() ) : - Periodic_4_hyperbolic_triangulation_2(domain, gt) { } + Periodic_4_hyperbolic_triangulation_2(domain, gt) { n_dpt = 14; } Periodic_4_hyperbolic_Delaunay_triangulation_2(const Periodic_4_hyperbolic_Delaunay_triangulation_2& tr) : Periodic_4_hyperbolic_triangulation_2(tr) { } @@ -156,6 +188,8 @@ namespace CGAL { void remove(Vertex_handle v); + int number_of_dummy_points() { return n_dpt; } + bool _side_of_octagon( const Face_handle& fh, const Offset& offset) const { int cnt = 0; typename GT::Side_of_fundamental_octagon side; @@ -163,7 +197,7 @@ namespace CGAL { Offset o = offset.inverse().append(fh->vertex(j)->get_offset()); Point p = o.apply( fh->vertex(j)->point() ); if ( side(p) == CGAL::ON_UNBOUNDED_SIDE ) { - if ( p.y() + (CGAL_PI / FT(8))*p.x() > 0 ) { + if ( p.y() + tan(CGAL_PI / FT(8))*p.x() > 0 ) { cnt++; } else { } @@ -175,6 +209,70 @@ namespace CGAL { }; // class Periodic_4_hyperbolic_Delaunay_triangulation_2 +template +typename Gt::FT dist(typename Gt::Point_2 p, typename Gt::Point_2 q) { + typename Gt::FT r = (p.x() - q.x())*(p.x() - q.x()) + (p.y() - q.y())*(p.y() - q.y()); + return sqrt(r); +} + + +template +typename Gt::FT +hyperbolic_diameter(typename Gt::Circle_2 c) { + typedef typename Gt::FT FT; + typedef typename Gt::Kernel K; + typedef typename Gt::Point_2 Point; + typedef typename Gt::Line_2 Line; + typedef typename Gt::Circle_2 Circle; + + + typedef CGAL::Exact_circular_kernel_2 CircK; + typedef CGAL::Point_2 Pt2; + typedef CGAL::Circle_2 Circ2; + typedef CGAL::Line_2 Line2; + typedef std::pair, unsigned> IsectOutput; + + typedef CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple, + CGAL::cpp0x::tuple< std::back_insert_iterator > > > Dispatcher; + + std::vector res0, res1; + Dispatcher disp1 = CGAL::dispatch_output( std::back_inserter(res1) ); + Dispatcher disp0 = CGAL::dispatch_output( std::back_inserter(res0) ); + + Pt2 ct(to_double(c.center().x()), to_double(c.center().y())); + Pt2 p0(0, 0); + double r = to_double(c.squared_radius()); + + Line2 ell(p0, ct); + Circ2 c2(ct, r); + Circ2 c0(p0, 1); + + if (ell.is_degenerate()) { + //cout << "\tThis is degenerate case!" << endl; + return 5.; + } else { + CGAL::intersection(c0, ell, disp0); + CGAL::intersection(c2, ell, disp1); + } + Point a(to_double(res0[0].first.x()), to_double(res0[0].first.y())); + Point b(to_double(res0[1].first.x()), to_double(res0[1].first.y())); + + Point p(to_double(res1[0].first.x()), to_double(res1[0].first.y())); + Point q(to_double(res1[1].first.x()), to_double(res1[1].first.y())); + + FT aq = dist(a, q); + FT pb = dist(p, b); + FT ap = dist(a, p); + FT qb = dist(q, b); + + //cout << "aq = " << aq << ", pb = " << pb << " | ap = " << ap << ", qb = " << qb << endl; + + double hyperdist = fabs(log(to_double((aq*pb)/(ap*qb)))); + + return hyperdist; +} + template bool @@ -187,7 +285,10 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2& dt, std::map(lc)/FT(2); // This is the exact value of the limit (1/4 of the squared systole) //FT lim = FT( pow(acosh(1. + sqrt(2.)),2.)/4. ); @@ -207,7 +308,7 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2& dt, std::map& dt, std::mapvertex(0)->point(), fit->vertex(1)->point(), fit->vertex(2)->point()); - if (max_sq_radius < c.squared_radius()) { - max_sq_radius = c.squared_radius(); + FT diam = hyperbolic_diameter(c); + if (max_diam < diam) { + max_diam = diam; } } } - - if (max_sq_radius < lim) { + if (max_diam < lim) { return true; } else { return false; @@ -287,6 +388,26 @@ insert(const Point &p, Face_handle start) { CGAL_triangulation_assertion(this->is_valid()); + for (int i = 0; i < dummy_points.size(); i++) { + if (dummy_points[i].is_inserted()) { + typedef Delaunay_triangulation_2 Delaunay; + Delaunay dt; + std::map vmap; + + if (is_removable(dummy_points[i].vertex(), dt, vmap)) { + //cout << "Removing dummy point " << i << endl; + remove(dummy_points[i].vertex()); + dummy_points[i].set_inserted(false); + } + } + } + + n_dpt = 0; + for (int i = 0; i < dummy_points.size(); i++) { + if (dummy_points[i].is_inserted()) + n_dpt++; + } + return v; } @@ -426,7 +547,7 @@ remove(Vertex_handle v) { CGAL_triangulation_assertion(this->is_valid()); } else { // is not removable - cout << "Vertex " << v->idx() << " cannot be removed!" << endl; + //cout << "Vertex " << v->idx() << " cannot be removed!" << endl; } } diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h index 26406442396..978a3389197 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h @@ -65,32 +65,30 @@ namespace CGAL { f_cnt = fcount; v_cnt = vcount; - std::vector pts; - if (rational) { // Push back the origin - pts.push_back(Point( FT(0), FT(0) )); + dummy_points.push_back(Dummy_point( FT(0), FT(0) )); // Compute rational approximation of internal points - pts.push_back( Point( FT( 1)/FT(2), FT(-4)/FT(19) ) ); // 0 - pts.push_back( Point( FT( 1)/FT(2), FT( 4)/FT(19) ) ); // 1 - pts.push_back( Point( FT( 4)/FT(19), FT( 1)/FT(2) ) ); // 2 - pts.push_back( Point( FT(-4)/FT(19), FT( 1)/FT(2) ) ); // 3 - pts.push_back( Point( FT(-1)/FT(2), FT( 4)/FT(19) ) ); // 4 - pts.push_back( Point( FT(-1)/FT(2), FT(-4)/FT(19) ) ); // 5 - pts.push_back( Point( FT(-4)/FT(19), FT(-1)/FT(2) ) ); // 6 - pts.push_back( Point( FT( 4)/FT(19), FT(-1)/FT(2) ) ); // 7 + dummy_points.push_back( Dummy_point( FT( 1)/FT(2), FT(-4)/FT(19) ) ); // 0 + dummy_points.push_back( Dummy_point( FT( 1)/FT(2), FT( 4)/FT(19) ) ); // 1 + dummy_points.push_back( Dummy_point( FT( 4)/FT(19), FT( 1)/FT(2) ) ); // 2 + dummy_points.push_back( Dummy_point( FT(-4)/FT(19), FT( 1)/FT(2) ) ); // 3 + dummy_points.push_back( Dummy_point( FT(-1)/FT(2), FT( 4)/FT(19) ) ); // 4 + dummy_points.push_back( Dummy_point( FT(-1)/FT(2), FT(-4)/FT(19) ) ); // 5 + dummy_points.push_back( Dummy_point( FT(-4)/FT(19), FT(-1)/FT(2) ) ); // 6 + dummy_points.push_back( Dummy_point( FT( 4)/FT(19), FT(-1)/FT(2) ) ); // 7 // Compute rational approximations of the midpoints - pts.push_back( Point( FT(-9)/FT(14), FT(0) ) ); // 4 - pts.push_back( Point( FT(-5)/FT(11), FT(-5)/FT(11) ) ); // 5 - pts.push_back( Point( FT(0), FT(-9)/FT(14) ) ); // 6 - pts.push_back( Point( FT(5)/FT(11), FT(-5)/FT(11) ) ); // 7 + dummy_points.push_back( Dummy_point( FT(-9)/FT(14), FT(0) ) ); // 4 + dummy_points.push_back( Dummy_point( FT(-5)/FT(11), FT(-5)/FT(11) ) ); // 5 + dummy_points.push_back( Dummy_point( FT(0), FT(-9)/FT(14) ) ); // 6 + dummy_points.push_back( Dummy_point( FT(5)/FT(11), FT(-5)/FT(11) ) ); // 7 // The vertex v_0 - pts.push_back( Point( FT(97)/FT(125), FT(-26)/FT(81) ) ); // 0 + dummy_points.push_back( Dummy_point( FT(97)/FT(125), FT(-26)/FT(81) ) ); // 0 } else { // Algebraic dummy points @@ -99,7 +97,7 @@ namespace CGAL { FT F1 = FT(1); FT F2 = FT(2); - pts.push_back(Point( FT(0), FT(0) )); // origin + dummy_points.push_back(Dummy_point( FT(0), FT(0) )); // origin FT i1(CGAL::sqrt(FT(2) - sqrt(FT(2)))); FT i2(FT(2)*i1); @@ -107,17 +105,17 @@ namespace CGAL { FT i4(-1 + i3); Point a0(FT( CGAL::sqrt( i4 )), FT(0)); // internal point a0 = rot8(a0); - pts.push_back(a0); + dummy_points.push_back(Dummy_point(a0)); for (int k = 1; k < 8; k++) { a0 = rotate(a0); - pts.push_back(a0); + dummy_points.push_back(Dummy_point(a0)); } - pts.push_back(Point(FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)), F0)); // μ(s_0) - pts.push_back(Point(FT( -CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2) ), FT(-CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_1) - pts.push_back(Point(F0, FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)))); // μ(s_2) - pts.push_back(Point(FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)), -FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_3) - pts.push_back(Point(FT( CGAL::sqrt(CGAL::sqrt(F2) + F1)/ F2 ), - FT( CGAL::sqrt(CGAL::sqrt(F2) - F1)/ F2) )); // v_0 + dummy_points.push_back(Dummy_point(FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)), F0)); // μ(s_0) + dummy_points.push_back(Dummy_point(FT( -CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2) ), FT(-CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_1) + dummy_points.push_back(Dummy_point(F0, FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)))); // μ(s_2) + dummy_points.push_back(Dummy_point(FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)), -FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_3) + dummy_points.push_back(Dummy_point(FT( CGAL::sqrt(CGAL::sqrt(F2) + F1)/ F2 ), - FT( CGAL::sqrt(CGAL::sqrt(F2) - F1)/ F2) )); // v_0 } @@ -199,7 +197,8 @@ namespace CGAL { Vertex_handle vertices[14]; for (int i = 0; i < vcount; i++) { vertices[i] = tds().create_vertex(); - vertices[i]->set_point(pts[i]); + vertices[i]->set_point(dummy_points[i]()); + dummy_points[i].set_vertex(vertices[i]); vertices[i]->set_idx(i); } @@ -212,11 +211,7 @@ namespace CGAL { faces[i]->set_number(i); } - - cout << "Setting dimension!" << endl; tds().set_dimension(2); - cout << "Dimension is now: " << tds().dimension() << endl; - for (int i = 0; i < 4; i++) { int fn = i; diff --git a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt index ce3dcc02d05..dcca31bfd95 100644 --- a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt +++ b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt @@ -16,10 +16,11 @@ if ( CGAL_FOUND ) include_directories (BEFORE "../../include") - #create_single_source_cgal_program( "test_p4ht2_dummy_points.cpp" ) - #create_single_source_cgal_program( "test_p4ht2_locate.cpp" ) - #create_single_source_cgal_program( "test_p4ht2_insertion.cpp" ) + create_single_source_cgal_program( "test_p4ht2_dummy_points.cpp" ) + create_single_source_cgal_program( "test_p4ht2_locate.cpp" ) + create_single_source_cgal_program( "test_p4ht2_insertion.cpp" ) create_single_source_cgal_program( "test_p4ht2_removal.cpp" ) + create_single_source_cgal_program( "test_p4ht2_remove_dummy_points.cpp" ) else() diff --git a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_remove_dummy_points.cpp b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_remove_dummy_points.cpp new file mode 100644 index 00000000000..b0f8edc2340 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_remove_dummy_points.cpp @@ -0,0 +1,73 @@ + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +typedef CORE::Expr NT; +typedef CGAL::Cartesian Kernel; +typedef Kernel::FT FT; +typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_traits_2 Traits; +typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_2 Triangulation; +typedef Triangulation::Face_iterator Face_iterator; +typedef Triangulation::Vertex_handle Vertex_handle; +typedef Triangulation::Point Point; +typedef Traits::Side_of_fundamental_octagon Side_of_fundamental_octagon; +typedef CGAL::Creator_uniform_2 Creator; + +int main(int argc, char** argv) { + + if (argc < 2) { + cout << "usage: " << argv[0] << " [number_of_iterations]" << endl; + return -1; + } + + int iter = atoi(argv[1]); + CGAL::Random_points_in_disc_2 g( 1.0 ); + Side_of_fundamental_octagon pred; + + int min = 1000; + int max = -1; + double mean = 0.0; + for (int j = 0; j < iter; j++) { + cout << "Iteration " << (j+1) << "/" << iter << "..." << endl; + + Triangulation tr; + tr.insert_dummy_points(); + assert(tr.is_valid(true)); + + int cnt = 0; + do { + Point pt = *g; + ++g; + if (pred(pt) != CGAL::ON_UNBOUNDED_SIDE) { + tr.insert(pt); + cnt++; + } + } while (tr.number_of_dummy_points() > 0); + assert(tr.is_valid()); + + if (cnt > max) + max = cnt; + if (cnt < min) + min = cnt; + mean += cnt; + } + mean /= (double)iter; + + cout << "Finished " << iter << " iterations!" << endl; + cout << "Minimum number of points inserted: " << min << endl; + cout << "Maximum number of points inserted: " << max << endl; + cout << "Average number of points inserted: " << mean << endl; + return 0; +} \ No newline at end of file