Merge pull request #5827 from MaelRL/ToS2-VD2_AP-GF

Add VD2's removal policy for ToS2 (+ caching)

Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h

@@ -147,6 +147,16 @@ public:
   {
   }
 
+  // Assignement
+  Delaunay_triangulation_on_sphere_2& operator=(Delaunay_triangulation_on_sphere_2 other) // intentional copy
+  {
+    Base::swap(static_cast<Base&>(other));
+    return *this;
+  }
+
+  // Destructor
+  ~Delaunay_triangulation_on_sphere_2() = default;
+
   // Predicates & Constructions
   Oriented_side side_of_oriented_circle(const Point& p, const Point& q, const Point& r, const Point& s, bool perturb = false) const;
   Oriented_side side_of_oriented_circle(const Face_handle f, const Point& p, bool perturb = false) const;

Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h

@@ -127,6 +127,9 @@ public:
   void swap(Triangulation_on_sphere_2& tr);
   Triangulation_on_sphere_2& operator=(Triangulation_on_sphere_2 tr); // intentional copy
 
+  // Destructor
+  ~Triangulation_on_sphere_2() = default;
+
 public:
   // Members
   const Geom_traits& geom_traits() const { return _gt; }

Voronoi_diagram_2/doc/Voronoi_diagram_2/Concepts/AdaptationPolicy_2.h

@@ -65,14 +65,12 @@ typedef Delaunay_graph::Edge Delaunay_edge;
 /*!
 
 */
-typedef Delaunay_graph::All_edges_iterator
+typedef Delaunay_graph::All_edges_iterator All_Delaunay_edges_iterator;
-All_Delaunay_edges_iterator;
 
 /*!
 
 */
-typedef Delaunay_graph::Finite_edges_iterator
+typedef Delaunay_graph::Finite_edges_iterator Finite_Delaunay_edges_iterator;
-Finite_Delaunay_edges_iterator;
 
 /*!
 
@@ -93,13 +91,9 @@ model of the concepts `DefaultConstructible`,
 
 `bool operator()(Delaunay_graph dg, Delaunay_edge_circulator ec)`
 
-`bool operator()(Delaunay_graph dg,`
+`bool operator()(Delaunay_graph dg, All_Delaunay_edges_iterator eit)`
 
-`All_Delaunay_edges_iterator eit)`
+`bool operator()(Delaunay_graph dg, Finite_Delaunay_edges_iterator eit)`
-
-`bool operator()(Delaunay_graph dg,`
-
-`Finite_Delaunay_edges_iterator eit)`
 
 The functor returns `true` iff the edge is rejected.
 */

Voronoi_diagram_2/include/CGAL/Delaunay_triangulation_on_sphere_adaptation_policies_2.h

@@ -0,0 +1,65 @@
+// Copyright (c) 2021 GeometryFactory
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Menelaos Karavelas <mkaravel@iacm.forth.gr>
+//                 Mael Rouxel-Labbé
+
+#ifndef CGAL_DELAUNAY_TRIANGULATION_ON_SPHERE_ADAPTATION_POLICIES_2_H
+#define CGAL_DELAUNAY_TRIANGULATION_ON_SPHERE_ADAPTATION_POLICIES_2_H 1
+
+#include <CGAL/license/Voronoi_diagram_2.h>
+
+#include <CGAL/Voronoi_diagram_2/basic.h>
+#include <CGAL/Voronoi_diagram_2/Delaunay_triangulation_on_sphere_degeneracy_testers.h>
+#include <CGAL/Voronoi_diagram_2/Policy_base.h>
+#include <CGAL/Voronoi_diagram_2/Default_site_inserters.h>
+#include <CGAL/Voronoi_diagram_2/Identity_rejectors.h>
+
+#include <CGAL/Identity_policy_2.h>
+
+namespace CGAL {
+
+//=========================================================================
+//=========================================================================
+
+template<class DToS2>
+struct Delaunay_triangulation_on_sphere_degeneracy_removal_policy_2
+  : public CGAL_VORONOI_DIAGRAM_2_INS::Policy_base
+  <DToS2,
+   CGAL_VORONOI_DIAGRAM_2_INS::Delaunay_triangulation_on_sphere_edge_tester_2<DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Identity_face_rejector<DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Default_site_inserter<typename DToS2::Point,DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Default_site_remover<DToS2> >
+{
+  typedef typename DToS2::Point      Site_2;
+};
+
+
+//=========================================================================
+//=========================================================================
+
+template<class DToS2>
+struct Delaunay_triangulation_on_sphere_caching_degeneracy_removal_policy_2
+  : public CGAL_VORONOI_DIAGRAM_2_INS::Caching_policy_base
+  <DToS2,
+   CGAL_VORONOI_DIAGRAM_2_INS::Delaunay_triangulation_on_sphere_edge_tester_2<DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Identity_face_rejector<DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Default_site_inserter<typename DToS2::Point,DToS2>,
+   CGAL_VORONOI_DIAGRAM_2_INS::Default_site_remover<DToS2> >
+{
+  typedef typename DToS2::Point      Site_2;
+};
+
+//=========================================================================
+//=========================================================================
+
+} // namespace CGAL
+
+#endif // CGAL_DELAUNAY_TRIANGULATION_ON_SPHERE_ADAPTATION_POLICIES_2_H

