diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/triangulation.cpp b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/triangulation.cpp index a074e32ae8c..d07a52fbcd3 100644 --- a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/triangulation.cpp +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/triangulation.cpp @@ -24,6 +24,7 @@ typedef Triangulation::Point_2 Point; typedef Triangulation::Face_handle Face_handle; typedef Triangulation::Vertex_handle Vertex_handle; typedef Triangulation::Locate_type Locate_type; +typedef Triangulation::Edge Edge; typedef Traits::FT FT; typedef std::set::iterator face_iterator; @@ -42,7 +43,7 @@ int main(void) { // This is the barycenter of face 10 Point query = Point(FT(-0.59841), FT(-0.15901)); - Face_handle fh = tr.locate( query, lt, li ); + Face_handle fh = tr.euclidean_visibility_locate( query, lt, li ); cout << "Point " << query << " is located in face " << fh->get_number(); cout << (lt == Triangulation::EDGE ? " [EDGE]" : (lt == Triangulation::VERTEX ? " [VERTEX]" : " [FACE]")) << endl; @@ -54,6 +55,13 @@ int main(void) { } cout << endl; + cout << "Number of faces before: " << tr.number_of_faces() << endl; + + tr.insert(query, fh); + + cout << "Number of faces after: " << tr.number_of_faces() << endl; + +/* cout << endl; faces_in_conflict.clear(); @@ -99,6 +107,8 @@ int main(void) { } cout << endl; +*/ + cout << endl << "Triangulation is valid: " << (tr.is_valid() ? "TRUE" : "FALSE") << endl; return 0; 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 6155a8db1e7..74ce1fa1f8c 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 @@ -164,7 +164,6 @@ public: InputIterator last, bool is_large_point_set = true); - }; // class Periodic_4_hyperbolic_Delaunay_triangulation_2 @@ -174,8 +173,18 @@ typename Periodic_4_hyperbolic_Delaunay_triangulation_2::Vertex_handle Periodic_4_hyperbolic_Delaunay_triangulation_2:: insert(const Point &p, Face_handle start) { - std::cout << "Inseritng now! Yeah!" << std::endl; - return Vertex_handle(); + + typedef typename Gt::Side_of_fundamental_octagon Side_of_fundamental_octagon; + + Side_of_fundamental_octagon check = Side_of_fundamental_octagon(); + CGAL::Bounded_side side = check(p); + + if (side != CGAL::ON_UNBOUNDED_SIDE) { + return this->insert_in_face(p, start); + } else { + return Vertex_handle(); + } + } diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h index 311ee908a05..b7cdad77adf 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h @@ -294,41 +294,69 @@ public: CGAL::Bounded_side operator()(Point_2 p) { + // Rotation by pi/4 CGAL::Aff_transformation_2 rotate(CGAL::ROTATION, std::sqrt(0.5), std::sqrt(0.5)); + + // The center of the Euclidean circle corresponding to the side s_1 (east) Point_2 CenterA ( FT( sqrt((sqrt(2.) + 1.) / 2.) ), FT(0.) ); + + // The squared radius of the Eucliden circle corresponding to the side s_1 FT Radius2 ( (sqrt(2.) - 1.) / 2. ); - Circle_2 Poincare( Point(0, 0), 1*1 ); - Circle_2 DomainA ( CenterA, Radius2 ); - Circle_2 DomainBb( rotate(CenterA), Radius2 ); + // Poincare disk (i.e., unit Euclidean disk) + Circle_2 Poincare ( Point(0, 0), 1*1 ); + // Euclidean circle corresponding to s_1 + Circle_2 EuclidCircA ( CenterA, Radius2 ); + + // Euclidean circle corresponding to s_2 (just rotate the center, radius is the same) + Circle_2 EuclidCircBb( rotate(CenterA), Radius2 ); + + // This transformation brings the point in the first quadrant (positive x, positive y) FT x(FT(p.x()) > FT(0) ? p.x() : -p.x()); FT y(FT(p.y()) > FT(0) ? p.y() : -p.y()); + // This brings the point in the first octant (positive x and y < x) if (y > x) { FT tmp = x; x = y; y = tmp; } + // This tells us whether the point is in the side of the open boundary + bool on_open_side = ( ( p.y() + (CGAL_PI / 8.) * p.x() ) < 0.0 ); + Point t(x, y); - CGAL::Bounded_side bs1 = Poincare.bounded_side(t); - CGAL::Bounded_side bs2 = DomainA.bounded_side(t); - CGAL::Bounded_side bs3 = DomainBb.bounded_side(t); + CGAL::Bounded_side PoincareSide = Poincare.bounded_side(t); + CGAL::Bounded_side CircASide = EuclidCircA.bounded_side(t); + CGAL::Bounded_side CircBbSide = EuclidCircBb.bounded_side(t); - if ( (bs1 != CGAL::ON_UNBOUNDED_SIDE) && - (bs2 != CGAL::ON_BOUNDED_SIDE) && - (bs3 != CGAL::ON_BOUNDED_SIDE) ) { - return CGAL::ON_BOUNDED_SIDE; - } - else if ( (bs1 == CGAL::ON_UNBOUNDED_SIDE) || - (bs2 == CGAL::ON_BOUNDED_SIDE) || - (bs3 == CGAL::ON_BOUNDED_SIDE) ) { + // First off, the point needs to be inside the Poincare disk. if not, there's no hope. + if ( PoincareSide == CGAL::ON_BOUNDED_SIDE ) { + + // Inside the Poincare disk, but still outside the fundamental domain + if ( CircASide == CGAL::ON_BOUNDED_SIDE || + CircBbSide == CGAL::ON_BOUNDED_SIDE ) { + return CGAL::ON_UNBOUNDED_SIDE; + } + + // Inside the Poincare disk and inside the fundamental domain + if ( CircASide == CGAL::ON_UNBOUNDED_SIDE && + CircBbSide == CGAL::ON_UNBOUNDED_SIDE ) { + return CGAL::ON_BOUNDED_SIDE; + } + + // This is boundary, but we only consider the upper half. The lower half means outside. + if (on_open_side) { + return CGAL::ON_UNBOUNDED_SIDE; + } else { + return CGAL::ON_BOUNDARY; + } + + } else { return CGAL::ON_UNBOUNDED_SIDE; } - - return CGAL::ON_BOUNDARY; } diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h index 87e6a06aa2a..64ce33c0db2 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h @@ -320,6 +320,31 @@ public: + template + Vertex_handle star_hole(const Point& p, EdgeIt edge_begin, EdgeIt edge_end) + { + std::list empty_list; + return star_hole(p, edge_begin, edge_end, + empty_list.begin(), empty_list.end()); + } + + + template + Vertex_handle star_hole(const Point& p, EdgeIt edge_begin, EdgeIt edge_end, + FaceIt face_begin, FaceIt face_end) + { + Vertex_handle v = _tds.star_hole(edge_begin, edge_end, face_begin, face_end); + v->set_point(p); + return v; + } + + + Vertex_handle insert_in_face(const Point& p, Face_handle fh) { + Vertex_handle vh = _tds.insert_in_face(fh); + vh->set_point(p); + return vh; + } + bool is_infinite(Vertex_handle v) const { @@ -854,7 +879,7 @@ public: public: std::vector insert_dummy_points(); - Face_handle locate(const Point& p, Locate_type& lt, int& li, Face_handle f = Face_handle()) const; + Face_handle euclidean_visibility_locate(const Point& p, Locate_type& lt, int& li, Face_handle f = Face_handle()) const; protected: // COMMON INSERTION for DELAUNAY and REGULAR TRIANGULATION @@ -1436,10 +1461,18 @@ is_valid(bool verbose, int level) const { } if (orientation( *p[0], *p[1], *p[2], off[0], off[1], off[2] ) != POSITIVE) { - if (verbose) { - std::cerr << "Orientation problem in face " << cit->get_number() << endl; - } - error = true; + if (verbose) { + cit->restore_orientation(); + for (int i=0; i<3; i++) { + p[i] = &cit->vertex(i)->point(); + off[i] = cit->offset(i); + } + if (orientation( *p[0], *p[1], *p[2], + off[0], off[1], off[2] ) != POSITIVE) { + std::cerr << "Orientation problem in face " << cit->get_number() << ", adjustment attempt failed!" << endl; + error = true; + } + } } } @@ -1515,7 +1548,7 @@ find_in_conflict( const Point& p, template typename TDS::Face_handle Periodic_4_hyperbolic_triangulation_2:: -locate(const Point& p, Locate_type& lt, int& li, Face_handle f) const +euclidean_visibility_locate(const Point& p, Locate_type& lt, int& li, Face_handle f) const { // Random generator (used to introduce a small perturbation to the choice of the starting vertex each time) diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_ds_face_base_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_ds_face_base_2.h index b9bf03fd941..c44c13b48bf 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_ds_face_base_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_ds_face_base_2.h @@ -217,6 +217,10 @@ public: V[0] = V[1] = V[2] = Vertex_handle(); } + void set_vertex(int k, Vertex_handle& vh) { + V[k] = vh; + } + void set_vertices( const Vertex_handle& v0, const Vertex_handle v1, const Vertex_handle& v2) @@ -239,6 +243,10 @@ public: N[2] = n2; } + void set_neighbor(int k, Face_handle& nfh) { + N[k] = nfh; + } + // CHECKING // the following trivial is_valid allows @@ -249,6 +257,38 @@ public: return true; } + int dimension() { + int d = 2; + for (int i = 0; i < 3; i++) + if (V[0] == Vertex_handle()) + d--; + + if (d < 0) + d = 0; + + return d; + } + + + void restore_orientation() { + // N(eighbors), V(ertices), o(ffsets), no (neighbor offsets) + + swap(N[1], N[2]); + swap(V[1], V[2]); + swap(o[1], o[2]); + swap(no[1], no[2]); + + } + + template + void swap(T& a, T& b) { + T tmp; + tmp = a; + a = b; + b = tmp; + } + + // For use by Compact_container. void * for_compact_container() const { return N[0].for_compact_container(); } void * & for_compact_container() { return N[0].for_compact_container(); } 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 79e48e09fa0..1bd1ce4308f 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 @@ -17,6 +17,7 @@ if ( CGAL_FOUND ) include_directories (BEFORE "../../include") create_single_source_cgal_program( "test_p4ht2_dummy_points.cpp" ) + create_single_source_cgal_program( "test_p4ht2_insertion.cpp" ) else() diff --git a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_dummy_points.cpp b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_dummy_points.cpp index 81db10c391e..8c90a4942c5 100644 --- a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_dummy_points.cpp +++ b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_dummy_points.cpp @@ -15,20 +15,10 @@ 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::Point_2 Point_2; -typedef Triangulation::Face_handle Face_handle; -typedef Triangulation::Vertex_handle Vertex_handle; -typedef Triangulation::Locate_type Locate_type; - -typedef std::pair Edge; -typedef unsigned short int Int; -typedef CGAL::Hyperbolic_word_4 Offset; -typedef std::pair OEdge; -typedef std::set OEdgeSet; -typedef OEdgeSet::iterator OEI; int ccw(int i) { return (i+1)%3; @@ -39,7 +29,7 @@ int main(void) { Triangulation tr; tr.insert_dummy_points(); - cout << "Triangulation successfully initialized with dummy points!" << endl; + cout << "Triangulation successfully initialized with dummy points!" << endl << "---------------------------------------------" << endl; cout << "Number of vertices: " << tr.number_of_vertices() << endl; cout << "Number of faces: " << tr.number_of_faces() << endl; cout << "Number of edges: " << tr.number_of_edges() << endl; diff --git a/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_insertion.cpp b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_insertion.cpp new file mode 100644 index 00000000000..84d83d7977f --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_insertion.cpp @@ -0,0 +1,69 @@ + + +#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 Kernel::Point_2 Point; +typedef CGAL::Creator_uniform_2 Creator; +typedef std::vector Vector; +typedef Triangulation::Face_handle Face_handle; +typedef Triangulation::Vertex_handle Vertex_handle; +typedef Triangulation::Locate_type Locate_type; + +int ccw(int i) { + return (i+1)%3; +} + +int main(void) { + + Triangulation tr; + tr.insert_dummy_points(); + + int N = 100; + Vector pts; + pts.reserve(N); + CGAL::Random_points_in_disc_2 g( 1.0 ); + CGAL::cpp11::copy_n( g, N, std::back_inserter(pts)); + + Locate_type lt; + int li; + + int bad = 0; + for (int i = 0; i < N; i++) { + Face_handle fh = tr.euclidean_visibility_locate( pts[i], lt, li ); + Vertex_handle vh = tr.insert(pts[i], fh); + if (vh == Vertex_handle()) + bad++; + } + + cout << "Tried to insert " << N << " random points, " << bad << " were rejected." << endl; + cout << "Number of vertices: " << tr.number_of_vertices() << endl; + cout << "Number of faces: " << tr.number_of_faces() << endl; + cout << "Number of edges: " << tr.number_of_edges() << endl; + cout << "Expected edges (by Euler relation): " << tr.number_of_vertices() + tr.number_of_faces() + 2 << endl; + + cout << "Triangulation is valid: " << (tr.is_valid(true) ? "YES" : "NO") << endl; + + assert(tr.is_valid()); + + + + return 0; +} \ No newline at end of file