Various fixes for location, (range) insertions

Missing "lt = ..."
Broken spatial sorting
Broken SFINAE for range insertion
Remove ghost edge for dimension 1 (useless, just use the check in locate_edge)
etc.

Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h

@@ -238,6 +238,18 @@ public:
 
 public:
   Vertex_handle insert(const Point& p, Face_handle f = Face_handle());
+
+  // For convenience, when P3 != PoS2
+  template <typename P>
+  Vertex_handle insert(const P& p,
+                       Face_handle f = Face_handle(),
+                       typename std::enable_if<!std::is_same<P, Point>::value>::type* = nullptr)
+  {
+    CGAL_triangulation_assertion((std::is_same<P, Point_3>::value));
+
+    return insert(geom_traits().construct_point_on_sphere_2_object()(p), f);
+  }
+
   Vertex_handle push_back(const Point& p, Face_handle f = Face_handle()) { return insert(p, f); }
   Vertex_handle insert(const Point& p, Locate_type lt, Face_handle loc, int li);
   Vertex_handle insert_first(const Point& p);
@@ -246,37 +258,25 @@ public:
   Vertex_handle insert_outside_affine_hull_regular(const Point& p);
   Vertex_handle insert_cocircular(const Point& p, Locate_type lt, Face_handle loc);
 
-  template <typename InputIterator>
+  // Range insertion
-  int insert(InputIterator first, InputIterator beyond,
-             typename std::enable_if<
-                        std::is_same<typename std::iterator_traits<InputIterator>::value_type,
-                                     Point>::value>::type* = nullptr);
-
-  // For convenience when P3 != PoS2
-  template <typename P>
-  Vertex_handle insert(const P& p,
-                       Face_handle f = Face_handle(),
-                       typename std::enable_if<!std::is_same<P, Point>::value>::type* = nullptr)
-  {
-    return insert(geom_traits().construct_point_on_sphere_2_object()(p), f);
-  }
 
+  // Input range has value type Point, with Point != Point_3
   template <typename InputIterator>
   int insert(InputIterator first, InputIterator beyond,
              typename std::enable_if<
                         !std::is_same<typename std::iterator_traits<InputIterator>::value_type,
-                                      Point>::value>::type* = nullptr)
+                                      Point_3>::value>::type* = nullptr);
-  {
-    typename Geom_traits::Construct_point_on_sphere_2 cst = geom_traits().construct_point_on_sphere_2_object();
 
-    return insert(boost::make_transform_iterator(first, cst),
+  // Input range has value type Point_3, possibly with Point == Point_3
-                  boost::make_transform_iterator(beyond, cst));
+  template <typename InputIterator>
-  }
+  int insert(InputIterator first, InputIterator beyond,
+             typename std::enable_if<
+                        std::is_same<typename std::iterator_traits<InputIterator>::value_type,
+                                     Point_3>::value>::type* = nullptr);
 
-  bool update_ghost_faces(Vertex_handle v = Vertex_handle(), bool first = false);
+  bool update_ghost_faces(Vertex_handle v, bool first = false);
 
-  //REMOVAL
+  // Removal
-  void remove_degree_3(Vertex_handle v, Face_handle f = Face_handle());
   void remove_1D(Vertex_handle v);
   void remove_2D(Vertex_handle v);
   void remove(Vertex_handle v);
@@ -367,36 +367,6 @@ test_conflict(const Point& p, Face_handle fh) const
 }
 
 // ------------------------ INSERTION --------------------------------//
-template <typename Gt, typename Tds>
-typename Delaunay_triangulation_on_sphere_2<Gt, Tds>::Vertex_handle
-Delaunay_triangulation_on_sphere_2<Gt, Tds>::
-insert(const Point& p, Face_handle start)
-{
-  Locate_type lt;
-  int li;
-  Face_handle loc = Base::locate(p, lt, li, start);
-
-  switch(lt)
-  {
-    case NOT_ON_SPHERE:
-      return Vertex_handle();
-    case TOO_CLOSE:
-    {
-      CGAL_assertion(loc != Face_handle());
-      return loc->vertex(li); // @todo test this for dim < 2
-    }
-    case VERTEX:
-    {
-      if(number_of_vertices() == 1)
-        return vertices_begin();
-
-      return (loc->vertex(li));
-    }
-    default: // the point can be inserted
-      return insert(p, lt, loc, li);
-  }
-}
-
 // inserts a new point to a 1d triangulation, the new point is also coplanar with the existing points.
 template <typename Gt, typename Tds>
 typename Delaunay_triangulation_on_sphere_2<Gt, Tds>::Vertex_handle
