From c40d7124d07780fcf8d64d0a8529cc8ebc3e4ad7 Mon Sep 17 00:00:00 2001
From: Laurent Rineau <laurent.rineau@cgal.org>
Date: Fri, 18 May 2018 17:47:01 +0200
Subject: [PATCH] Fix CDT_2 errors, using snapping of intersection points
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

When `Constrained_triangulation_2` computes the intersection of two
segments, with a floating-point number type, and with
`Exact_predicates_tag`, the computed intersection point is snapped to an
extremity of the two segments, if it is closest to 4 ulp (with the l-inf
distance).

That value `4` can be changed by defining the macro
`CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE` to another value.
---
 .../CGAL/Constrained_triangulation_2.h        | 64 ++++++++++++++++++-
 1 file changed, 63 insertions(+), 1 deletion(-)

diff --git a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
index 6ca828a3adb..6f60a7ec484 100644
--- a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
+++ b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h
@@ -37,6 +37,9 @@
 
 #include <boost/mpl/if.hpp>
 #include <boost/iterator/filter_iterator.hpp>
+#include <boost/math/special_functions/next.hpp>
+#include <boost/type_traits/is_floating_point.hpp>
+
 namespace CGAL {
 
 struct No_intersection_tag{};
@@ -1422,11 +1425,70 @@ 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;
+  }
+
+  using boost::math::float_advance;
+#ifdef CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE
+  const int dist = CGAL_CDT_2_INTERSECTION_SNAPPING_ULP_DISTANCE;
+#else
+  const int dist = 4;
+#endif
+  const Bbox_2 bbox(float_advance(pi.x(), -dist), float_advance(pi.y(), -dist),
+                    float_advance(pi.x(), +dist), float_advance(pi.y(), +dist));
+  if(do_overlap(bbox, pa.bbox())) pi = pa;
+  if(do_overlap(bbox, pb.bbox())) pi = pb;
+  if(do_overlap(bbox, pc.bbox())) pi = pc;
+  if(do_overlap(bbox, pd.bbox())) 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