Triangulation on hyperbolic surface 2: Delaunay triangulations and epsilon-nets

Documentation/doc/biblio/cgal_manual.bib

@@ -714,6 +714,23 @@ @article{cgal:dljjaw-sisp-15
   publisher={ACM New York, NY, USA}
 }
 
+@InProceedings{cgal:dlpt-naihs-25,
+  author =	{Despr\'{e}, Vincent and Lanuel, Camille and Pouget, Marc and Teillaud, Monique},
+  title =	{{\epsilon-Net Algorithm Implementation on Hyperbolic Surfaces}},
+  booktitle =	{33rd Annual European Symposium on Algorithms (ESA 2025)},
+  pages =	{61:1--61:18},
+  series =	{Leibniz International Proceedings in Informatics (LIPIcs)},
+  ISBN =	{978-3-95977-395-9},
+  ISSN =	{1868-8969},
+  year =	{2025},
+  volume =	{351},
+  editor =	{Benoit, Anne and Kaplan, Haim and Wild, Sebastian and Herman, Grzegorz},
+  publisher =	{Schloss Dagstuhl -- Leibniz-Zentrum f{\"u}r Informatik},
+  address =	{Dagstuhl, Germany},
+  doi =		{10.4230/LIPIcs.ESA.2025.61},
+  annote =	{Keywords: Hyperbolic surface, Delaunay triangulation, Data structure, Combinatorial map, Implementation, CGAL}
+}
+
 @article{cgal:dfg-cvtaa-99t,
   title={{Centroidal Voronoi Tessellations: Applications and Algorithms}},
   author={Du, Q. and Faber, V. and Gunzburger, M.},

Hyperbolic_triangulation_2/doc/Hyperbolic_triangulation_2/Concepts/HyperbolicDelaunayTriangulationTraits_2.h

@@ -80,7 +80,7 @@ class HyperbolicDelaunayTriangulationTraits_2 {
                 A predicate object. Must provide the function operator
 
                 `Oriented_side operator()(Hyperbolic_point_2 p, Hyperbolic_point_2 q, Hyperbolic_point_2 query),`
-                which returns the position of the point `query` relative to the oriented hyperbolic
+                which returns the position of the point `query` relative to the circle supporting the hyperbolic
                 segment with vertices `p` and `q`.
         */
         typedef unspecified_type                         Side_of_oriented_hyperbolic_segment_2;
@@ -105,6 +105,17 @@ class HyperbolicDelaunayTriangulationTraits_2 {
         */
         typedef unspecified_type                         Is_Delaunay_hyperbolic;
 
+         /*!
+                A predicate object. Must provide the function operator
+
+                `Orientation operator()(Hyperbolic_point_2 p, Hyperbolic_point_2 q, Hyperbolic_point_2 query),`
+
+                which returns the position of the point `query` relative to the oriented hyperbolic
+                segment with vertices `p` and `q`.
+
+        */
+        typedef unspecified_type                         Hyperbolic_orientation_2;
+
   /// @}
 
   /// \name Construction Types
@@ -167,6 +178,7 @@ class HyperbolicDelaunayTriangulationTraits_2 {
         Side_of_oriented_circle_2              side_of_oriented_circle_2_object();
         Side_of_oriented_hyperbolic_segment_2          side_of_oriented_hyperbolic_segment_2_object();
         Is_Delaunay_hyperbolic                                               is_Delaunay_hyperbolic_object();
+        Hyperbolic_orientation_2                hyperbolic_orientation_2();
   /// @}
 
   /// \name

Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_CK_traits_2.h

@@ -50,7 +50,6 @@ namespace internal {
                                           const Hyperbolic_point_2& q,
                                           const Hyperbolic_point_2& r) const
     {
-
       Construct_circle_or_line_supporting_bisector<Traits> cclsb(_gt);
 
       Hyperbolic_point_2 po(CGAL::ORIGIN);
@@ -76,6 +75,7 @@ namespace internal {
           typedef typename CK2_Intersection_traits<Traits, Circle_2, Circle_2>::type Intersection_result;
           std::vector< Intersection_result > inters;
           intersection(*c_pq, *c_qr, std::back_inserter(inters));
+          CGAL_assertion(inters.size() != 0);
 
           CGAL_assertion_code(bool ok=)
           assign(pair, inters[0]);
@@ -107,6 +107,7 @@ namespace internal {
       typedef typename CK2_Intersection_traits<Traits, Euclidean_line_2, Circle_2>::type Intersection_result;
       std::vector< Intersection_result > inters;
       intersection(*l, *c, std::back_inserter(inters));
+      CGAL_assertion(inters.size() != 0);
 
       CGAL_assertion_code(bool ok=)
       assign(pair,inters[0]);
@@ -366,6 +367,7 @@ class Hyperbolic_Delaunay_triangulation_CK_traits_2
   typedef typename internal::Construct_circle_or_line_supporting_bisector<Self> Construct_circle_or_line_supporting_bisector;
   typedef internal::Construct_hyperbolic_segment_2<Self>        Construct_hyperbolic_segment_2;
   typedef typename Base::Construct_segment_2                    Construct_segment_2;
+  typedef typename internal::Hyperbolic_orientation_2<Self>     Hyperbolic_orientation_2;
 
 public:
   Hyperbolic_Delaunay_triangulation_CK_traits_2(const Base& kernel = Base()) : Base(kernel) {}
@@ -411,6 +413,10 @@ class Hyperbolic_Delaunay_triangulation_CK_traits_2
   compute_squared_Euclidean_distance_2_object() const
   { return this->Base::compute_squared_distance_2_object(); }
 
+  Hyperbolic_orientation_2
+  hyperbolic_orientation_2() const
+  { return Hyperbolic_orientation_2(*this); }
+
 };
 
 // Take out the code below to some separate file

Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h

@@ -492,6 +492,7 @@ class Hyperbolic_Delaunay_triangulation_traits_2
   typedef internal::Construct_hyperbolic_bisector_2<Self>       Construct_hyperbolic_bisector_2;
   typedef internal::Is_Delaunay_hyperbolic<Self>                Is_Delaunay_hyperbolic;
   typedef internal::Side_of_oriented_hyperbolic_segment_2<Self> Side_of_oriented_hyperbolic_segment_2;
+  typedef typename internal::Hyperbolic_orientation_2<Self>     Hyperbolic_orientation_2;
 
   // Needed for P4HT2
   typedef typename Kernel::Construct_bisector_2                 Construct_Euclidean_bisector_2;
@@ -549,6 +550,10 @@ class Hyperbolic_Delaunay_triangulation_traits_2
   compute_squared_Euclidean_distance_2_object() const
   { return this->Base::compute_squared_distance_2_object(); }
 
+  Hyperbolic_orientation_2
+  hyperbolic_orientation_2() const
+  { return Hyperbolic_orientation_2(*this); }
+
 }; // class Hyperbolic_Delaunay_triangulation_traits_2
 
 } // namespace CGAL

Hyperbolic_triangulation_2/include/CGAL/Hyperbolic_triangulation_2/internal/Hyperbolic_Delaunay_triangulation_traits_2_functions.h

@@ -258,6 +258,45 @@ class Side_of_oriented_hyperbolic_segment_2
   const Traits& _gt;
 };
 
+template <typename Traits>
+class Hyperbolic_orientation_2
+{
+  typedef typename Traits::FT                                 FT;
+  typedef typename Traits::Hyperbolic_point_2                 Hyperbolic_point_2;
+  typedef typename Traits::Side_of_oriented_hyperbolic_segment_2 Side_of_oriented_hyperbolic_segment_2;
+
+public:
+  Hyperbolic_orientation_2(const Traits& gt = Traits()) : _gt(gt) {}
+
+  Orientation operator()(const Hyperbolic_point_2& p,
+                           const Hyperbolic_point_2& q,
+                           const Hyperbolic_point_2& query) const
+  {
+    Hyperbolic_point_2 origin = Hyperbolic_point_2(FT(0), FT(0));
+    Orientation orientation_to_origin = orientation(p, origin, q);
+    if (orientation_to_origin == COLLINEAR) {
+      return orientation(p, q, query);
+    }
+
+    Side_of_oriented_hyperbolic_segment_2 orientation_test = _gt.side_of_oriented_hyperbolic_segment_2_object();
+    Oriented_side orientation_to_disk = orientation_test(p, q, query);
+    if (orientation_to_disk == ON_POSITIVE_SIDE) {
+      return orientation_to_origin;
+    } else if (orientation_to_disk == ON_NEGATIVE_SIDE) {
+      if (orientation_to_origin == LEFT_TURN) {
+        return RIGHT_TURN;
+      } else {
+        return LEFT_TURN;
+      }
+    } else {
+      return COLLINEAR;
+    }
+  }
+
+private:
+  const Traits& _gt;
+};
+
 } // end namespace internal
 
 } // end namespace CGAL

Triangulation_on_hyperbolic_surface_2/demo/Triangulation_on_hyperbolic_surface_2/CMakeLists.txt

@@ -12,15 +12,32 @@ set(CMAKE_AUTOMOC ON)
 # CGAL and its components
 find_package(CGAL REQUIRED COMPONENTS Core Qt6)
 find_package(Qt6 QUIET COMPONENTS Widgets)
+find_package(Qt6 QUIET COMPONENTS Svg)
 
 if (NOT CGAL_Qt6_FOUND  OR NOT Qt6_FOUND)
   message(STATUS "This project requires the Qt6 library, and will not be compiled.")
   return()
 endif()
 
+# Boost and its components
+find_package( Boost REQUIRED )
+
+if ( NOT Boost_FOUND )
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+  return()
+endif()
+
 # ui files, created with Qt Designer
 qt6_wrap_ui(UIS drawing_window_description.ui)
 
 add_executable( Triangulation_on_hyperbolic_surface_2_demo Triangulation_on_hyperbolic_surface_2_demo.cpp window.cpp ${UIS})
-target_link_libraries(Triangulation_on_hyperbolic_surface_2_demo PRIVATE CGAL::CGAL CGAL::CGAL_Qt6 Qt6::Widgets )
+target_link_libraries(Triangulation_on_hyperbolic_surface_2_demo PRIVATE CGAL::CGAL CGAL::CGAL_Qt6 Qt6::Widgets Qt6::Svg)
 cgal_add_compilation_test(Triangulation_on_hyperbolic_surface_2_demo)
+
+add_executable( Delaunay_triangulation_on_hyperbolic_surface_2_demo Delaunay_triangulation_on_hyperbolic_surface_2_demo.cpp window.cpp ${UIS})
+target_link_libraries(Delaunay_triangulation_on_hyperbolic_surface_2_demo PRIVATE CGAL::CGAL CGAL::CGAL_Qt6 Qt6::Widgets Qt6::Svg)
+
+add_executable( Hyperbolic_Dirichlet_domain_2_demo Hyperbolic_Dirichlet_domain_2_demo.cpp window.cpp ${UIS})
+target_link_libraries(Hyperbolic_Dirichlet_domain_2_demo PRIVATE CGAL::CGAL CGAL::CGAL_Qt6 Qt6::Widgets Qt6::Svg)
+
+set(CMAKE_BUILD_TYPE "Release")

Triangulation_on_hyperbolic_surface_2/demo/Triangulation_on_hyperbolic_surface_2/Delaunay_triangulation_on_hyperbolic_surface_2_demo.cpp

@@ -0,0 +1,110 @@
+#include "window.h"
+
+#include <CGAL/Gmpq.h>
+#include <CGAL/Circular_kernel_2.h>
+#include <CGAL/Algebraic_kernel_for_circles_2_2.h>
+#include <CGAL/Hyperbolic_Delaunay_triangulation_CK_traits_2.h>
+#include <CGAL/Hyperbolic_surface_traits_2.h>
+#include <CGAL/Hyperbolic_fundamental_domain_factory_2.h>
+#include <CGAL/Delaunay_triangulation_on_hyperbolic_surface_2.h>
+
+#include <CGAL/Timer.h>
+
+using namespace CGAL;
+
+typedef Gmpq                                                      NumberType;
+typedef Circular_kernel_2<Simple_cartesian<NumberType>,Algebraic_kernel_for_circles_2_2<NumberType>> Kernel;
+typedef Hyperbolic_Delaunay_triangulation_CK_traits_2<Kernel>     ParentTraits;
+typedef Hyperbolic_surface_traits_2<ParentTraits>                 Traits;
+
+typedef Hyperbolic_fundamental_domain_2<Traits>                   Domain;
+typedef Hyperbolic_isometry_2<Traits>                             Isometry;
+typedef Hyperbolic_fundamental_domain_factory_2<Traits>           Factory;
+typedef Delaunay_triangulation_on_hyperbolic_surface_2<Traits>    Delaunay_triangulation;
+
+/*
+  HOW TO USE THIS DEMO
+  ./Delaunay_triangulation_on_hyperbolic_surface_2_demo [epsilon] [surface seed] [precision]
+  Without arguments, uses the default values defined below.
+*/
+
+// DEFAULT VALUES
+double epsilon = 0.25;
+int seed = time(NULL);
+int p = 1;
+
+int main(int argc, char *argv[])
+{
+  if(argc == 1) {
+    std::cout << "HOW TO USE THIS DEMO:\n"
+              << "./Delaunay_triangulation_on_hyperbolic_surface_2_demo [epsilon] [surface seed] [precision]\n"
+              << "Without arguments, uses default values: epsilon = 0.25, random seed, precision = 1.\n"
+              << "WARNING: when not using the CGAL::Gmpq number type, precision is ignored and, instead, to_double is used to round coordinates of circumcenters.\n"
+              << "--------------------"
+              << std::endl;
+  }
+  // 1. Parse args and generate the input
+  if (argc > 1) {
+    epsilon = std::stod(argv[1]);
+  }
+
+  Domain domain;
+  if (argc <= 2) {
+    std::cout << "Using random seed " << seed << std::endl;
+  } else {
+    seed = atoi(argv[2]);
+  }
+  Factory factory;
+  std::cout << "Generating surface with seed " << seed << "..." << std::endl;
+  domain = factory.make_hyperbolic_fundamental_domain_g2(seed);
+  Delaunay_triangulation dt = Delaunay_triangulation(domain);
+
+  if (argc > 3) {
+    p = atoi(argv[3]);
+  }
+
+  // 2. Get a vertex
+  // So that if you run the demo on a same surface but with different values of epsilon,
+  // the drawing will be centered at the same vertex and it will look similar.
+  Point v0 = dt.anchor().vertices[0];
+
+  // 3. Compute epsilon-net and display useful info
+  if constexpr(!std::is_same<NumberType, Gmpq>::value) {
+    std::cout << "WARNING: Not using the CGAL::Gmpq number type. Precision will be ignored and to_double approximation will be used instead." << std::endl;
+  }
+  std::cout << "Computing a " << epsilon << "-net with floating-point precision " << p*53 << "..." << std::endl;
+  Timer timer;
+  timer.start();
+  std::cout << "Is epsilon-net? " << dt.epsilon_net(epsilon, p) << std::endl;
+  timer.stop();
+  std::cout << "Done in " << timer.time() << " seconds." << std::endl;
+  dt.combinatorial_map().display_characteristics(std::cout) << std::endl;
+
+  // 4. SET THE FIRST ANCHOR OF THE DRAWING
+  Anchor anchor = dt.locate(v0);
+  int index = 0;
+  for (int i = 0; i < 3; i++) {
+    if (v0 == anchor.vertices[i]) {
+      index = i;
+    }
+  }
+
+  Anchor start = Anchor();
+  start.dart = anchor.dart;
+  for (int i = 0; i < 3; i++) {
+    start.vertices[i] = anchor.vertices[(i + index) % 3];
+    if (i < index) {
+      start.dart = dt.Base::ccw(start.dart);
+    }
+  }
+
+  // 5. DRAW
+  QApplication app(argc, argv);
+  app.setApplicationName("Delaunay triangulation on hyperbolic surface 2 Demo");
+  DemoWindow window;
+  window.item().draw_triangulation(dt, start);
+  window.show();
+  QStringList args = app.arguments();
+  args.removeAt(0);
+  return app.exec();
+}

Triangulation_on_hyperbolic_surface_2/demo/Triangulation_on_hyperbolic_surface_2/Hyperbolic_Dirichlet_domain_2_demo.cpp

@@ -0,0 +1,35 @@
+#include "window.h"
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Circular_kernel_2.h>
+#include <CGAL/Algebraic_kernel_for_circles_2_2.h>
+#include <CGAL/Hyperbolic_Delaunay_triangulation_CK_traits_2.h>
+#include <CGAL/Hyperbolic_surface_traits_2.h>
+
+#include <CGAL/Hyperbolic_Dirichlet_domain_2.h>
+#include <CGAL/Hyperbolic_fundamental_domain_factory_2.h>
+
+typedef CGAL::Gmpq                                                                                              NumberType;
+typedef CGAL::Circular_kernel_2<CGAL::Simple_cartesian<NumberType>,CGAL::Algebraic_kernel_for_circles_2_2<NumberType>> Kernel;
+typedef CGAL::Hyperbolic_Delaunay_triangulation_CK_traits_2<Kernel>                                             ParentTraits;
+typedef CGAL::Hyperbolic_surface_traits_2<ParentTraits>                                                         Traits;
+
+typedef CGAL::Hyperbolic_fundamental_domain_2<Traits>                                                           Domain;
+typedef CGAL::Hyperbolic_fundamental_domain_factory_2<Traits>                                                   Factory;
+
+int main(int argc, char **argv)
+{
+    int seed = time(NULL);
+    std::cout << "Using seed " << seed << std::endl;
+    Factory factory;
+    Domain domain = factory.make_hyperbolic_fundamental_domain_g2(seed);
+
+    QApplication app(argc, argv);
+    app.setApplicationName("Hyperbolic Dirichlet domain 2 Demo");
+    DemoWindow window;
+    window.item().draw_Dirichlet(domain);
+    window.show();
+    QStringList args = app.arguments();
+    args.removeAt(0);
+    return app.exec();
+}

Triangulation_on_hyperbolic_surface_2/demo/Triangulation_on_hyperbolic_surface_2/Triangulation_on_hyperbolic_surface_2_demo.cpp

@@ -1,17 +1,22 @@
 #include "window.h"
 
 #include <CGAL/Exact_rational.h>
+#include <CGAL/Gmpq.h>
 #include <CGAL/Simple_cartesian.h>
 
-#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
+#include <CGAL/Hyperbolic_Delaunay_triangulation_CK_traits_2.h>
 #include <CGAL/Hyperbolic_surface_traits_2.h>
 #include <CGAL/Hyperbolic_fundamental_domain_factory_2.h>
 #include <CGAL/Triangulation_on_hyperbolic_surface_2.h>
 
 using namespace CGAL;
 
-typedef Simple_cartesian<Exact_rational>                                Kernel;
-typedef Hyperbolic_Delaunay_triangulation_traits_2<Kernel>              ParentTraits;
+typedef Gmpq                                                      NumberType;
+typedef Circular_kernel_2<Simple_cartesian<NumberType>, Algebraic_kernel_for_circles_2_2<NumberType>> Kernel;
+typedef Hyperbolic_Delaunay_triangulation_CK_traits_2<Kernel>     ParentTraits;
+
+// typedef Simple_cartesian<NumberType>                                Kernel;
+// typedef Hyperbolic_Delaunay_triangulation_traits_2<Kernel>              ParentTraits;
 typedef Hyperbolic_surface_traits_2<ParentTraits>                       Traits;
 typedef Hyperbolic_fundamental_domain_2<Traits>                         Domain;
 typedef Hyperbolic_fundamental_domain_factory_2<Traits>                 Factory;
@@ -31,7 +36,7 @@ int main(int argc, char** argv)
 
   // 2. Draw the triangulation
   QApplication app(argc, argv);
-  app.setApplicationName("Hyperbolic surfaces triangulation 2 Demo");
+  app.setApplicationName("Triangulation on hyperbolic surface 2 Demo");
 
   DemoWindow window;
   window.item().draw_triangulation(triangulation);