Mesh_2: fix bug with Mesh_2/test/Mesh_2/triwild-10070-min1.obj

The bug was reproduced by

```shell
/path/to/build/test/Mesh_2/test_mesh_obj Mesh_2/test/Mesh_2/triwild-10070-min1.obj
```

The log of the error, with all the debug options of Mesh_2 was:

```plain
Before conforming Gabriel: 4 vertices.
  Smallest squared distance between constraint endpoints: 0.010522898628224919
  edge #1= 453.57378188125 445.0501073328125 -- #3= 454.05836694375 445.48749872343745 is encroached by #2= 453.77453855000005 445.23131233749996
  add_constrained_edge_to_be_conformed(#1= 453.57378188125 445.0501073328125, #3= 454.05836694375 445.48749872343745)
  edge #2= 453.77453855000005 445.23131233749996 -- #0= 453.49763280000002 444.98137439999999 is encroached by #1= 453.57378188125 445.0501073328125
  add_constrained_edge_to_be_conformed(#2= 453.77453855000005 445.23131233749996, #0= 453.49763280000002 444.98137439999999)
split_cluster_point(454.05836694375 445.48749872343745 , 453.57378188125 445.0501073328125)
  reduced: 0
result:  453.77453855000005 445.23131233750001
  edge #1= 453.57378188125 445.0501073328125 -- #3= 454.05836694375 445.48749872343745 is encroached by #2= 453.77453855000005 445.23131233749996
Refine_edges_with_clusters::refinement_point_impl(Edge: (#1= 453.57378188125 445.0501073328125, #3= 454.05836694375 445.48749872343745) = 453.77453855000005 445.23131233750001
(453.77453855000005 445.23131233750001) accepted
on edge, insert(453.77453855000005 445.23131233750001): 3 boundary edges in the zone
 inserted new vertex #4= 453.77453855000005 445.23131233750001
update_clusters
va_has_a_cluster=0
vb_has_a_cluster=1
clusters.size()=2
E  edge #3= 454.05836694375 445.48749872343745 -- #2= 453.77453855000005 445.23131233749996 is encroached by #4= 453.77453855000005 445.23131233750001
  add_constrained_edge_to_be_conformed(#3= 454.05836694375 445.48749872343745, #2= 453.77453855000005 445.23131233749996)
  edge #1= 453.57378188125 445.0501073328125 -- #4= 453.77453855000005 445.23131233750001 is encroached by #2= 453.77453855000005 445.23131233749996
  add_constrained_edge_to_be_conformed(#1= 453.57378188125 445.0501073328125, #4= 453.77453855000005 445.23131233750001)
Cluster at #3= 454.05836694375 445.48749872343745 is updated.
  vm: #4= 453.77453855000005 445.23131233750001
  reduction: 1
  min_sq_len: 0.1461900214383674
clusters.size() after update_cluster=2
split_cluster_point(453.49763280000002 444.98137439999999 , 453.77453855000005 445.23131233749996)
  reduced: 0
result:  453.64993096249998 445.11884026562501
  edge #2= 453.77453855000005 445.23131233749996 -- #0= 453.49763280000002 444.98137439999999 is encroached by #1= 453.57378188125 445.0501073328125
Refine_edges_with_clusters::refinement_point_impl(Edge: (#2= 453.77453855000005 445.23131233749996, #0= 453.49763280000002 444.98137439999999) = 453.64993096249998 445.11884026562501
(453.64993096249998 445.11884026562501) accepted
on edge, insert(453.64993096249998 445.11884026562501): 0 boundary edges in the zone
terminate called after throwing an instance of 'CGAL::Precondition_exception'
  what():  CGAL ERROR: precondition violation!
Expr: i >= 0 && i < 3
```

The following edge `(#1, #3)` is encroached by the vertex `#2`:

```plain
  edge #1= 453.57378188125 445.0501073328125 -- #3= 454.05836694375 445.48749872343745 is encroached by #2= 453.77453855000005 445.23131233749996
```

