Merge pull request #4735 from lrineau/Triangulation_2-fix_CDT_2_degenerate_case-GF

Fix a bug in <CGAL/Constrained_triangulation_2.h>

STL_Extension/include/CGAL/Compact_container.h

@@ -1119,7 +1119,7 @@ namespace handle {
   template <class H> struct Hash_functor;
 
   template<class DSC, bool Const>
-  struct Hash_functor<CC_iterator<DSC, Const>>{
+  struct Hash_functor<CC_iterator<DSC, Const> >{
     std::size_t
     operator()(const CC_iterator<DSC, Const>& i)
     {

Triangulation_2/include/CGAL/Constrained_triangulation_2.h

@@ -51,6 +51,25 @@ struct CGAL_DEPRECATED No_intersection_tag :
 
 namespace internal {
 
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  struct Indentation_level {
+    int n;
+    Indentation_level() : n(0) {}
+    friend std::ostream& operator<<(std::ostream& os, Indentation_level level) {
+      return os << std::string(2*level.n, ' ');
+    }
+    Indentation_level& operator++() { ++n; return *this; }
+    Indentation_level& operator--() { --n; return *this; }
+    struct Exit_guard {
+      Exit_guard(Indentation_level& level): level(level) { ++level; }
+      Exit_guard(const Exit_guard& other) : level(other.level) { ++level; }
+      Indentation_level& level;
+      ~Exit_guard() { --level; }
+    };
+    Exit_guard open_new_scope() { return Exit_guard(*this); }
+  } cdt_2_indent_level;
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
+
 template <typename K>
 struct Itag {
   typedef typename boost::mpl::if_<typename Algebraic_structure_traits<typename K::FT>::Is_exact,
@@ -728,10 +747,25 @@ insert_constraint(Vertex_handle  vaa, Vertex_handle vbb)
   std::stack<std::pair<Vertex_handle, Vertex_handle> > stack;
   stack.push(std::make_pair(vaa,vbb));
 
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::insert_constraint( #" << vaa->time_stamp() << "= " << vaa->point()
+            << " , #" << vbb->time_stamp() << "= " << vbb->point()
+            << " )\n";
+  internal::Indentation_level::Exit_guard exit_guard = CGAL::internal::cdt_2_indent_level.open_new_scope();
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
   while(! stack.empty()){
     boost::tie(vaa,vbb) = stack.top();
     stack.pop();
     CGAL_triangulation_precondition( vaa != vbb);
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+    std::cerr << CGAL::internal::cdt_2_indent_level
+              << "CT_2::insert_constraint, stack pop=( #" << vaa->time_stamp() << "= " << vaa->point()
+              << " , #" << vbb->time_stamp() << "= " << vbb->point()
+              << " ) remaining stack size: "
+              << stack.size() << '\n';
+    CGAL_assertion(this->is_valid());
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
     Vertex_handle vi;
 
     Face_handle fr;
@@ -767,10 +801,26 @@ insert_constraint(Vertex_handle  vaa, Vertex_handle vbb)
                                                  vi);
     if ( intersection) {
       if (vi != vaa && vi != vbb) {
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::insert_constraint stask push [vaa, vi] ( #" << vaa->time_stamp() << "= " << vaa->point()
+            << " , #" << vi->time_stamp() << "= " << vi->point()
+            << " )\n";
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::insert_constraint stask push [vi, vbb] ( #" << vi->time_stamp() << "= " << vi->point()
+            << " , #" << vbb->time_stamp() << "= " << vbb->point()
+            << " )\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
         stack.push(std::make_pair(vaa,vi));
         stack.push(std::make_pair(vi,vbb));
       }
       else{
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::insert_constraint stask push [vaa, vbb]( #" << vaa->time_stamp() << "= " << vaa->point()
+            << " , #" << vbb->time_stamp() << "= " << vbb->point()
+            << " )\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
         stack.push(std::make_pair(vaa,vbb));
       }
       continue;
@@ -814,6 +864,30 @@ find_intersected_faces(Vertex_handle vaa,
 
   // to deal with the case where the first crossed edge
   // is constrained
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::find_intersected_faces ( #" << vaa->time_stamp() << "= " << vaa->point()
+            << " , #" << vbb->time_stamp() << "= " << vbb->point()
+            << " )\n"
+            << CGAL::internal::cdt_2_indent_level
+            << "> current constrained edges are:\n";
+  for(Constrained_edges_iterator edge_it = this->constrained_edges_begin(),
+        end = this->constrained_edges_end();
+      edge_it != end; ++edge_it)
+  {
+    std::cerr <<CGAL::internal::cdt_2_indent_level
+              << "> (#"
+              << edge_it->first->vertex(cw(edge_it->second))->time_stamp()
+              << ", #"
+              << edge_it->first->vertex(ccw(edge_it->second))->time_stamp()
+              << ")\n";
+  }
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "> current face is ( #" << current_face->vertex(0)->time_stamp()
+            << " #" << current_face->vertex(1)->time_stamp()
+            << " #" << current_face->vertex(2)->time_stamp() << " )\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
+
   if(current_face->is_constrained(ind)) {
     vi=intersect(current_face, ind, vaa, vbb);
     return true;
@@ -971,14 +1045,22 @@ intersect(Face_handle f, int i,
   const Point& pc = vcc->point();
   const Point& pd = vdd->point();
 
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::intersect segment ( #" << vaa->time_stamp() << "= " << vaa->point()
+            << " , #" << vbb->time_stamp() << "= " << vbb->point()
+            << " ) with edge ( #"<< vcc->time_stamp() << "= " << vcc->point()
+            << " , #" << vdd->time_stamp() << "= " << vdd->point()
+            << " )\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
   Point pi; //creator for point is required here
   Itag itag = Itag();
   bool ok  = intersection(geom_traits(), pa, pb, pc, pd, pi, itag );
 
   Vertex_handle vi;
   if ( !ok) {  //intersection detected but not computed
-    int i = limit_intersection(geom_traits(), pa, pb, pc, pd, itag);
-    switch(i){
+    int int_index = limit_intersection(geom_traits(), pa, pb, pc, pd, itag);
+    switch(int_index){
     case 0 : vi = vaa; break;
     case 1 : vi = vbb; break;
     case 2 : vi = vcc; break;
@@ -992,6 +1074,11 @@ intersect(Face_handle f, int i,
     remove_constrained_edge(f, i);
     vi = virtual_insert(pi, f);
   }
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::intersect, `vi` is ( #" << vi->time_stamp() << "= " << vi->point()
+            << " )\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
 
   // vi == vc or vi == vd may happen even if intersection==true
   // due to approximate construction of the intersection
@@ -1247,6 +1334,14 @@ void
 Constrained_triangulation_2<Gt,Tds,Itag>::
 remove_constrained_edge(Face_handle f, int i)
 {
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  std::cerr << CGAL::internal::cdt_2_indent_level
+            << "CT_2::remove_constrained_edge ( #"
+            << f->vertex(cw(i))->time_stamp()
+            << ", #"
+            << f->vertex(ccw(i))->time_stamp()
+            << ")\n";
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
   f->set_constraint(i, false);
   if (dimension() == 2)
     (f->neighbor(i))->set_constraint(mirror_index(f,i), false);
@@ -1518,7 +1613,8 @@ intersection(const Gt& gt,
   if(!result) return result;
   if(pi == pa || pi == pb || pi == pc || pi == pd) {
 #ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
-    std::cerr << "  CT_2::intersection: intersection is an existing point "
+    std::cerr << CGAL::internal::cdt_2_indent_level
+              << "  CT_2::intersection: intersection is an existing point "
               << pi << std::endl;
 #endif
     return result;
@@ -1540,7 +1636,8 @@ intersection(const Gt& gt,
   if(do_overlap(bb, bbox(pd))) pi = pd;
 #ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
   if(pi == pa || pi == pb || pi == pc || pi == pd) {
-    std::cerr << "  CT_2::intersection: intersection SNAPPED to an existing point "
+    std::cerr << CGAL::internal::cdt_2_indent_level
+              << "  CT_2::intersection: intersection SNAPPED to an existing point "
               << pi << std::endl;
   }
 #endif
@@ -1578,6 +1675,17 @@ compute_intersection(const Gt& gt,
     construct_segment=gt.construct_segment_2_object();
   Object result = compute_intersec(construct_segment(pa,pb),
                                    construct_segment(pc,pd));
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+  typename Gt::Segment_2 s;
+  if(assign(s, result)) {
+    std::cerr << CGAL::internal::cdt_2_indent_level
+              << "compute_intersection: " << s << '\n';
+  }
+  if(assign(pi, result)) {
+    std::cerr << CGAL::internal::cdt_2_indent_level
+              << "compute_intersection: " << pi << '\n';
+  }
+#endif // CGAL_CDT_2_DEBUG_INTERSECTIONS
   return assign(pi, result);
 }
 

Triangulation_2/test/Triangulation_2/issue_4405.cpp

@@ -9,13 +9,43 @@ typedef CGAL::Epick Kernel;
 typedef Kernel::FT FieldNumberType;
 typedef Kernel::Point_2 Point2;
 typedef Kernel::Point_3 Point3;
+
+template <typename Vb>
+class My_vertex_base : public Vb {
+  std::size_t time_stamp_;
+public:
+  My_vertex_base() : Vb(), time_stamp_(-1) {
+  }
+
+  My_vertex_base(const My_vertex_base& other) :
+    Vb(other),
+    time_stamp_(other.time_stamp_)
+  {}
+
+  typedef CGAL::Tag_true Has_timestamp;
+
+  std::size_t time_stamp() const {
+    return time_stamp_;
+  }
+  void set_time_stamp(const std::size_t& ts) {
+    time_stamp_ = ts;
+  }
+
+  template < class TDS >
+  struct Rebind_TDS {
+    typedef typename Vb::template Rebind_TDS<TDS>::Other Vb2;
+    typedef My_vertex_base<Vb2> Other;
+  };
+};
+
 struct FaceInfo2
 {
   unsigned long long m_id;
 };
 typedef CGAL::Projection_traits_xy_3<Kernel> TriangulationTraits;
 typedef CGAL::Triangulation_vertex_base_with_id_2<TriangulationTraits> VertexBaseWithId;
-typedef CGAL::Triangulation_vertex_base_2<TriangulationTraits, VertexBaseWithId> VertexBase;
+typedef My_vertex_base<VertexBaseWithId> Vb2;
+typedef CGAL::Triangulation_vertex_base_2<TriangulationTraits, Vb2> VertexBase;
 typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, Kernel> FaceBaseWithInfo;
 typedef CGAL::Constrained_triangulation_face_base_2<TriangulationTraits, FaceBaseWithInfo> FaceBase;
 typedef CGAL::Triangulation_data_structure_2<VertexBase, FaceBase> TriangulationData;

Triangulation_2/test/Triangulation_2/test_cdt_degenerate_case.cpp

@@ -1,3 +1,4 @@
+#define CGAL_CDT_2_DEBUG_INTERSECTIONS 1
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Constrained_Delaunay_triangulation_2.h>
 #include <CGAL/Constrained_triangulation_plus_2.h>
@@ -5,7 +6,41 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel EPIC;
 typedef EPIC::Point_2 Point_2;
-typedef CGAL::Triangulation_vertex_base_2<EPIC> Vb;
+
+template <typename Vb>
+class My_vertex_base : public Vb {
+  std::size_t time_stamp_;
+public:
+  My_vertex_base() : Vb(), time_stamp_(-1) {
+  }
+
+  My_vertex_base(const My_vertex_base& other) :
+    Vb(other),
+    time_stamp_(other.time_stamp_)
+  {}
+
+  typedef CGAL::Tag_true Has_timestamp;
+
+  std::size_t time_stamp() const {
+    return time_stamp_;
+  }
+  void set_time_stamp(const std::size_t& ts) {
+    time_stamp_ = ts;
+  }
+
+  template < class TDS >
+  struct Rebind_TDS {
+    typedef typename Vb::template Rebind_TDS<TDS>::Other Vb2;
+    typedef My_vertex_base<Vb2> Other;
+  };
+};
+
+#ifdef CGAL_CDT_2_DEBUG_INTERSECTIONS
+typedef My_vertex_base<CGAL::Triangulation_vertex_base_2<EPIC> > Vb;
+#else
+typedef CGAL::Triangulation_vertex_base_2<K> Vb;
+#endif
+
 typedef CGAL::Constrained_triangulation_face_base_2<EPIC> Fb;
 typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
 typedef CGAL::Exact_predicates_tag Itag;
@@ -14,6 +49,7 @@ typedef CGAL::Constrained_triangulation_plus_2<CDT> CDTp2;
 
 template <class CDT>
 void test() {
+  std::cerr.precision(17);
   CDT cdt;
   cdt.insert_constraint(Point_2(  48.0923419883269,   299.7232779774145  ),
                         Point_2(  66.05373710316852,  434.231770798343   ));