Merge pull request #8201 from afabri/ToS-IO-GF

Triangulation on Sphere: Fix I/O

# Conflicts:
#	Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h
#	Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h
#	Triangulation_on_sphere_2/test/Triangulation_on_sphere_2/CMakeLists.txt

Triangulation_on_sphere_2/doc/Triangulation_on_sphere_2/Concepts/TriangulationOnSphereFaceBase_2.h

@@ -32,4 +32,16 @@ public:
 
   /// provides write access to a Boolean used to indicate if the face is a ghost face.
   void set_ghost(const bool b);
+
+  /*!
+  inputs the non-combinatorial information given by the face:
+  its ghost status and other possible information.
+  */
+  std::istream& operator>>(std::istream& is, TriangulationOnSphereFaceBase_2& v);
+
+  /*!
+  outputs the non combinatorial operation given by the face:
+  its ghost status and other possible information.
+  */
+  std::ostream& operator<<(std::ostream& os, const TriangulationOnSphereFaceBase_2& v);
 };

Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h

@@ -327,9 +327,8 @@ public:
   Arc_on_sphere_2 dual_on_sphere(const All_edges_iterator ei) const { return dual_on_sphere(*ei); }
 
   // Validity
-  bool is_plane() const;
-  bool is_valid(bool verbose = false, int level = 0) const;
   bool is_valid_face(Face_handle fh, bool verbose = false, int level = 0) const;
+  bool is_valid(bool verbose = false, int level = 0) const;
 };
 
 // ------------------------ PREDICATES / CONSTRUCTIONS --------------------------------//
@@ -1090,39 +1089,6 @@ dual_on_sphere(const Edge& e) const
 
 //-------------------------------------------CHECK------------------------------------------------//
 
-// checks whether a given triangulation is plane (all points are coplanar)
-template <typename Gt, typename Tds>
-bool
-Delaunay_triangulation_on_sphere_2<Gt, Tds>::
-is_plane() const
-{
-  if(number_of_vertices() <= 3)
-    return true;
-
-  bool plane = true;
-
-  Vertices_iterator it1 = vertices_begin(), it2(it1), it3(it1), it4(it1);
-  std::advance(it2, 1);
-  std::advance(it3, 2);
-  std::advance(it4, 3);
-
-  while(it4 != vertices_end())
-  {
-    Orientation s = side_of_oriented_circle(point(it1), point(it2), point(it3), point(it4));
-    plane = plane && s == ON_ORIENTED_BOUNDARY;
-
-    if(!plane)
-      return false;
-
-    ++it1;
-    ++it2;
-    ++it3;
-    ++it4;
-  }
-
-  return true;
-}
-
 template <typename Gt, typename Tds>
 bool
 Delaunay_triangulation_on_sphere_2<Gt, Tds>::
@@ -1154,47 +1120,19 @@ bool
 Delaunay_triangulation_on_sphere_2<Gt, Tds>::
 is_valid(bool verbose, int level) const
 {
-  bool result = true;
-
-  if(!tds().is_valid(verbose, level))
-  {
+  // in any dimension
   if(verbose)
-      std::cerr << "invalid data structure" << std::endl;
-
-    CGAL_assertion(false);
-    return false;
+  {
+    std::cerr << " number of vertices " << number_of_vertices() << "\t" << std::endl;
+    std::cerr << " number of faces " << number_of_faces() << "\t" << std::endl;
   }
 
+  bool result = Base::is_valid(verbose, level);
+  CGAL_assertion(result);
+
   for(All_faces_iterator fit=all_faces_begin(); fit!=all_faces_end(); ++fit)
     result = result && is_valid_face(fit, verbose, level);
 
-  for(Vertices_iterator vit=vertices_begin(); vit!=vertices_end(); ++vit)
-    result = result && Base::is_valid_vertex(vit, verbose, level);
-
-  switch(dimension())
-  {
-    case 0:
-      break;
-    case 1:
-      CGAL_assertion(this->is_plane());
-      break;
-    case 2:
-      for(All_faces_iterator it=all_faces_begin(); it!=all_faces_end(); ++it)
-      {
-        Orientation s = orientation_on_sphere(point(it, 0), point(it, 1), point(it, 2));
-        result = result && (s != NEGATIVE || it->is_ghost());
-        CGAL_assertion(result);
-      }
-
-      result = result && (number_of_faces() == 2 * number_of_vertices() - 4);
-      CGAL_assertion(result);
-      break;
-  }
-
-  // in any dimension
-  if(verbose)
-    std::cerr << " number of vertices " << number_of_vertices() << "\t" << std::endl;
-
   CGAL_assertion(result);
   return result;
 }

Triangulation_on_sphere_2/include/CGAL/Projection_on_sphere_traits_3.h

@@ -37,7 +37,7 @@ private:
   mutable Base_point proj_pt;
 
 public:
-  Point_with_scale() : Base_point() { } // vertex base wants a default constructor
+  Point_with_scale() : Base_point(), cached(false) { } // vertex base wants a default constructor
   Point_with_scale(const Base_point& p) : Base_point(p), cached(false) { }
 
   const Base_point& get_projection(const Base_point& center,

Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h

@@ -482,7 +482,7 @@ public:
 
   //-----------------------DEBUG--------------------------------------------------------------------
   void check_neighboring_faces() const;
-
+  bool is_plane() const;
   bool is_valid_vertex(Vertex_handle fh, bool verbose = false, int level = 0) const;
   bool is_valid(bool verbose = false, int level = 0) const;
 
@@ -1239,6 +1239,39 @@ check_neighboring_faces() const
   }
 }
 
