diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h index 4a0e3f131dee..2424d928ad1f 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h +++ b/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 bool Protect_edges_sizing_field:: -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