Voronoi_diagram_2/include/CGAL/Voronoi_diagram_2/Cached_degeneracy_testers.h

@@ -99,18 +99,9 @@ public:
     return operator()(dual, Edge(f,i));
   }
 
+  template <typename EdgeIterator>
   bool operator()(const Delaunay_graph& dual,
-                  const Edge_circulator& ec) const {
+                  const EdgeIterator eit) const {
-    return operator()(dual, *ec);
-  }
-
-  bool operator()(const Delaunay_graph& dual,
-                  const All_edges_iterator& eit) const {
-    return operator()(dual, *eit);
-  }
-
-  bool operator()(const Delaunay_graph& dual,
-                  const Finite_edges_iterator& eit) const {
     return operator()(dual, *eit);
   }
 
@@ -214,18 +205,9 @@ public:
     return operator()(dual, Edge(f,i));
   }
 
+  template <typename EdgeIterator>
   bool operator()(const Delaunay_graph& dual,
-                  const Edge_circulator& ec) const {
+                  const EdgeIterator eit) const {
-    return operator()(dual, *ec);
-  }
-
-  bool operator()(const Delaunay_graph& dual,
-                  const All_edges_iterator& eit) const {
-    return operator()(dual, *eit);
-  }
-
-  bool operator()(const Delaunay_graph& dual,
-                  const Finite_edges_iterator& eit) const {
     return operator()(dual, *eit);
   }
 
@@ -250,11 +232,11 @@ public:
 
   bool is_valid() const { return true; }
 