+// checks whether a given triangulation is plane (all points are coplanar)
+template <typename Gt, typename Tds>
+bool
+Triangulation_on_sphere_2<Gt, Tds>::
+is_plane() const
+{
+  if(number_of_vertices() <= 3)
+    return true;
+
+  bool is_plane = true;
+
+  Vertices_iterator it1 = vertices_begin(), it2(it1), it3(it1), it4(it1);
+  std::advance(it2, 1);
+  std::advance(it3, 2);
+  std::advance(it4, 3);
+
+  while(it4 != vertices_end())
+  {
+    Orientation s = orientation(point(it1), point(it2), point(it3), point(it4));
+    is_plane = is_plane && s == COPLANAR;
+
+    if(!is_plane)
+      return false;
+
+    ++it1;
+    ++it2;
+    ++it3;
+    ++it4;
+  }
+
+  return true;
+}
+
 template <typename Gt, typename Tds>
 bool
 Triangulation_on_sphere_2<Gt, Tds>::
@@ -1274,27 +1307,21 @@ is_valid(bool verbose,
   if(dimension() <= 0 || (dimension() == 1 && number_of_vertices() == 2))
     return result;
 
+  for(Vertices_iterator vit=vertices_begin(); vit!=vertices_end(); ++vit)
+    result = result && is_valid_vertex(vit, verbose, level);
+
   if(dimension() == 1)
   {
-    Vertices_iterator it1 = vertices_begin(), it2(it1), it3(it1);
-    std::advance(it2, 1);
-    std::advance(it3, 2);
-    while(it3 != vertices_end())
-    {
-      const Orientation s = orientation(point(it1), point(it2), point(it3));
-      result = result && (s == COLLINEAR);
+    result = result && this->is_plane();
     CGAL_assertion(result);
-      ++it1; ++it2; ++it3;
-    }
   }
   else // dimension() == 2
   {
     for(All_faces_iterator it=all_faces_begin(); it!=all_faces_end(); ++it)
     {
       const Orientation s = orientation_on_sphere(point(it, 0), point(it, 1), point(it, 2));
-      CGAL_assertion(s == LEFT_TURN || is_ghost(it));
-
       result = result && (s == LEFT_TURN || is_ghost(it));
+      CGAL_assertion(result);
     }
 
     // check number of faces. This cannot be done by the TDS,
@@ -1314,6 +1341,7 @@ void
 Triangulation_on_sphere_2<Gt, Tds>::
 file_output(std::ostream& os) const
 {
+  os << _gt.center() << " " << _gt.radius() << "\n";
   _tds.file_output(os, Vertex_handle(), true);
 }
 
@@ -1323,7 +1351,15 @@ Triangulation_on_sphere_2<Gt, Tds>::
 file_input(std::istream& is)
 {
   clear();
-  Vertex_handle v = _tds.file_input(is, true);
+
+  Point_3 center;
+  FT radius;
+  is >> center >> radius;
+  _gt.set_center(center);
+  _gt.set_radius(radius);
+
+  Vertex_handle v = _tds.file_input(is, false);
+  CGAL_assertion(is_valid());
   return v;
 }
 
@@ -1340,7 +1376,6 @@ std::istream&
 operator>>(std::istream& is, Triangulation_on_sphere_2<Gt, Tds>& tr)
 {
   tr.file_input(is);
-  CGAL_assertion(tr.is_valid());
   return is;
 }
 

Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_face_base_2.h

@@ -64,6 +64,27 @@ protected:
   bool _ghost_flag;
 };
 
+template < class Gt, class Fb >
+std::ostream&
+operator<<(std::ostream &os, const Triangulation_on_sphere_face_base_2<Gt, Fb> &f)
+{
+  // non combinatorial information. Default = point
+    os << static_cast<const Fb&>(f);
+    os << (f.is_ghost() ? " 1" : " 0");
+    return os;
+}
+
+template < class Gt, class Fb >
+std::istream&
+operator>>(std::istream &is, Triangulation_on_sphere_face_base_2<Gt, Fb> &f)
+{
+  int g = -1;
+  is >>  static_cast<Fb&>(f);
+  is >> g;
+  f.set_ghost(g == 1);
+  return is;
+}
+
 } // namespace CGAL
 
 #endif //CGAL_TRIANGULATION_ON_SPHERE_FACE_BASE_2_H