@@ -426,7 +396,8 @@ insert_cocircular(const Point& p,
 
   tds().delete_face(loc);
 
-  update_ghost_faces(v);
+  CGAL_postcondition(dimension() == 1);
+
   return v;
 }
 
@@ -484,7 +455,6 @@ insert_third(const Point& p)
                                        point(f, 1),
                                        point(f->neighbor(0), 1)) != RIGHT_TURN);
 
-  update_ghost_faces(nv);
   return nv;
 }
 
@@ -579,6 +549,128 @@ insert(const Point& p, Locate_type lt, Face_handle loc, int /*li*/)
   return v;
 }
 
+template <typename Gt, typename Tds>
+typename Delaunay_triangulation_on_sphere_2<Gt, Tds>::Vertex_handle
+Delaunay_triangulation_on_sphere_2<Gt, Tds>::
+insert(const Point& p, Face_handle start)
+{
+  Locate_type lt;
+  int li;
+  Face_handle loc = Base::locate(p, lt, li, start);
+
+  switch(lt)
+  {
+    case NOT_ON_SPHERE:
+      return Vertex_handle();
+    case TOO_CLOSE:
+    {
+      CGAL_triangulation_assertion(loc != Face_handle());
+      return loc->vertex(li);
+    }
+    case VERTEX:
+    {
+      if(number_of_vertices() == 1)
+        return vertices_begin();
+
+      return loc->vertex(li);
+    }
+    default: // the point can be inserted
+      return insert(p, lt, loc, li);
+  }
+}
+
+template <typename Gt, typename Tds>
+template <typename InputIterator>
+int
+Delaunay_triangulation_on_sphere_2<Gt, Tds>::
+insert(InputIterator first, InputIterator beyond,
+       typename std::enable_if<
+                  !std::is_same<typename std::iterator_traits<InputIterator>::value_type,
+                                Point_3>::value>::type*)
+{
+  CGAL_static_assertion((std::is_same<typename std::iterator_traits<InputIterator>::value_type, Point>::value));
+  CGAL_static_assertion(!(std::is_same<Point, Point_3>::value));
+
+  const int n = number_of_vertices();
+
+  // On paper, Spatial_sort_traits_adapter_3 should be used. However, it is not compatible
+  // with spatial_sort_on_sphere() because of the way Transform_coordinates_traits_3 were written:
+  // there are hard calls to x(), y(), and z() whereas it should have been calls to Compute_x_3 & such.
+
+  std::vector<Point> input_points(first, beyond);
+
+  struct Point_3_with_iterator
+    : public Point_3
+  {
+    Point_3_with_iterator(const Point& p, Self const * const tr)
+      : Point_3(tr->geom_traits().construct_point_3_object()(p)), input_point_ptr(&p)
+    { }
+
+    const Point* input_point_ptr;
+  };
+
+  std::vector<Point_3_with_iterator> p3_points;
+  p3_points.reserve(std::distance(first, beyond));
+  for(const Point& p : input_points)
+    p3_points.push_back(Point_3_with_iterator(p, this));
+
+  CGAL::cpp98::random_shuffle(p3_points.begin(), p3_points.end());
+
+  // @todo:
+  // - bench this
+  // - should points not on the sphere already be filtered because it might mess up the sort?
+  spatial_sort_on_sphere(p3_points.begin(), p3_points.end(),
+                         geom_traits(), square(geom_traits().radius()), geom_traits().center());
+
+  // Reordering could be done in-place, but not worth the hassle for now
+  std::vector<Point> points;
+  points.reserve(input_points.size());
+  for(const Point_3_with_iterator& p3wi : p3_points)
+  {
+    CGAL_triangulation_assertion(p3wi.input_point_ptr != nullptr);
+    points.push_back(*(p3wi.input_point_ptr));
+  }
+
+  Face_handle hint;
+  Vertex_handle v;
+  for(const Point& p : points)
+  {
+    v = insert(p, hint);
+    if(v != Vertex_handle()) // can happen e.g. if the point is not on the sphere
+      hint = v->face();
+  }
+
+  return number_of_vertices() - n;
+}
+
+template <typename Gt, typename Tds>
+template <typename InputIterator>
+int
+Delaunay_triangulation_on_sphere_2<Gt, Tds>::
+insert(InputIterator first, InputIterator beyond,
+       typename std::enable_if<
+                  std::is_same<typename std::iterator_traits<InputIterator>::value_type,
+                               Point_3>::value>::type*)
+{
+  const int n = number_of_vertices();
+
+  std::vector<Point_3> points(first, beyond);
+  CGAL::cpp98::random_shuffle(points.begin(), points.end());
+
+  spatial_sort_on_sphere(points.begin(), points.end(), geom_traits(),
+                         square(geom_traits().radius()), geom_traits().center());
+
+  Face_handle hint;
+  for(const Point_3& p : points)
+  {
+    Vertex_handle v = insert(p, hint);
+    if(v != Vertex_handle()) // could happen if the point is not on the sphere
+      hint = v->face();
+  }
+
+  return number_of_vertices() - n;
+}
+
 /*
  * method to test and mark faces incident to v as ghost-faces or solid-faces.
  * the Boolean 'first' defines whether the dimension of the triangulation increased
@@ -589,37 +681,13 @@ bool
 Delaunay_triangulation_on_sphere_2<Gt, Tds>::
 update_ghost_faces(Vertex_handle v, bool first)
 {
-  if(number_of_vertices() < 3)
+  CGAL_triangulation_assertion(dimension() == 2);
-    return false;
-
-  CGAL_assertion(dimension() >= 1);
 
   bool ghost_found = false;
-  if(dimension() == 1)
+
-  {
-    All_edges_iterator eit = all_edges_begin();
-    for(; eit!=all_edges_end(); ++eit)
-    {
-      Face_handle f = eit->first;
-      Face_handle fn = f->neighbor(0);
-      const Point& q = point(fn, 1);
-      if(this->collinear_between(point(f, 0), point(f, 1), q))
-      {
-        f->set_ghost(true);
-        ghost_found = true;
-      }
-      else
-      {
-        f->set_ghost(false);
-      }
-    }
-  }
-  else // dimension==2
-  {
   if(first) // first time dimension 2
   {
-      All_faces_iterator fi = all_faces_begin();
+    for(All_faces_iterator fi=all_faces_begin(); fi!=all_faces_end(); ++fi)
-      for(; fi!=all_faces_end(); fi++)
     {
       if(orientation_on_sphere(fi) != POSITIVE)
       {
@@ -634,7 +702,7 @@ update_ghost_faces(Vertex_handle v, bool first)
   }
   else // not first
   {
-      CGAL_assertion(v != Vertex_handle());
+    CGAL_triangulation_assertion(v != Vertex_handle());
     Face_circulator fc = this->incident_faces(v, v->face());
     Face_circulator done(fc);
     do
@@ -651,90 +719,19 @@ update_ghost_faces(Vertex_handle v, bool first)
     }
     while(++fc != done);
   }
-  }
 
   return ghost_found;
 }
 
-template <typename Gt, typename Tds>
-template <typename InputIterator>
-int
-Delaunay_triangulation_on_sphere_2<Gt, Tds>::
-insert(InputIterator first, InputIterator beyond,
-       typename std::enable_if<
-                  std::is_same<typename std::iterator_traits<InputIterator>::value_type, Point>::value>::type*)
-{
-  int n = number_of_vertices();
-
-  std::vector<Point> points(first, beyond);
-  CGAL::cpp98::random_shuffle(points.begin(), points.end());
-
-  typename Geom_traits::Construct_point_3 cp3 = geom_traits().construct_point_3_object();
-
-  typedef typename Geom_traits::Construct_point_3 Construct_point_3;
-  typedef typename boost::result_of<const Construct_point_3(const Point&)>::type Ret;
-  typedef boost::function_property_map<Construct_point_3,Point, Ret> fpmap;
-  typedef CGAL::Spatial_sort_traits_adapter_3<Geom_traits, fpmap> Search_traits_3;
-
-  // @todo:
-  // - bench this
-  // - should points not on the sphere already be filtered because it might mess up the sort?
-  spatial_sort_on_sphere(points.begin(), points.end(),
-                         Search_traits_3(
-                           boost::make_function_property_map<Point, Ret, Construct_point_3>(
-                             geom_traits().construct_point_3_object()), geom_traits()),
-                         square(geom_traits().radius()), geom_traits().center());
-
-  Face_handle hint;
-  Vertex_handle v;
-  for(typename std::vector<Point>::const_iterator p=points.begin(), end=points.end(); p!=end; ++p)
-  {
-    v = insert(*p, hint);
-    if(v != Vertex_handle()) // could happen if the point is not on the sphere
-      hint = v->face();
-  }
-
-  return number_of_vertices() - n;
-}
-
 //-------------------------------------REMOVAL----------------------------------------------------//
-template <typename Gt, typename Tds>
-void
-Delaunay_triangulation_on_sphere_2<Gt, Tds>::
-remove_degree_3(Vertex_handle v, Face_handle f)
-{
-  CGAL_precondition(v != Vertex_handle());
-
-  if(f == Face_handle())
-    f = v->face();
-
-  tds().remove_degree_3(v, f);
-}
-
-template <typename Gt, typename Tds>
-void
-Delaunay_triangulation_on_sphere_2<Gt, Tds>::
-remove(Vertex_handle v)
-{
-  CGAL_precondition(v != Vertex_handle());
-
-  if(number_of_vertices() <= 3)
-    tds().remove_dim_down(v);
-  else if(dimension() == 2)
-    remove_2D(v);
-  else
-    remove_1D(v);
-}
-
 template <typename Gt, typename Tds>
 void
 Delaunay_triangulation_on_sphere_2<Gt, Tds>::
 remove_1D(Vertex_handle v)
 {
-  CGAL_precondition(v != Vertex_handle());
+  CGAL_triangulation_precondition(v != Vertex_handle());
 
   tds().remove_1D(v);
-  update_ghost_faces();
 }
 
 template <typename Gt, typename Tds>
@@ -747,7 +744,6 @@ remove_2D(Vertex_handle v)
   if(test_dim_down(v)) // resulting triangulation has dim 1
   {
     tds().remove_dim_down(v);
-    update_ghost_faces(); //1d triangulation, no vertex needed to update ghost-faces
   }
   else
   {
@@ -757,6 +753,21 @@ remove_2D(Vertex_handle v)
   }
 }
 
+template <typename Gt, typename Tds>
+void
+Delaunay_triangulation_on_sphere_2<Gt, Tds>::
+remove(Vertex_handle v)
+{
+  CGAL_triangulation_precondition(v != Vertex_handle());
+
+  if(number_of_vertices() <= 3)
+    tds().remove_dim_down(v);
+  else if(dimension() == 2)
+    remove_2D(v);
+  else
+    remove_1D(v);
+}
+
 // tests whether the dimension will decrease when removing v from the triangulation.
 template <typename Gt, typename Tds>
 bool

Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h

@@ -695,22 +695,99 @@ test_distance(const Point& p,
   }
 }
 
+/*
+ * Location for degenerated cases: locates the conflicting edge in a 1-dimensional triangulation.
+ * This method is used when the new point is coplanar with the existing vertices.
+ * The Boolean 'on_diametral_plane' indicates whether the points are also coplanar with the
+ * center of the sphere (true) or not (false).
+ */
+template <class Gt, class Tds>
+typename Triangulation_on_sphere_2<Gt, Tds>::Face_handle
+Triangulation_on_sphere_2<Gt, Tds>::
+locate_edge(const Point& p,
+            Locate_type& lt,
+            int& li,
+            const bool on_diametral_plane) const
+{
+  CGAL_precondition(dimension() == 1);
+
+  Face_handle loc;
+  if(on_diametral_plane)
+  {
+    All_edges_iterator eit;
+    for(eit=all_edges_begin(); eit!=all_edges_end(); ++eit)
+    {
+      const Face_handle f = eit->first;
+      const Face_handle fn = f->neighbor(0);
+      const Point& q = point(fn, 1);
+
+      // Check if *eit is an a "ghost" edge, that is if in the diametral plane, the center
+      // of the sphere is to the right of (eit->source(); eit->target()).
+      // Note that this doesn't check if p is actually "in conflict" with that edge,
+      // so that edge is just kept in memory, and if no solid edge is found, then it has to be this edge
+      if(collinear_between(point(f, 0), point(f, 1), q))
+      {
+        // The new point is on the same 3D plane as the existing vertices, but not on an existing edge
+        loc = eit->first;
+      }
+      else
+      {
+        // not ghost, check the cone
+        if(collinear_between(point(eit->first, 0), point(eit->first, 1), p))
+        {
+          loc = eit->first;
+          lt = EDGE;
+          li = 2;
+          test_distance(p, loc, lt, li);
+          return loc;
+        }
+      }
+    }
+
+    // Couldn't find a solid edge that is in conflict with p
+    lt = OUTSIDE_CONVEX_HULL;
+    li = 4;
+
+    test_distance(p, loc, lt, li);
+    return loc;
+  }
+  else // not coplanar with the center of the sphere
+  {
+    for(All_edges_iterator eit = all_edges_begin(); eit!=all_edges_end(); ++eit)
+    {
+      // The plane going through the center of the sphere, v1, and v2 cuts the non-diametral circle
+      // in two, with the negative side of the plane corresponding to the shorter part, so
+      // if p is coplanar with all the points on that circle and on the negative side, it splits
+      // the edge v1-v2
+      if(orientation_on_sphere(point(eit->first, 0), point(eit->first, 1), p) == ON_NEGATIVE_SIDE)
+      {
+        loc = eit->first;
+        lt = EDGE;
+        li = 2;
+        test_distance(p, loc, lt, li);
+        return loc;
+      }
+    }
+
+    CGAL_assertion(false);
+    return loc;
+  }
+}
+
 template <typename Gt, typename Tds>
 typename Triangulation_on_sphere_2<Gt, Tds>::Face_handle
 Triangulation_on_sphere_2<Gt, Tds>::
 march_locate_1D(const Point& p, Locate_type& lt, int& li) const
 {
-  CGAL_assertion(dimension() == 1);
+  CGAL_triangulation_assertion(dimension() == 1);
 
-
+  // Check if p is coplanar with the existing points first three points of the triangulation
-  // check if p is coplanar with existing points
-  // first three points of triangulation
   Face_handle f = all_edges_begin()->first;
-  Vertex_handle v1 = f->vertex(0);
+  const Vertex_handle v1 = f->vertex(0);
-  Vertex_handle v2 = f->vertex(1);
+  const Vertex_handle v2 = f->vertex(1);
-  Vertex_handle v3 = f->neighbor(0)->vertex(1);
+  const Vertex_handle v3 = f->neighbor(0)->vertex(1);
 
-  Orientation orient = orientation(point(v1), point(v2), point(v3), p);
+  const Orientation orient = orientation(point(v1), point(v2), point(v3), p);
   if(orient != ON_ORIENTED_BOUNDARY)
   {
     lt = OUTSIDE_AFFINE_HULL;
@@ -719,22 +796,23 @@ march_locate_1D(const Point& p, Locate_type& lt, int& li) const
     return f;
   }
 
-  // from then on, p is coplanar with all the triangulation's points
+  // From then on, p is coplanar with all the triangulation's points
 
-  // check if p is coradial with one existing point
+  // Check if p is coradial with one existing point
   Vertices_iterator vi;
   for(vi=vertices_begin(); vi!=vertices_end(); ++vi) // @todo turns insertion into O(n)
   {
     if(are_equal(point(vi), p))
     {
       lt = VERTEX;
-      li = 1;
+      f = vi->face();
+      li = f->index(vi); // could be simply '1'
       return f;
     }
   }
 
   // v1, v2, and v3 are coplanar so this is just checking if they are coplanar with the center of the sphere
-  Orientation pqr = orientation_on_sphere(point(v1), point(v2), point(v3));
+  const Orientation pqr = orientation_on_sphere(point(v1), point(v2), point(v3));
   if(pqr == ON_ORIENTED_BOUNDARY)
     return locate_edge(p, lt, li, true /*on_diametral_plane*/);
   else
@@ -749,7 +827,8 @@ march_locate_2D(Face_handle f,
                 Locate_type& lt,
                 int& li) const
 {
-  CGAL_assertion(!f->is_ghost());
+  CGAL_triangulation_precondition(dimension() == 2);
+  CGAL_triangulation_precondition(!is_ghost(f));
 
   boost::rand48 rng;
   boost::uniform_smallint<> two(0, 1);
@@ -760,7 +839,7 @@ march_locate_2D(Face_handle f,
 
   for(;;)
   {
-    if(f->is_ghost())
+    if(is_ghost(f))
     {
       if(orientation(f, t) == ON_POSITIVE_SIDE) // conflict with the corresponding face
       {
@@ -1004,15 +1083,12 @@ locate(const Point& p,
     return Face_handle();
   }
 
-  //std::cout << "Locate " << p << ", dimension " << dimension() << std::endl;
-
   switch(dimension())
   {
     case -2: // empty triangulation
     {
       lt = OUTSIDE_AFFINE_HULL;
-      li = 4; // li should not be used in this case
+      li = 4;
-
       return Face_handle();
     }
     case -1: // 1 vertex
@@ -1026,7 +1102,6 @@ locate(const Point& p,
           lt = VERTEX;
           Face_handle f = vi->face();
           li = f->index(vi);
-
           return f;
         }
       }
@@ -1039,22 +1114,21 @@ locate(const Point& p,
 
       return f;
     }
-    case 1:
+    case 1: // 3+ coplanar vertices
     {
       return march_locate_1D(p, lt, li);
     }
   }
 
   // below is dimension() == 2
-
   if(start == Face_handle())
     start = all_faces_begin();
 
-  if(start->is_ghost())
+  if(is_ghost(start))
   {
-    for(All_faces_iterator it = this->_tds.face_iterator_base_begin(); it !=all_faces_end(); it++)
+    for(All_faces_iterator it=all_faces_begin(); it !=all_faces_end(); ++it)
     {
-      if(!it->is_ghost())
+      if(!is_ghost(it))
       {
         start = it;
         break;
@@ -1122,8 +1196,7 @@ show_all() const
   if(dimension() == 1)
   {
     std::cerr << " all edges dim 1 " << std::endl;
-    All_edges_iterator aeit;
+    for(All_edges_iterator aeit=all_edges_begin(); aeit!=all_edges_end(); ++aeit)
-    for(aeit=all_edges_begin(); aeit!=all_edges_end(); ++aeit)
     {
       show_face(aeit->first);
       std::cerr << "   ------------   " << std::endl;
@@ -1133,8 +1206,7 @@ show_all() const
   }
 
   std::cerr << " faces " << std::endl;
-  All_faces_iterator fi;
+  for(All_faces_iterator fi = all_faces_begin(); fi !=all_faces_end(); ++fi)
-  for(fi = all_faces_begin(); fi !=all_faces_end(); ++fi)
   {
     show_face(fi);
     std::cerr << "   ------------   " << std::endl;
@@ -1143,8 +1215,7 @@ show_all() const
   if(number_of_vertices() > 1)
   {
     std::cerr << "print triangulation vertices:" << std::endl;
-    Vertices_iterator vi;
+    for(Vertices_iterator vi=vertices_begin(); vi!=vertices_end(); ++vi)
-    for(vi=vertices_begin(); vi!=vertices_end(); ++vi)
     {
       show_vertex(vi);
       std::cerr << "  / associated face: " << (void*)(&(*(vi->face()))) << std::endl;;