-  bool is_valid(const Delaunay_graph& dual) const {
+  bool is_valid(const Delaunay_graph& dual) const
+  {
     bool valid = true;
-    All_edges_iterator eit;
+    for (auto eit = dual.all_edges_begin(); eit != dual.all_edges_end(); ++eit) {
-    for (eit = dual.all_edges_begin(); eit != dual.all_edges_end(); ++eit) {
+      const Edge& e = *eit;
-      Edge e = *eit;
       bool b = !emap.is_defined(e) || (emap[e] != UNDEFINED);
       valid = valid && b;
     }

Voronoi_diagram_2/include/CGAL/Voronoi_diagram_2/Delaunay_triangulation_on_sphere_degeneracy_testers.h

@@ -0,0 +1,92 @@
+// Copyright (c) 2021 GeometryFactory.
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Menelaos Karavelas <mkaravel@iacm.forth.gr>
+//                 Mael Rouxel-Labbé
+
+#ifndef CGAL_VORONOI_DIAGRAM_2_DELAUNAY_TRIANGULATION_ON_SPHERE_DEGENERACY_TESTERS_H
+#define CGAL_VORONOI_DIAGRAM_2_DELAUNAY_TRIANGULATION_ON_SPHERE_DEGENERACY_TESTERS_H 1
+
+#include <CGAL/license/Voronoi_diagram_2.h>
+
+#include <CGAL/Voronoi_diagram_2/basic.h>
+#include <CGAL/Voronoi_diagram_2/Adaptation_traits_base_2.h>
+#include <CGAL/Voronoi_diagram_2/Identity_rejectors.h>
+
+namespace CGAL {
+namespace VoronoiDiagram_2 {
+namespace Internal {
+
+//=========================================================================
+//=========================================================================
+
+// tests whether a dual edge has zero length
+template <class DG>
+class Delaunay_triangulation_on_sphere_edge_tester_2
+  : public CGAL_VORONOI_DIAGRAM_2_INS::Rejector_base
+{
+public:
+  typedef DG                                             Delaunay_graph;
+
+  typedef typename Delaunay_graph::Edge                  Edge;
+  typedef typename Delaunay_graph::Face_handle           Face_handle;
+  typedef typename Delaunay_graph::Edge_circulator       Edge_circulator;
+  typedef typename Delaunay_graph::All_edges_iterator    All_edges_iterator;
+  typedef typename Delaunay_graph::Finite_edges_iterator Finite_edges_iterator;
+  typedef bool                                           result_type;
+
+ private:
+  typedef typename Delaunay_graph::Geom_traits           Geom_traits;
+  typedef typename Delaunay_graph::Vertex_handle         Vertex_handle;
+  typedef typename Delaunay_graph::Point                 Point;
+
+ public:
+  bool operator()(const Delaunay_graph& dual,
+                  const Face_handle f, int i) const
+  {
+    if(dual.dimension() == 1)
+      return false;
+
+    const Vertex_handle v3 = f->vertex(i);
+    const Vertex_handle v4 = dual.tds().mirror_vertex(f, i);
+
+    const Vertex_handle v1 = f->vertex(dual.ccw(i));
+    const Vertex_handle v2 = f->vertex(dual.cw(i));
+
+    const Point& p1 = v1->point();
+    const Point& p2 = v2->point();
+    const Point& p3 = v3->point();
+    const Point& p4 = v4->point();
+    const Oriented_side os = dual.geom_traits().side_of_oriented_circle_on_sphere_2_object()(p1,p2,p3,p4);
+
+    return (os == ON_ORIENTED_BOUNDARY);
+  }
+
+  bool operator()(const Delaunay_graph& dual, const Edge& e) const
+  {
+    return operator()(dual, e.first, e.second);
+  }
+
+  template <typename EdgeIterator>
+  bool operator()(const Delaunay_graph& dual,
+                  const EdgeIterator eit) const
+  {
+    return operator()(dual, *eit);
+  }
+};
+
+//=========================================================================
+//=========================================================================
+
+} // namespace Internal
+} // namespace VoronoiDiagram_2
+} // namespace CGAL
+
+#endif // CGAL_VORONOI_DIAGRAM_2_DELAUNAY_TRIANGULATION_ON_SPHERE_DEGENERACY_TESTERS_H

Voronoi_diagram_2/test/Voronoi_diagram_2/CMakeLists.txt

@@ -8,11 +8,16 @@ find_package(CGAL REQUIRED)
 
 include_directories(BEFORE "include")
 
-# create a target per cppfile
+create_single_source_cgal_program("vda_ag.cpp")
-file(
+create_single_source_cgal_program("vda_dt.cpp")
-  GLOB cppfiles
+create_single_source_cgal_program("vda_pt.cpp")
-  RELATIVE ${CMAKE_CURRENT_SOURCE_DIR}
+create_single_source_cgal_program("vda_rt.cpp")
-  ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+create_single_source_cgal_program("vda_sdg.cpp")
-foreach(cppfile ${cppfiles})
+
-  create_single_source_cgal_program("${cppfile}")
+find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater)
-endforeach()
+include(CGAL_Eigen3_support)
+
+if(TARGET CGAL::Eigen3_support)
+  create_single_source_cgal_program("vda_tos2.cpp")
+  target_link_libraries(vda_tos2 PUBLIC CGAL::Eigen3_support)
+endif()

Voronoi_diagram_2/test/Voronoi_diagram_2/vda_tos2.cpp

@@ -0,0 +1,201 @@
+#include <CGAL/Voronoi_diagram_2.h>
+#include "vda_test.h"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_on_sphere_adaptation_traits_2.h>
+#include <CGAL/Projection_on_sphere_traits_3.h>
+#include <CGAL/Delaunay_triangulation_on_sphere_adaptation_policies_2.h>
+#include <CGAL/Delaunay_triangulation_on_sphere_2.h>
+
+#include <iostream>
+#include <iterator>
+#include <fstream>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::FT FT;
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+
+typedef CGAL::Projection_on_sphere_traits_3<K> Traits;
+typedef CGAL::Triangulation_on_sphere_face_base_2<Traits> Fb;
+typedef CGAL::Triangulation_on_sphere_vertex_base_2<Traits> Vb;
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
+typedef CGAL::Delaunay_triangulation_on_sphere_2<Traits, Tds> DToS2;
+
+typedef Traits::SK SK;
+typedef Traits::Point_on_sphere_2 Point_on_sphere_2;
+typedef Traits::Arc_on_sphere_2 Arc_on_sphere_2;
+typedef DToS2::Face_handle Face_handle;
+
+typedef CGAL::Delaunay_triangulation_on_sphere_adaptation_traits_2<DToS2> AT;
+typedef CGAL::Delaunay_triangulation_on_sphere_caching_degeneracy_removal_policy_2<DToS2> AP;
+typedef CGAL::Voronoi_diagram_2<DToS2, AT, AP> VD;
+
+int main(int argc, char** argv)
+{
+  const char* filename = (argc > 1) ? argv[1] : "data/radar.xyz";
+  const double radius = (argc > 2) ? std::stod(argv[2]) : 1;
+  const double step = (argc > 3) ? std::stod(argv[3]) : 0.01;
+
+  std::vector<Point_3> points;
+  double x, y, z;
+
+  std::ifstream in(filename);
+  if(!in)
+  {
+    std::cerr << "Invalid input file: " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  while(in >> x >> y >> z)
+    points.emplace_back(x, y, z);
+
+  std::cout << points.size() << " points" << std::endl;
+
+  Traits traits(Point_3(1,1,1), radius); // radius is 1 by default
+  DToS2 dtos(traits);
+
+  Traits::Construct_point_on_sphere_2 cst = traits.construct_point_on_sphere_2_object();
+
+  for(const auto& pt : points)
+  {
+//    std::cout << "----- Inserting (" << pt
+//              << ") at squared distance " << CGAL::squared_distance(pt, traits.center())
+//              << " from the center of the sphere" << std::endl;
+
+    // shenanigans just to get a OFF with pts on the sphere (it should just be 'cst(pt)')
+    dtos.insert(cst(pt).get_projection(traits.center(), traits.radius()));
+
+//    std::cout << "The triangulation now has dimension: " << dtos.dimension() << " and\n";
+//    std::cout << dtos.number_of_vertices() << " vertices" << std::endl;
+//    std::cout << dtos.number_of_edges() << " edges" << std::endl;
+//    std::cout << dtos.number_of_faces() << " solid faces" << std::endl;
+//    std::cout << dtos.number_of_ghost_faces() << " ghost faces" << std::endl;
+  }
+
+#ifdef CGAL_VD2_TEST_TOS2_OUTPUT
+  CGAL::write_OFF("result.off", dtos, CGAL::parameters::stream_precision(17));
+
+  std::ofstream out_primal("edges_primal.polylines.cgal");
+  out_primal.precision(17);
+  std::ofstream out_dual("edges_dual.polylines.cgal");
+  out_dual.precision(17);
+#endif
+
+  for(typename DToS2::All_faces_iterator fit = dtos.all_faces_begin(); fit!=dtos.all_faces_end(); ++fit)
+  {
+    Face_handle fh = fit;
+
+    const Point_3 c = dtos.dual(fh);
+    CGAL_USE(c);
+
+    const Point_on_sphere_2 cs = dtos.dual_on_sphere(fh);
+    const FT r = dtos.geom_traits().radius();
+    assert(CGAL::abs(CGAL::squared_distance(dtos.construct_point(cs),
+                                            dtos.geom_traits().center()) - r*r) <= 1e-10);
+
+    for(int i=0; i<3; ++i)
+    {
+      Face_handle nfh = fh->neighbor(i);
+      if(/*!dtos.is_ghost(nfh) &&*/ fh > nfh)
+        continue;
+
+      typename DToS2::Edge e(fh, i);
+      const bool diametral_edge = CGAL::collinear(dtos.construct_point(dtos.point(fh, (i+1)%3)),
+                                                  dtos.construct_point(dtos.point(fh, (i+2)%3)),
+                                                  dtos.geom_traits().center());
+
+      // primal
+      if(!diametral_edge)
+      {
+        const Arc_on_sphere_2 as = dtos.segment_on_sphere(e);
+        std::vector<typename SK::Point_3> discretization_points;
+        CGAL::Triangulations_on_sphere_2::internal::subsample_arc_on_sphere_2<SK>(as, std::back_inserter(discretization_points), step);
+        assert(discretization_points.size() >= 2);
+
+#ifdef CGAL_VD2_TEST_TOS2_OUTPUT
+        for(std::size_t i=0; i<discretization_points.size()-1; ++i)
+          out_primal << "2 " << discretization_points[i] << " " << discretization_points[i+1] << "\n";
+#endif
+      }
+      else
+      {
+#ifdef CGAL_VD2_TEST_TOS2_OUTPUT
+        const Segment_3 s = dtos.segment(e);
+        out_primal << "2 " << s.source() << " " << s.target() << "\n";
+#endif
+      }
+
+      // Dual
+//      if(dtos.is_contour(e))
+//        continue;
+
+      const Point_on_sphere_2 c1 = dtos.circumcenter_on_sphere(e.first);
+      const Point_on_sphere_2 c2 = dtos.circumcenter_on_sphere(dtos.mirror_edge(e).first);
+
+      // That should never be possible, but with constructions...
+      const bool diametral_dual = CGAL::collinear(dtos.construct_point(c1),
+                                                  dtos.construct_point(c2),
+                                                  dtos.geom_traits().center());
+      if(!diametral_dual)
+      {
+        const Arc_on_sphere_2 ad = dtos.dual_on_sphere(e);
+        std::vector<typename SK::Point_3> discretization_points;
+        CGAL::Triangulations_on_sphere_2::internal::subsample_arc_on_sphere_2<SK>(ad, std::back_inserter(discretization_points), step);
+        assert(discretization_points.size() >= 2);
+
+#ifdef CGAL_VD2_TEST_TOS2_OUTPUT
+        for(std::size_t i=0; i<discretization_points.size()-1; ++i)
+          out_dual << "2 " << discretization_points[i] << " " << discretization_points[i+1] << "\n";
+#endif
+      }
+      else
+      {
+#ifdef CGAL_VD2_TEST_TOS2_OUTPUT
+        const Segment_3 d = dtos.dual(e);
+        out_dual << "2 " << d.source() << " " << d.target() << "\n";
+#endif
+      };
+    }
+  }
+
+  std::cout << "Voronoi diagram:" << std::endl;
+  AT adaptation_traits(dtos.geom_traits());
+
+  VD vd(dtos, false /*don't give the diagram ownership of the dt*/, adaptation_traits);
+
+  std::cout << vd.number_of_vertices() << " vertices" << std::endl;
+  std::cout << vd.number_of_halfedges() << " halfedges" << std::endl;
+  std::cout << vd.number_of_faces() << " faces" << std::endl;
+  std::cout << "dimension = " << vd.dual().dimension() << std::endl;
+
+  VD::Vertex_iterator vit = vd.vertices_begin(), vend = vd.vertices_end();
+  for(; vit!=vend; ++vit)
+  {
+    std::cout << "Voronoi vertex at " << vit->point() << std::endl;
+    std::cout << "valence: " << vit->degree() << std::endl;
+    DToS2::Face_handle dual_face = vit->dual();
+    CGAL_USE(dual_face);
+
+    VD::Halfedge_around_vertex_circulator cir = vd.incident_halfedges(vit), cend = cir;
+    do {
+      std::cout << "  Incident Halfedge with source: " << cir->source()->point() << std::endl;
+      VD::Face_handle incident_Voronoi_face = cir->face();
+      DToS2::Vertex_handle Voronoi_face_seed = incident_Voronoi_face->dual();
+      CGAL_USE(Voronoi_face_seed);
+    } while(cir != cend);
+  }
+
+  VD::Halfedge_iterator hit = vd.halfedges_begin(), hend = vd.halfedges_end();
+  for(; hit!=hend; ++hit)
+  {
+    std::cout << "Halfedge between " << hit->source()->point() << " and " << hit->target()->point() << std::endl;
+    VD::Halfedge_handle next_h = hit->next();
+    VD::Halfedge_handle opposite_h = hit->opposite();
+    CGAL_USE(next_h);
+    CGAL_USE(opposite_h);
+  }
+
+  return EXIT_SUCCESS;
+}