Merge pull request #3098 from lrineau/Triangulation_2-Debug_CDT2-lrineau

Fix CDT_2 errors, using snapping of intersection points

Kernel_23/doc/Kernel_23/CGAL/Bbox_2.h

@@ -95,6 +95,11 @@ updates `b` to be the bounding box of `b` and `c` and returns itself.
 */
 Bbox_2& operator+=(const Bbox_2 &c);
 
+/*!
+dilates the bounding box by a specified number of ULP.
+*/
+void dilate(int dist);
+  
 /// @}
 
 }; /* end Bbox_2 */

Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h

@@ -109,6 +109,10 @@ updates `b` to be the bounding box of `b` and `c` and returns itself.
 */
 Bbox_3& operator+=(const Bbox_3 &c);
 
+/*!
+dilates the bounding box by a specified number of ULP.
+*/
+void dilate(int dist);
 /// @}
 
 }; /* end Bbox_3 */

Kernel_23/include/CGAL/Bbox_2.h

@@ -30,6 +30,7 @@
 #include <CGAL/IO/io.h>
 #include <CGAL/Dimension.h>
 #include <CGAL/array.h>
+#include <boost/math/special_functions/next.hpp>
 
 namespace CGAL {
 
@@ -76,6 +77,7 @@ public:
   inline Bbox_2     operator+(const Bbox_2 &b) const;
   inline Bbox_2&     operator+=(const Bbox_2 &b);
 
+  inline void dilate(int dist);
 };
 
 inline
@@ -156,7 +158,17 @@ Bbox_2::operator+=(const Bbox_2& b)
   rep[3] = (std::max)(ymax(), b.ymax());
   return *this;
 }
-
+inline
+void
+Bbox_2::dilate(int dist)
+{
+  using boost::math::float_advance;
+  float_advance(rep[0],-dist);
+  float_advance(rep[1],-dist);
+  float_advance(rep[2],dist);
+  float_advance(rep[3],dist);
+}
+  
 inline
 bool
 do_overlap(const Bbox_2 &bb1, const Bbox_2 &bb2)

Kernel_23/include/CGAL/Bbox_3.h

@@ -31,6 +31,7 @@
 #include <CGAL/IO/io.h>
 #include <CGAL/Dimension.h>
 #include <CGAL/array.h>
+#include <boost/math/special_functions/next.hpp>
 
 namespace CGAL {
 
@@ -78,6 +79,8 @@ public:
 
   Bbox_3  operator+(const Bbox_3& b) const;
   Bbox_3& operator+=(const Bbox_3& b);
+
+  void dilate(int dist);
 };
 
 inline
@@ -176,6 +179,20 @@ Bbox_3::operator+=(const Bbox_3& b)
   return *this;
 }
 
+inline
+void
+Bbox_3::dilate(int dist)
+{
+  using boost::math::float_advance;
+  float_advance(rep[0],-dist);
+  float_advance(rep[1],-dist);
+  float_advance(rep[2],-dist);
+  float_advance(rep[3],dist);
+  float_advance(rep[4],dist);
+  float_advance(rep[5],dist);
+}
+
+
 inline
 bool
 do_overlap(const Bbox_3& bb1, const Bbox_3& bb2)

Kernel_23/include/CGAL/internal/Projection_traits_3.h

@@ -805,6 +805,7 @@ public:
   typedef typename Rp::Construct_scaled_vector_3              Construct_scaled_vector_2;
   typedef typename Rp::Construct_triangle_3                   Construct_triangle_2;
   typedef typename Rp::Construct_line_3                       Construct_line_2;
