fix approx_is_too_large(edge) that needs the corresponding curve_inde…

…x to be valid

Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h

@@ -273,8 +273,7 @@ class Protect_edges_sizing_field
 
   /// Returns `true` if the edge `(va,vb)` is a not good enough approximation
   /// of its feature.
-  bool approx_is_too_large(const Vertex_handle& va,
-                           const Vertex_handle& vb,
+  bool approx_is_too_large(const Edge& e,
                            const bool is_edge_in_complex) const;
 
   /// Returns `true` if the balls of `va` and `vb` intersect.
@@ -458,10 +457,12 @@ class Protect_edges_sizing_field
       dim = -1 - dim;
 
     const FT s = field(p, dim, index);
-    if(s <= FT(0)) {
+    if(s < minimal_size_)
+    {
       std::stringstream msg;
       msg.precision(17);
-      msg << "Error: the field is null at ";
+      msg << "Error: the field is smaller than minimal size ("
+        << minimal_size_ << ") at ";
       if(dim == 0) msg << "corner (";
       else msg << "point (";
       msg << p << ")";
@@ -645,7 +646,7 @@ insert_corners()
 #endif
 
     // Get weight (the ball radius is given by the 'query_size' function)
-    FT w = CGAL::square(query_size(p, 0, p_index));
+    FT w = (std::min)(minimal_weight_, CGAL::square(query_size(p, 0, p_index)));
 
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
       std::cerr << "Weight from sizing field: " << w << std::endl;
@@ -1428,7 +1429,7 @@ refine_balls()
       if( // topology condition
           non_adjacent_but_intersect(va, vb, is_edge_in_complex)
           // approximation condition
-       || (use_distance_field() && approx_is_too_large(va, vb, is_edge_in_complex)))
+       || (use_distance_field() && approx_is_too_large(*eit, is_edge_in_complex)))
       {
         using CGAL::Mesh_3::internal::distance_divisor;
 
@@ -1557,34 +1558,29 @@ do_balls_intersect(const Vertex_handle& va, const Vertex_handle& vb) const
 template <typename C3T3, typename MD, typename Sf, typename Df>
 bool
 Protect_edges_sizing_field<C3T3, MD, Sf, Df>::
-approx_is_too_large(const Vertex_handle& va, const Vertex_handle& vb, const bool is_edge_in_complex) const
+approx_is_too_large(const Edge& e, const bool is_edge_in_complex) const
 {
   if ( ! is_edge_in_complex )
-  {
     return false;
-  }
-  typedef typename Kernel::Point_3 Point_3;
 
-  Curve_index curve_index = domain_.curve_index((va->in_dimension() < vb->in_dimension()) ? vb->index() : va->index());
+  const auto& [va, vb] = c3t3_.triangulation().vertices(e);
 
-  const Point_3& pa = va->point().point();
-  const Point_3& pb = vb->point().point();
-  const Point_3& segment_middle = CGAL::midpoint(pa, pb);
-  // Obtain the geodesic middle point
-  FT signed_geodesic_distance = domain_.signed_geodesic_distance(pa, pb, curve_index);
-  Point_3 geodesic_middle;
-  if (signed_geodesic_distance >= FT(0))
-  {
-    geodesic_middle = domain_.construct_point_on_curve(pa, curve_index, signed_geodesic_distance / 2);
-  }
-  else
-  {
-    geodesic_middle = domain_.construct_point_on_curve(pb, curve_index, -signed_geodesic_distance / 2);
-  }
-  // Compare distance to the parameter's distance
-  FT squared_evaluated_distance = CGAL::squared_distance(segment_middle, geodesic_middle);
-  FT segment_middle_edge_distance = query_distance(segment_middle, 1, curve_index);
-  return squared_evaluated_distance > CGAL::square(segment_middle_edge_distance);
+  const Bare_point& pa = va->point().point();
+  const Bare_point& pb = vb->point().point();
+
+  // Construct the geodesic middle point
+  const Curve_index curve_index = c3t3_.curve_index(e);
+  const FT signed_geodesic_distance = domain_.signed_geodesic_distance(pa, pb, curve_index);
+  const Bare_point geodesic_middle = (signed_geodesic_distance >= FT(0))
+      ? domain_.construct_point_on_curve(pa, curve_index, signed_geodesic_distance / 2)
+      : domain_.construct_point_on_curve(pb, curve_index, -signed_geodesic_distance / 2);
+
+  const Bare_point edge_middle = CGAL::midpoint(pa, pb);
+  const FT squared_evaluated_distance = CGAL::squared_distance(edge_middle, geodesic_middle);
+
+  // Compare distance to the distance field from criteria
+  const FT max_distance_to_curve = query_distance(edge_middle, 1, curve_index);
+  return squared_evaluated_distance > CGAL::square(max_distance_to_curve);
 }
 
 template <typename C3T3, typename MD, typename Sf, typename Df>