The inserted vertex `#4= 453.77453855000005 445.23131233750001` is the projection of the vertex `#2= 453.77453855000005 445.23131233749996` (rules to refine a segment part of a cluster.

The distance between `#2` and its projection `#4` is:

- `0` on the x-axis,
- `ulp(v2.y()))` on the y-axis.

... and then, Mesh_2 cannot do anything with that.

The fix is to implement a snapping strategy: if the encroaching vertex is close enough, then chose it be the refinement point of the encroached segment. "Close enough" means:

- create tiny `Bbox_2` centered on the encroaching vertex, with snaps equal to 8 ulp on both axis,
- if that tiny bbox intersects the encroached segment, then the refinement is the encroaching vertex.

Author: lrineau

Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h

@@ -12,7 +12,10 @@
 #ifndef CGAL_INTERNAL_PROJECTION_TRAITS_3_H
 #define CGAL_INTERNAL_PROJECTION_TRAITS_3_H
 
+#include <CGAL/Bbox_2.h>
+#include <CGAL/Bbox_3.h>
 #include <CGAL/assertions.h>
+#include <CGAL/enum.h>
 #include <CGAL/tags.h>
 #include <CGAL/Point_3.h>
 #include <CGAL/Segment_3.h>
@@ -22,6 +25,11 @@
 #include <CGAL/Kernel_23/internal/Has_boolean_tags.h>
 #include <CGAL/Triangulation_structural_filtering_traits.h>
 
+#include <array>
+#include <limits>
+#include <optional>
+#include <variant>
+
 namespace CGAL {
 
 namespace internal {
@@ -409,6 +417,24 @@ class Compute_determinant_projected_3
   }
 };
 
+template <class R,int dim>
+struct Do_intersect_projected_3
+{
+  template <typename T> bool operator()(const Bbox_2& bbox, const T& obj) const {
+    static constexpr double infinity = std::numeric_limits<double>::infinity();
+    std::array<double, 3> bbox_min = {-infinity, -infinity, -infinity};
+    std::array<double, 3> bbox_max = {infinity, infinity, infinity};
+
+    using Proj = Projector<R, dim>;
+    bbox_min[Proj::x_index] = bbox.xmin();
+    bbox_min[Proj::y_index] = bbox.ymin();
+    bbox_max[Proj::x_index] = bbox.xmax();
+    bbox_max[Proj::y_index] = bbox.ymax();
+    Bbox_3 bbox3(bbox_min[0], bbox_min[1], bbox_min[2], bbox_max[0], bbox_max[1], bbox_max[2]);
+    return do_intersect(bbox3, obj);
+  }
+};
+
 template <class R,int dim>
 class  Intersect_projected_3
 {
@@ -951,6 +977,7 @@ class Projection_traits_3
   typedef Collinear_are_ordered_along_line_projected_3<Rp,dim> Collinear_are_ordered_along_line_2;
   typedef Squared_distance_projected_3<Rp,dim>                Compute_squared_distance_2;
   typedef Intersect_projected_3<Rp,dim>                       Intersect_2;
+  typedef Do_intersect_projected_3<Rp,dim>                    Do_intersect_2;
   typedef Compute_squared_radius_projected<Rp,dim>            Compute_squared_radius_2;
   typedef Compute_scalar_product_projected_3<Rp,dim>          Compute_scalar_product_2;
   typedef Compute_squared_length_projected_3<Rp,dim>          Compute_squared_length_2;
@@ -1151,6 +1178,12 @@ class Projection_traits_3
     return Intersect_2();
   }
 
+  Do_intersect_2
+  do_intersect_2_object() const
+  {
+    return Do_intersect_2();
+  }
+
   Construct_point_2 construct_point_2_object() const
   { return Construct_point_2();}
 

Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_base_3.h

@@ -503,6 +503,7 @@ class Projection_traits_base_3
 
   typedef typename K::Compute_area_3                Compute_area_2;
   typedef typename K::Construct_bbox_3              Construct_bbox_2;
+  typedef typename K::Do_intersect_3                Do_intersect_2; // for do_intersect(bbox_2, segment_2)
 
   Less_x_2
   less_x_2_object() const
@@ -607,6 +608,9 @@ class Projection_traits_base_3
   Construct_bbox_2  construct_bbox_2_object() const
     {return Construct_bbox_2();}
 
+  Do_intersect_2 do_intersect_2_object() const
+    {return Do_intersect_2();}
+
 
   // Special functor, not in the Kernel concept
   class Projection_to_plan {

Mesh_2/include/CGAL/Mesh_2/Refine_edges.h

@@ -29,10 +29,9 @@
 #include <CGAL/tags.h>
 
 #include <boost/iterator/filter_iterator.hpp>
-
 #include <boost/iterator/transform_iterator.hpp>
+
 #include <iterator>
-#include <limits>
 #include <optional>
 #include <stack>
 #include <type_traits>
@@ -625,43 +624,68 @@ class Refine_edges_base :
               << IO::oformat(vb, With_point_tag{}) << ") = ";
     std::cerr << p << '\n';
 #endif // CGAL_MESH_2_DEBUG_REFINEMENT_POINTS
+    return maybe_snap_to_existing_vertex(edge, p);
+  }
+
+  auto get_closest_encroaching_vertex(Edge e) const {
+    struct Result {
+      Vertex_handle vertex;
+      FT squared_distance;
+    };
+    std::optional<Result> result;
+    Is_locally_conforming_Gabriel<Tr> is_locally_conforming_Gabriel{};
+
+    for(int i = 0; i < 2; ++i, e = tr.mirror_edge(e)) {
+      auto [c, index] = e;
+      if(tr.is_infinite(c))
+        continue;
+      auto v = c->vertex(index);
+      auto p = v->point();
+      if(is_locally_conforming_Gabriel(tr, c, index, p, v))
+        continue;
+      auto va = e.first->vertex(tr.cw (e.second));
+      auto vb = e.first->vertex(tr.ccw(e.second));
+      typename Tr::Segment ab(va->point(), vb->point());
+      auto height_squared = squared_distance(p, ab);
+      if(!result || height_squared < result->squared_distance) {
+        result = Result{c->vertex(index), height_squared};
+      }
+    }
+    return result;
+  }
+
+  Point maybe_snap_to_existing_vertex(const Edge& edge, Point p) const {
     if constexpr(std::is_floating_point_v<FT>) {
-      Is_locally_conforming_Gabriel<Tr> is_locally_conforming_Gabriel{};
-      constexpr FT max_FT = (std::numeric_limits<FT>::max)();
-      auto ulp = [&](const FT& x) { return std::nextafter(x, max_FT) - x; };
-      auto sq_delta = 256 * CGAL::square((std::min)(ulp(p.x()), ulp(p.y())));
-
-      typename Tr::Segment ab(this->va->point(), this->vb->point());
-      std::optional<FT> min_height_squared;
-      Vertex_handle closest_encroaching_vertex;
-      auto e = edge;
-      for(int i = 0; i < 2; ++i, e = triangulation_ref_impl().mirror_edge(e)) {
-        auto [c, index] = e;
-        if(triangulation_ref_impl().is_infinite(c))
-          continue;
-        auto cp = c->vertex(index)->point();
-        if(is_locally_conforming_Gabriel(triangulation_ref_impl(), c, index, cp))
-          continue;
-        auto height_squared = squared_distance(cp, ab);
-        if(!min_height_squared || height_squared < *min_height_squared) {
-          closest_encroaching_vertex = c->vertex(index);
-          min_height_squared = height_squared;
-        }
+      auto cstr_bbox = tr.geom_traits().construct_bbox_2_object();
+      auto do_intersect = tr.geom_traits().do_intersect_2_object();
+      auto va = edge.first->vertex(tr.cw (edge.second));
+      auto vb = edge.first->vertex(tr.ccw(edge.second));
+
+      auto closest_encroaching_vertex = get_closest_encroaching_vertex(edge);
+      if(!closest_encroaching_vertex) {
+        return p;
+      }
+      FT edge_length_squared = squared_distance(va->point(), vb->point());
+      if(closest_encroaching_vertex->squared_distance * FT(10000) > edge_length_squared) {
+        return p;
       }
-      if(min_height_squared) {
-        FT edge_length_squared = squared_distance(this->va->point(), this->vb->point());
-        if(*min_height_squared * FT(10000) < edge_length_squared && *min_height_squared < sq_delta) {
-          // the angle is smaller than approx 0.5 degree, and the distance is very small
-          // so we consider that the point is aligned with ab
+
+      // here, the angle is smaller than approx 0.5 degree, and the distance is very small
+      // so we consider that the point is aligned with ab
+
+      auto bbox = cstr_bbox(closest_encroaching_vertex->vertex->point());
+      bbox.dilate(4);
+      if(do_intersect(bbox, typename Tr::Segment(va->point(), vb->point()))
+         )
+      {
 #ifdef CGAL_MESH_2_DEBUG_REFINEMENT_POINTS
-          std::cerr << "return point of existing vertex " << IO::oformat(closest_encroaching_vertex, With_point_tag{})
-                    << std::endl;
-#endif // CGAL_MESH_2_DEBUG_REFINEMENT_POINTS
-          return closest_encroaching_vertex->point();
-        }
+        std::cerr << "return point of existing vertex "
+                  << IO::oformat(closest_encroaching_vertex->vertex, With_point_tag{}) << std::endl;
+  #endif // CGAL_MESH_2_DEBUG_REFINEMENT_POINTS
+        p = closest_encroaching_vertex->vertex->point();
+        return p;
       }
     }
-
     return p;
   }
 

Mesh_2/include/CGAL/Mesh_2/Refine_edges_with_clusters.h

@@ -120,19 +120,21 @@ class Refine_edges_base_with_clusters :
           std::cerr << "midpoint(" << IO::oformat(this->va, With_point_tag{})
                     << " , " << IO::oformat(this->vb, With_point_tag{}) << ")\n";
 #endif // CGAL_MESH_2_DEBUG_CLUSTERS
-          return display_return_point(midpoint(this->va->point(), this->vb->point()));
-        }
-      else {
+          return display_return_point(
+              this->maybe_snap_to_existing_vertex(edge, midpoint(this->va->point(), this->vb->point())));
+      } else {
         // va only is a cluster
         va_has_a_cluster = true;
-        return display_return_point(split_cluster_point(this->va,this->vb,ca));
+        return display_return_point(
+            this->maybe_snap_to_existing_vertex(edge, split_cluster_point(this->va, this->vb, ca)));
       }
     } else
     if( clusters.get_cluster(this->vb,this->va,cb,cb_it) ){
       // vb only is a cluster
       vb_has_a_cluster = true;
-      return display_return_point(split_cluster_point(this->vb,this->va,cb));
-    }else{
+      return display_return_point(
+          this->maybe_snap_to_existing_vertex(edge, split_cluster_point(this->vb, this->va, cb)));
+    } else {
       // no cluster
 
 #ifdef CGAL_MESH_2_DEBUG_CLUSTERS

Mesh_2/test/Mesh_2/triwild-10070-min1.obj

@@ -0,0 +1,8 @@
+v 453.49763280000002 444.98137439999999 0
+v 453.57378188125 445.0501073328125 0
+v 453.77453855000005 445.23131233749996 0
+v 454.05836694375 445.48749872343745 0
+l 1 2
+l 3 4
+l 1 3
+l 4 2