Triangulation_on_sphere_2/test/Triangulation_on_sphere_2/CMakeLists.txt

@@ -10,6 +10,7 @@ create_single_source_cgal_program( "test_dtos_degenerate_cases.cpp" )
 create_single_source_cgal_program( "test_dtos_illegal_points.cpp" )
 create_single_source_cgal_program( "test_dtos_projection_traits.cpp" )
 create_single_source_cgal_program( "test_dtos_traits.cpp" )
+create_single_source_cgal_program( "issue_8200.cpp" )
 
 find_package(Eigen3 3.1.0 QUIET) #(requires 3.1.0 or greater)
 include(CGAL_Eigen3_support)

Triangulation_on_sphere_2/test/Triangulation_on_sphere_2/issue_8200.cpp

@@ -0,0 +1,73 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Delaunay_triangulation_on_sphere_2.h>
+#include <CGAL/Projection_on_sphere_traits_3.h>
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel          K;
+
+typedef CGAL::Projection_on_sphere_traits_3<K>                       Traits;
+typedef CGAL::Delaunay_triangulation_on_sphere_2<Traits>             DToS2;
+
+typedef Traits::Point_3                                              Point_3;
+
+int main(int, char**)
+{
+  std::vector<Point_3> points;
+  points.emplace_back( 3,  1,  1);
+  points.emplace_back(-8,  1,  1);
+  points.emplace_back( 1,  2,  1);
+  points.emplace_back( 1, -2,  1);
+  points.emplace_back( 1,  1, 10);
+
+  Traits traits(Point_3(4,1,1)); // 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;
+    dtos.insert(cst(pt));
+
+    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_solid_faces() << " solid faces" << std::endl;
+    std::cout << dtos.number_of_ghost_faces() << " ghost faces" << std::endl;
+  }
+
+  assert(dtos.is_valid());
+  std::ofstream out("dtos.txt");
+  out.precision(17);
+  out << dtos << std::endl;
+  out.close();
+
+  CGAL::IO::write_OFF("dtos.off", dtos, CGAL::parameters::stream_precision(17));
+
+  DToS2 dtos2;
+  std::ifstream in("dtos.txt");
+  in >> dtos2;
+  in.close();
+  assert(dtos2.is_valid());
+
+  std::cout << "DTOS2 center: " << dtos2.geom_traits().center() << " radius: " << dtos2.geom_traits().radius() << std::endl;
+  std::cout << "DTOS2 has dimension: " << dtos2.dimension() << " and\n";
+  std::cout << dtos2.number_of_vertices() << " vertices" << std::endl;
+  std::cout << dtos2.number_of_edges() << " edges" << std::endl;
+  std::cout << dtos2.number_of_solid_faces() << " solid faces" << std::endl;
+  std::cout << dtos2.number_of_ghost_faces() << " ghost faces" << std::endl;
+
+  CGAL::IO::write_OFF("dtos2.off", dtos2, CGAL::parameters::stream_precision(17));
+
+  assert(dtos.number_of_vertices() == dtos2.number_of_vertices());
+  assert(dtos.number_of_edges() == dtos2.number_of_edges());
+  assert(dtos.number_of_solid_faces() == dtos2.number_of_solid_faces());
+  assert(dtos.number_of_ghost_faces() == dtos2.number_of_ghost_faces());
+
+  std::cout << "Done." << std::endl;
+
+  return EXIT_SUCCESS;
+}