First results with insert(). A new test has been added, it inserts 100 random points (using the predicate to see whether they are valid or not).

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<Face_handle>::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;

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<Gt, Tds>::Vertex_handle
 Periodic_4_hyperbolic_Delaunay_triangulation_2<Gt, Tds>::
 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();
+  }
+
 }
 
 

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<Kernel> 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 );
+        // Poincare disk (i.e., unit Euclidean disk)
-        Circle_2 DomainA ( CenterA,           Radius2 );
+        Circle_2 Poincare    ( Point(0, 0),       1*1 );
-        Circle_2 DomainBb( rotate(CenterA),   Radius2 );
 
+        // 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 PoincareSide = Poincare.bounded_side(t);
-        CGAL::Bounded_side bs2 = DomainA.bounded_side(t);
+        CGAL::Bounded_side CircASide    = EuclidCircA.bounded_side(t);
-        CGAL::Bounded_side bs3 = DomainBb.bounded_side(t);
+        CGAL::Bounded_side CircBbSide   = EuclidCircBb.bounded_side(t);
 
-        if (  (bs1 != CGAL::ON_UNBOUNDED_SIDE) &&
+        // First off, the point needs to be inside the Poincare disk. if not, there's no hope.
-              (bs2 != CGAL::ON_BOUNDED_SIDE)   &&
+        if ( PoincareSide == CGAL::ON_BOUNDED_SIDE ) {
-              (bs3 != CGAL::ON_BOUNDED_SIDE)      ) {
+          
-          return CGAL::ON_BOUNDED_SIDE;
+          // Inside the Poincare disk, but still outside the fundamental domain
-        } 
+          if ( CircASide  == CGAL::ON_BOUNDED_SIDE || 
-        else if ( (bs1 == CGAL::ON_UNBOUNDED_SIDE) ||
+               CircBbSide == CGAL::ON_BOUNDED_SIDE   ) {
-                  (bs2 == CGAL::ON_BOUNDED_SIDE)   ||
+            return CGAL::ON_UNBOUNDED_SIDE;
-                  (bs3 == CGAL::ON_BOUNDED_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;
 
       }
 

Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_2.h

@@ -320,6 +320,31 @@ public:
 
 
 
+  template<class EdgeIt>
+  Vertex_handle star_hole(const Point& p, EdgeIt edge_begin, EdgeIt edge_end)
+  {
+    std::list<Face_handle> empty_list;
+    return star_hole(p, edge_begin,         edge_end, 
+                        empty_list.begin(), empty_list.end());
+  }
+
+
+  template<class EdgeIt, class FaceIt>
+  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<Vertex_handle> 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) {
+      	if (verbose) {
-            std::cerr << "Orientation problem in face " << cit->get_number() << endl;
+          cit->restore_orientation();
-      		}
+          for (int i=0; i<3; i++) {
-      		error = true;
+            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 <class GT, class TDS>
 typename TDS::Face_handle Periodic_4_hyperbolic_triangulation_2<GT, TDS>::
-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)

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 <class T>
+  	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(); }

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()
   

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<NT>                                     					Kernel;
+typedef Kernel::FT                                                                  FT;
 typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_traits_2<Kernel> 		Traits;
 typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_2<Traits>				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<Vertex_handle, Vertex_handle>                                     Edge;
-typedef unsigned short int                                                          Int;
-typedef CGAL::Hyperbolic_word_4<Int, Traits>                                        Offset;
-typedef std::pair<Edge, bool>                                                       OEdge;
-typedef std::set<OEdge>                                                             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;

Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_insertion.cpp

@@ -0,0 +1,69 @@
+
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/uniform_smallint.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <CGAL/point_generators_2.h>
+
+#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h>
+#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h>
+#include <CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h>
+#include <CGAL/CORE_Expr.h>
+#include <CGAL/Cartesian.h>
+#include <CGAL/determinant.h>
+
+
+typedef CORE::Expr                                              					NT;					
+typedef CGAL::Cartesian<NT>                                     					Kernel;
+typedef Kernel::FT                                                                  FT;
+typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_traits_2<Kernel> 		Traits;
+typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_2<Traits>				Triangulation;
+
+typedef Kernel::Point_2                                                             Point;
+typedef CGAL::Creator_uniform_2<double, Point>                                      Creator;
+typedef std::vector<Point>                                                          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<Point, Creator> 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;
+}