+  typedef typename Rp::Construct_bbox_3                       Construct_bbox_2;
 
   struct Less_xy_2 {
     typedef bool result_type;
@@ -989,6 +990,9 @@ public:
   Construct_line_2  construct_line_2_object() const
     {return Construct_line_2();}
 
+  Construct_bbox_2  construct_bbox_2_object() const
+    {return Construct_bbox_2();}
+
   Compute_scalar_product_2 compute_scalar_product_2_object() const
     {return Compute_scalar_product_2();}
 

Triangulation_2/doc/Triangulation_2/Concepts/ConstrainedTriangulationTraits_2.h

@@ -76,6 +76,15 @@ between `p` and `l`.
 */ 
 typedef unspecified_type Compute_squared_distance_2; 
 
+/*!
+A function object whose
+`operator()` computes the bounding box of a point.
+
+`unspecified_type operator()(Point_2 p);` Returns the bounding box of `p`.
+The result type is either `Bbox_2` or `Bbox_3` (for projection traits classes).
+*/ 
+typedef unspecified_type Compute_bounding_box_2; 
+
 /// @} 
 
 /// \name Access to Constructor Objects 

Triangulation_2/include/CGAL/Constrained_triangulation_2.h

@@ -40,6 +40,10 @@
 
 #include <boost/mpl/if.hpp>
 #include <boost/iterator/filter_iterator.hpp>
+
+#include <boost/utility/result_of.hpp>
+#include <boost/type_traits/is_floating_point.hpp>
+
 namespace CGAL {
 
 struct No_intersection_tag{};
@@ -1440,11 +1444,72 @@ intersection(const Gt& gt,
 	     const typename Gt::Point_2& pc, 
 	     const typename Gt::Point_2& pd,
 	     typename Gt::Point_2& pi,
-	     Exact_predicates_tag)
+	     Exact_predicates_tag,
+	     CGAL::Tag_false /* not a FT is not floating-point */)
 {
   return compute_intersection(gt,pa,pb,pc,pd,pi);
 }
 
+template<class Gt>
+inline bool
+intersection(const Gt& gt,
+	     const typename Gt::Point_2& pa,
+	     const typename Gt::Point_2& pb,
+	     const typename Gt::Point_2& pc,
+	     const typename Gt::Point_2& pd,
+	     typename Gt::Point_2& pi,
+	     Exact_predicates_tag,
+	     CGAL::Tag_true /* FT is a floating-point type */)
+{
+  const bool result = compute_intersection(gt,pa,pb,pc,pd,pi);
+  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 "
+              << pi << std::endl;
+#endif
+    return result;
+  }
+
+
+#ifdef CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE
+  const int dist = CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE;
+#else
+  const int dist = 4;
+#endif
+  typedef typename Gt::Construct_bbox_2 Construct_bbox_2;
+  Construct_bbox_2 bbox = gt.construct_bbox_2_object();
+  typename boost::result_of<const Construct_bbox_2(const typename Gt::Point_2&)>::type bb(bbox(pi));
+  bb.dilate(dist);
+  if(do_overlap(bb, bbox(pa))) pi = pa;
+  if(do_overlap(bb, bbox(pb))) pi = pb;
+  if(do_overlap(bb, bbox(pc))) pi = pc;
+  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 "
+              << pi << std::endl;
+  }
+#endif
+  return result;
+}
+
+template<class Gt>
+inline bool
+intersection(const Gt& gt,
+	     const typename Gt::Point_2& pa,
+	     const typename Gt::Point_2& pb,
+	     const typename Gt::Point_2& pc,
+	     const typename Gt::Point_2& pd,
+	     typename Gt::Point_2& pi,
+	     Exact_predicates_tag exact_predicates_tag)
+{
+  typedef typename Gt::FT FT;
+  return intersection(gt,pa,pb,pc,pd,pi,
+                      exact_predicates_tag,
+                      Boolean_tag<boost::is_floating_point<FT>::value>());
+}
+
 
 template<class Gt>
 bool

Triangulation_2/include/CGAL/internal/Triangulation_2_projection_traits_base_3.h

@@ -380,6 +380,7 @@ public:
   typedef typename K::Construct_circumcenter_3      Construct_circumcenter_2;
 
   typedef typename K::Compute_area_3                Compute_area_2;
+  typedef typename K::Construct_bbox_3              Construct_bbox_2;
 
   Less_x_2
   less_x_2_object() const
@@ -466,6 +467,9 @@ public:
   {return Compute_area_2();}
 
 
+  Construct_bbox_2  construct_bbox_2_object() const
+    {return Construct_bbox_2();}
+
 
   // Special functor, not in the Kernel concept
   class Projection_to_plan {