From da2b0edcd17b25dfc3883485c0eb680bb058d115 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Wed, 27 Aug 2025 16:50:52 +0200 Subject: [PATCH 01/19] cache polyline length and use it in loop case --- .../Mesh_domain_with_polyline_features_3.h | 47 ++++++++++++------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index d31ac841756..d4bc6a3c1e6 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -213,17 +213,18 @@ public: /// returns the length of the polyline FT length() const { - //TODO: cache result - FT result (0); - const_iterator it = points_.begin(); - const_iterator previous = it++; - - for ( const_iterator end = points_.end() ; it != end ; ++it, ++previous ) + if(length_ < 0.) { - result += distance(*previous,*it); - } + FT result(0); + const_iterator it = points_.begin(); + const_iterator previous = it++; - return result; + for(const_iterator end = points_.end(); it != end; ++it, ++previous) { + result += distance(*previous, *it); + } + length_ = result; + } + return length_; } /// returns the signed geodesic distance between `p` and `q`. @@ -244,13 +245,23 @@ public: else { return -result; } } - if(is_loop()) { - const FT positive_distance = curve_segment_length(p, q, CGAL::POSITIVE, pit, qit); - const FT negative_distance = curve_segment_length(p, q, CGAL::NEGATIVE, pit, qit); - return (positive_distance < negative_distance) - ? positive_distance - : (- negative_distance); - } else { + if(is_loop()) + { + FT positive_distance, negative_distance; + if(pit <= qit) + { + positive_distance = curve_segment_length(p, q, CGAL::POSITIVE, pit, qit); + negative_distance = length() - positive_distance; + } + else + { + negative_distance = curve_segment_length(q, p, CGAL::POSITIVE, qit, pit); + positive_distance = length() - negative_distance; + } + return (positive_distance < negative_distance) ? positive_distance : (-negative_distance); + } + else + { return (pit <= qit) ? curve_segment_length(p, q, CGAL::POSITIVE, pit, qit) : ( - curve_segment_length(p, q, CGAL::NEGATIVE, pit, qit) ); @@ -448,6 +459,10 @@ private: public: Data points_; + +private: + mutable FT length_ = -1.; + }; // end class Polyline From c350da1ac74373db67ef13e2f2879666b96612d0 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 28 Aug 2025 12:40:04 +0200 Subject: [PATCH 02/19] make domain.construct_point_on_curve() return location on polyline and use it in domain.curve_segment_is_covered() --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 37 ++++++++++++++++--- .../Mesh_domain_with_polyline_features_3.h | 37 +++++++++++++------ 2 files changed, 57 insertions(+), 17 deletions(-) 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 4cead545328..f3043b95c6b 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 @@ -46,6 +46,7 @@ #include #include #include +#include #include @@ -142,6 +143,8 @@ public: using Distance_Function = typename CGAL::Default::Get::type; + using Polyline_iter = typename CGAL::Mesh_3::internal::Polyline::const_iterator; + private: typedef typename CGAL::Kernel_traits::Kernel Kernel; typedef Delaunay_triangulation_3 Dt; @@ -491,6 +494,17 @@ private: return use_edge_distance_; } + Polyline_iter polyline_locate(Vertex_handle vh) const + { + CGAL_assertion(vh->in_dimension() < 2); + return vertex_to_polyline_iterator_.at(vh); + } + void set_polyline_iterator(Vertex_handle vh, Polyline_iter it) + { + CGAL_assertion(vh->in_dimension() < 2); + vertex_to_polyline_iterator_[vh] = it; + } + private: C3T3& c3t3_; const MeshDomain& domain_; @@ -504,6 +518,7 @@ private: int refine_balls_iteration_nb; bool nonlinear_growth_of_balls; const std::size_t maximal_number_of_vertices_; + std::unordered_map vertex_to_polyline_iterator_; Mesh_error_code* const error_code_; #ifndef CGAL_NO_ATOMIC /// Pointer to the atomic Boolean that can stop the process @@ -1054,7 +1069,7 @@ insert_balls_on_edges() FT curve_length = domain_.curve_length(curve_index); - Bare_point other_point = + auto [other_point, polyline_iter] = domain_.construct_point_on_curve(p, curve_index, curve_length / 2); @@ -1066,7 +1081,9 @@ insert_balls_on_edges() p_index, Vertex_handle(), CGAL::Emptyset_iterator()).first; + set_polyline_iterator(vp, polyline_iter); } + // No 'else' because in that case 'is_vertex(..)' already filled // the variable 'vp'. vq = vp; @@ -1238,10 +1255,11 @@ insert_balls(const Vertex_handle& vp, curve_index, d_sign) << ")\n"; #endif - const Bare_point new_point = + const auto [bp, polyline_iter] = domain_.construct_point_on_curve(cp(vp_wp), curve_index, d_signF * d / 2); + const Bare_point new_point = bp; const int dim = 1; // new_point is on edge const Index index = domain_.index_from_curve_index(curve_index); const FT point_weight = CGAL::square(size_(new_point, dim, index)); @@ -1256,8 +1274,10 @@ insert_balls(const Vertex_handle& vp, index, Vertex_handle(), out); - if(forced_stop()) return out; const Vertex_handle new_vertex = pair.first; + set_polyline_iterator(new_vertex, polyline_iter); + + if(forced_stop()) return out; out = pair.second; const FT sn = get_radius(new_vertex); if(sp <= sn) { @@ -1332,8 +1352,9 @@ insert_balls(const Vertex_handle& vp, for ( int i = 1 ; i <= n ; ++i ) { // New point position - Bare_point new_point = + const auto [bp, polyline_iter] = domain_.construct_point_on_curve(p, curve_index, pt_dist); + const Bare_point& new_point = bp; // Weight (use as size the min between norm_step_size and linear interpolation) FT current_size = (std::min)(norm_step_size, sp + CGAL::abs(pt_dist)/d*(sq-sp)); @@ -1347,6 +1368,7 @@ insert_balls(const Vertex_handle& vp, std::pair pair = smart_insert_point(new_point, point_weight, dim, index, prev, out); Vertex_handle new_vertex = pair.first; + set_polyline_iterator(new_vertex, polyline_iter); out = pair.second; // Add edge to c3t3 @@ -1559,7 +1581,7 @@ approx_is_too_large(const Edge& e, const bool is_edge_in_complex) const // 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)) + const auto [geodesic_middle, _ /*polyline_iter*/] = (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); @@ -1882,7 +1904,10 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2, const bool cov = domain_.is_curve_segment_covered(curve_index, orientation, cp(v1_wp), cp(v2_wp), - cw(v1_wp), cw(v2_wp)); + cw(v1_wp), cw(v2_wp), + polyline_locate(v1), + polyline_locate(v2)); + #if CGAL_MESH_3_PROTECTION_DEBUG & 1 if(cov) { std::cerr << " But the curve is locally covered\n"; diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index d4bc6a3c1e6..a55c6f8058b 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -57,6 +57,7 @@ class Polyline public: typedef typename Data::const_iterator const_iterator; + typedef std::pair Point_and_location; Polyline() {} ~Polyline() {} @@ -119,7 +120,9 @@ public: bool is_curve_segment_covered(CGAL::Orientation orientation, const Point_3& c1, const Point_3& c2, - const FT sq_r1, const FT sq_r2) const + const FT sq_r1, const FT sq_r2, + const const_iterator cc1_it, + const const_iterator cc2_it) const { CGAL_assertion(orientation != CGAL::ZERO); typename Kernel::Has_on_bounded_side_3 cover_pred = @@ -129,8 +132,8 @@ public: const Sphere_3 s1(c1, sq_r1); const Sphere_3 s2(c2, sq_r2); - const_iterator c1_it = locate(c1); - const_iterator c2_it = locate(c2); + const_iterator c1_it = cc1_it; + const_iterator c2_it = cc2_it; if(orientation == CGAL::NEGATIVE) { ++c1_it; @@ -272,7 +275,7 @@ public: /// returns a point at geodesic distance `distance` from p along the /// polyline. The polyline is oriented from starting point to end point. /// The distance could be negative. - Point_3 point_at(const Point_3& p, FT distance) const + Point_and_location point_at(const Point_3& p, FT distance) const { // use first point of the polyline instead of p distance += curve_segment_length(start_point(),p,CGAL::POSITIVE); @@ -302,7 +305,7 @@ public: ++pit; if (pit == points_.end()) - return *previous; + return {*previous, previous}; segment_length = this->distance(*previous,*pit); } @@ -311,7 +314,8 @@ public: typedef typename Kernel::Vector_3 Vector_3; Vector_3 v (*previous, *pit); - return (*previous) + (distance / CGAL::sqrt(v.squared_length())) * v; + return {(*previous) + (distance / CGAL::sqrt(v.squared_length())) * v, + previous}; } bool are_ordered_along(const Point_3& p, const Point_3& q) const @@ -582,6 +586,9 @@ public: typedef GT R; typedef typename MD::Point_3 Point_3; + using Polyline_const_iterator = typename Mesh_3::internal::Polyline::const_iterator; + using Point_and_location = typename Mesh_3::internal::Polyline::Point_and_location; + /// \name Creation /// @{ @@ -732,7 +739,7 @@ public: FT curve_length(const Curve_index& curve_index) const; /// implements `MeshDomainWithFeatures_3::construct_point_on_curve()`. - Point_3 + Point_and_location construct_point_on_curve(const Point_3& starting_point, const Curve_index& curve_index, FT distance) const; @@ -753,7 +760,9 @@ public: bool is_curve_segment_covered(const Curve_index& index, CGAL::Orientation orientation, const Point_3& c1, const Point_3& c2, - const FT sq_r1, const FT sq_r2) const; + const FT sq_r1, const FT sq_r2, + const Polyline_const_iterator c1_it, + const Polyline_const_iterator c2_it) const; /** * Returns the index to be stored in a vertex lying on the surface identified @@ -866,6 +875,8 @@ private: typedef CGAL::AABB_traits_3 AABB_curves_traits; +// typedef typename Polyline::const_iterator Polyline_const_iterator; + Corners corners_; Corners_tmp_incidences corners_tmp_incidences_; Corner_index current_corner_index_; @@ -1030,7 +1041,7 @@ curve_length(const Curve_index& curve_index) const template -typename Mesh_domain_with_polyline_features_3::Point_3 +typename Mesh_domain_with_polyline_features_3::Point_and_location Mesh_domain_with_polyline_features_3:: construct_point_on_curve(const Point_3& starting_point, const Curve_index& curve_index, @@ -1550,13 +1561,17 @@ Mesh_domain_with_polyline_features_3:: is_curve_segment_covered(const Curve_index& index, CGAL::Orientation orientation, const Point_3& c1, const Point_3& c2, - const FT sq_r1, const FT sq_r2) const + const FT sq_r1, const FT sq_r2, + const Polyline_const_iterator c1_it, + const Polyline_const_iterator c2_it) const { typename Edges::const_iterator eit = edges_.find(index); CGAL_assertion(eit != edges_.end()); return eit->second.is_curve_segment_covered(orientation, - c1, c2, sq_r1, sq_r2); + c1, c2, + sq_r1, sq_r2, + c1_it, c2_it); } } //namespace CGAL From 17e3ee7189816d949ce1dbfa9c7566b310501b66 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 29 Aug 2025 15:24:35 +0200 Subject: [PATCH 03/19] use iterator in curve_segment_length() add a special locate_corner() because we don't want to store all incident polylines to a corner, and that it's only a few, cheap, locate operations --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 22 +++++-- .../Mesh_domain_with_polyline_features_3.h | 57 ++++++++++++++++++- 2 files changed, 70 insertions(+), 9 deletions(-) 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 f3043b95c6b..e0a3cade315 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 @@ -494,14 +494,22 @@ private: return use_edge_distance_; } - Polyline_iter polyline_locate(Vertex_handle vh) const + Polyline_iter locate_in_polyline(Vertex_handle vh, const Curve_index& index) const { CGAL_assertion(vh->in_dimension() < 2); - return vertex_to_polyline_iterator_.at(vh); + + if(vh->in_dimension() == 0) //corner + { + auto cp = c3t3_.triangulation().geom_traits().construct_point_3_object(); + return domain_.locate_corner(index, cp(vh->point())); + } + else + return vertex_to_polyline_iterator_.at(vh); } + void set_polyline_iterator(Vertex_handle vh, Polyline_iter it) { - CGAL_assertion(vh->in_dimension() < 2); + CGAL_assertion(vh->in_dimension() == 1); vertex_to_polyline_iterator_[vh] = it; } @@ -1368,8 +1376,8 @@ insert_balls(const Vertex_handle& vp, std::pair pair = smart_insert_point(new_point, point_weight, dim, index, prev, out); Vertex_handle new_vertex = pair.first; - set_polyline_iterator(new_vertex, polyline_iter); out = pair.second; + set_polyline_iterator(new_vertex, polyline_iter); // Add edge to c3t3 if(!c3t3_.is_in_complex(prev, new_vertex)) { @@ -1850,6 +1858,8 @@ curve_segment_length(const Vertex_handle v1, FT arc_length = (v1_valid_curve_index && v2_valid_curve_index) ? domain_.curve_segment_length(cp(v1_wp), cp(v2_wp), + locate_in_polyline(v1, curve_index), + locate_in_polyline(v2, curve_index), curve_index, orientation) : compute_distance(v1, v2); //curve polyline may not be consistent @@ -1905,8 +1915,8 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2, orientation, cp(v1_wp), cp(v2_wp), cw(v1_wp), cw(v2_wp), - polyline_locate(v1), - polyline_locate(v2)); + locate_in_polyline(v1, curve_index), + locate_in_polyline(v2, curve_index)); #if CGAL_MESH_3_PROTECTION_DEBUG & 1 if(cov) { diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index a55c6f8058b..bdc7dee58a4 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -305,7 +305,7 @@ public: ++pit; if (pit == points_.end()) - return {*previous, previous}; + return {*previous, last_segment_source()}; segment_length = this->distance(*previous,*pit); } @@ -318,6 +318,18 @@ public: previous}; } + const_iterator locate_corner(const Point_3& p) const + { + const_iterator res = points_.end(); + if(p == start_point()) + res = points_.begin(); + else if(p == end_point()) + res = last_segment_source(); + + CGAL_assertion(res != points_.end()); + return res; + } + bool are_ordered_along(const Point_3& p, const Point_3& q) const { CGAL_precondition(!is_loop()); @@ -586,8 +598,9 @@ public: typedef GT R; typedef typename MD::Point_3 Point_3; - using Polyline_const_iterator = typename Mesh_3::internal::Polyline::const_iterator; - using Point_and_location = typename Mesh_3::internal::Polyline::Point_and_location; + using Polyline = Mesh_3::internal::Polyline; + using Polyline_const_iterator = typename Polyline::const_iterator; + using Point_and_location = typename Polyline::Point_and_location; /// \name Creation /// @{ @@ -735,6 +748,13 @@ public: const Curve_index& curve_index, CGAL::Orientation orientation) const; + FT curve_segment_length(const Point_3& p, + const Point_3 q, + const Polyline_const_iterator p_it, + const Polyline_const_iterator q_it, + const Curve_index& curve_index, + CGAL::Orientation orientation) const; + /// implements `MeshDomainWithFeatures_3::curve_length()`. FT curve_length(const Curve_index& curve_index) const; @@ -764,6 +784,10 @@ public: const Polyline_const_iterator c1_it, const Polyline_const_iterator c2_it) const; + /// locates the corner point `p` on the curve identified by `curve_index` + Polyline_const_iterator locate_corner(const Curve_index& curve_index, + const Point_3& p) const; + /** * Returns the index to be stored in a vertex lying on the surface identified * by `index`. @@ -1026,6 +1050,22 @@ curve_segment_length(const Point_3& p, const Point_3 q, return eit->second.curve_segment_length(p, q, orientation); } +template +typename Mesh_domain_with_polyline_features_3::FT +Mesh_domain_with_polyline_features_3:: +curve_segment_length(const Point_3& p, + const Point_3 q, + const Polyline_const_iterator p_it, + const Polyline_const_iterator q_it, + const Curve_index& curve_index, + CGAL::Orientation orientation) const +{ + // Get corresponding polyline + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + + return eit->second.curve_segment_length(p, q, orientation, p_it, q_it); +} template typename Mesh_domain_with_polyline_features_3::FT @@ -1574,6 +1614,17 @@ is_curve_segment_covered(const Curve_index& index, c1_it, c2_it); } +template +typename Mesh_domain_with_polyline_features_3::Polyline_const_iterator +Mesh_domain_with_polyline_features_3:: +locate_corner(const Curve_index& curve_index, + const Point_3& p) const +{ + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + return eit->second.locate_corner(p); +} + } //namespace CGAL #endif // CGAL_MESH_DOMAIN_WITH_POLYLINE_FEATURES_3_H From dbcbff15fcd87f9eaa6457e07f8eeea560751946 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 29 Aug 2025 16:40:07 +0200 Subject: [PATCH 04/19] use iterator for construct_point_on_curve() --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 34 +++++++++++++------ .../Mesh_domain_with_polyline_features_3.h | 15 +++++--- 2 files changed, 33 insertions(+), 16 deletions(-) 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 e0a3cade315..a1807c4dd99 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 @@ -1032,6 +1032,7 @@ insert_balls_on_edges() struct Feature_tuple { Curve_index curve_index_; + Polyline_iter polyline_begin_; std::pair point_s_; std::pair point_t_; }; @@ -1049,16 +1050,16 @@ insert_balls_on_edges() std::cerr << "\n** treat curve #" << curve_index << std::endl; #endif const Bare_point& p = ft.point_s_.first; - const Bare_point& q = ft.point_t_.first; - const Index& p_index = ft.point_s_.second; - const Index& q_index = ft.point_t_.second; Vertex_handle vp,vq; if ( ! domain_.is_loop(curve_index) ) { vp = get_vertex_corner_from_point(p,p_index); - vq = get_vertex_corner_from_point(q,q_index); + + const Bare_point& q = ft.point_t_.first; + const Index& q_index = ft.point_t_.second; + vq = get_vertex_corner_from_point(q, q_index); } else { @@ -1080,7 +1081,8 @@ insert_balls_on_edges() auto [other_point, polyline_iter] = domain_.construct_point_on_curve(p, curve_index, - curve_length / 2); + curve_length / 2, + ft.polyline_begin_); p_size = (std::min)(p_size, compute_distance(p, other_point) / 3); vp = smart_insert_point(p, @@ -1266,7 +1268,8 @@ insert_balls(const Vertex_handle& vp, const auto [bp, polyline_iter] = domain_.construct_point_on_curve(cp(vp_wp), curve_index, - d_signF * d / 2); + d_signF * d / 2, + locate_in_polyline(vp, curve_index)); const Bare_point new_point = bp; const int dim = 1; // new_point is on edge const Index index = domain_.index_from_curve_index(curve_index); @@ -1361,7 +1364,7 @@ insert_balls(const Vertex_handle& vp, { // New point position const auto [bp, polyline_iter] = - domain_.construct_point_on_curve(p, curve_index, pt_dist); + domain_.construct_point_on_curve(p, curve_index, pt_dist, locate_in_polyline(vp, curve_index)); const Bare_point& new_point = bp; // Weight (use as size the min between norm_step_size and linear interpolation) @@ -1583,15 +1586,24 @@ approx_is_too_large(const Edge& e, const bool is_edge_in_complex) const if ( ! is_edge_in_complex ) return false; - const Bare_point& pa = e.first->vertex(e.second)->point().point(); - const Bare_point& pb = e.first->vertex(e.third)->point().point(); + Vertex_handle va = e.first->vertex(e.second); + Vertex_handle vb = e.first->vertex(e.third); + + 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 auto [geodesic_middle, _ /*polyline_iter*/] = (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); + ? domain_.construct_point_on_curve(pa, + curve_index, + signed_geodesic_distance / 2, + locate_in_polyline(va, curve_index)) + : domain_.construct_point_on_curve(pb, + curve_index, + -signed_geodesic_distance / 2, + locate_in_polyline(vb, curve_index)); const Bare_point edge_middle = CGAL::midpoint(pa, pb); const FT squared_evaluated_distance = CGAL::squared_distance(edge_middle, geodesic_middle); diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index bdc7dee58a4..83ba77ec009 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -275,10 +275,12 @@ public: /// returns a point at geodesic distance `distance` from p along the /// polyline. The polyline is oriented from starting point to end point. /// The distance could be negative. - Point_and_location point_at(const Point_3& p, FT distance) const + Point_and_location point_at(const Point_3& p, + FT distance, + const_iterator p_it) const { // use first point of the polyline instead of p - distance += curve_segment_length(start_point(),p,CGAL::POSITIVE); + distance += curve_segment_length(start_point(),p,CGAL::POSITIVE, points_.begin(), p_it); // If polyline is a loop, ensure that distance is given from start_point() if ( is_loop() ) @@ -762,7 +764,8 @@ public: Point_and_location construct_point_on_curve(const Point_3& starting_point, const Curve_index& curve_index, - FT distance) const; + FT distance, + Polyline_const_iterator starting_point_it) const; /// implements `MeshDomainWithFeatures_3::distance_sign_along_loop()`. CGAL::Sign distance_sign_along_loop(const Point_3& p, const Point_3& q, @@ -1012,6 +1015,7 @@ get_curves(OutputIterator out) const } *out++ = {eit->first, + eit->second.points_.begin(), std::make_pair(p,p_index), std::make_pair(q,q_index)}; } @@ -1085,14 +1089,15 @@ typename Mesh_domain_with_polyline_features_3::Point_and_location Mesh_domain_with_polyline_features_3:: construct_point_on_curve(const Point_3& starting_point, const Curve_index& curve_index, - FT distance) const + FT distance, + Polyline_const_iterator starting_point_it) const { // Get corresponding polyline typename Edges::const_iterator eit = edges_.find(curve_index); CGAL_assertion(eit != edges_.end()); // Return point at geodesic_distance distance from starting_point - return eit->second.point_at(starting_point,distance); + return eit->second.point_at(starting_point, distance, starting_point_it); } From 46d0d0799ef7e07a103dcc0a628269d12ab805ae Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 2 Sep 2025 09:35:52 +0200 Subject: [PATCH 05/19] use saved polyline_iterator everywhere mesh_polyhedral_domain_with_features works fine with this version --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 26 ++++-- .../Mesh_domain_with_polyline_features_3.h | 93 +++++++++++++++---- 2 files changed, 90 insertions(+), 29 deletions(-) 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 a1807c4dd99..bf292c68236 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 @@ -143,7 +143,7 @@ public: using Distance_Function = typename CGAL::Default::Get::type; - using Polyline_iter = typename CGAL::Mesh_3::internal::Polyline::const_iterator; + using Polyline_iterator = typename CGAL::Mesh_3::internal::Polyline::const_iterator; private: typedef typename CGAL::Kernel_traits::Kernel Kernel; @@ -494,7 +494,7 @@ private: return use_edge_distance_; } - Polyline_iter locate_in_polyline(Vertex_handle vh, const Curve_index& index) const + Polyline_iterator locate_in_polyline(Vertex_handle vh, const Curve_index& index) const { CGAL_assertion(vh->in_dimension() < 2); @@ -507,7 +507,7 @@ private: return vertex_to_polyline_iterator_.at(vh); } - void set_polyline_iterator(Vertex_handle vh, Polyline_iter it) + void set_polyline_iterator(Vertex_handle vh, Polyline_iterator it) { CGAL_assertion(vh->in_dimension() == 1); vertex_to_polyline_iterator_[vh] = it; @@ -526,7 +526,7 @@ private: int refine_balls_iteration_nb; bool nonlinear_growth_of_balls; const std::size_t maximal_number_of_vertices_; - std::unordered_map vertex_to_polyline_iterator_; + std::unordered_map vertex_to_polyline_iterator_; Mesh_error_code* const error_code_; #ifndef CGAL_NO_ATOMIC /// Pointer to the atomic Boolean that can stop the process @@ -1032,7 +1032,7 @@ insert_balls_on_edges() struct Feature_tuple { Curve_index curve_index_; - Polyline_iter polyline_begin_; + Polyline_iterator polyline_begin_; std::pair point_s_; std::pair point_t_; }; @@ -1592,18 +1592,22 @@ approx_is_too_large(const Edge& e, const bool is_edge_in_complex) const 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); + Polyline_iterator pa_it = locate_in_polyline(va, curve_index); + Polyline_iterator pb_it = locate_in_polyline(vb, curve_index); + + // Construct the geodesic middle point + const FT signed_geodesic_distance + = domain_.signed_geodesic_distance(pa, pb, pa_it, pb_it, curve_index); const auto [geodesic_middle, _ /*polyline_iter*/] = (signed_geodesic_distance >= FT(0)) ? domain_.construct_point_on_curve(pa, curve_index, signed_geodesic_distance / 2, - locate_in_polyline(va, curve_index)) + pa_it) : domain_.construct_point_on_curve(pb, curve_index, -signed_geodesic_distance / 2, - locate_in_polyline(vb, curve_index)); + pb_it); const Bare_point edge_middle = CGAL::midpoint(pa, pb); const FT squared_evaluated_distance = CGAL::squared_distance(edge_middle, geodesic_middle); @@ -1973,7 +1977,9 @@ orientation_of_walk(const Vertex_handle& start, cp(start_wp), cp(next_wp), cp(next_along_curve_wp), curve_index); } else { // otherwise, the sign is just the sign of the geodesic distance - return domain_.distance_sign(cp(start_wp), cp(next_wp), curve_index); + return domain_.distance_sign(cp(start_wp), cp(next_wp), curve_index, + locate_in_polyline(start, curve_index), + locate_in_polyline(next, curve_index)); } } diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 83ba77ec009..61974b4468a 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -237,6 +237,12 @@ public: const_iterator pit = locate(p); const_iterator qit = locate(q,false); + return signed_geodesic_distance(p, q, pit, qit); + } + + FT signed_geodesic_distance(const Point_3& p, const Point_3& q, + const_iterator pit, const_iterator qit) const + { // If p and q are in the same segment of the polyline if ( pit == qit ) { @@ -271,29 +277,63 @@ public: } } + const_iterator previous_segment_source(const_iterator it) const + { + CGAL_assertion(it != points_.end()); + if(it == first_segment_source()) + { + CGAL_assertion(is_loop()); + it = last_segment_source(); + } + else + { + --it; + } + return it; + } - /// returns a point at geodesic distance `distance` from p along the + /// returns a point at geodesic distance `distance` from `start_pt` along the /// polyline. The polyline is oriented from starting point to end point. /// The distance could be negative. - Point_and_location point_at(const Point_3& p, + Point_and_location point_at(const Point_3& start_pt, FT distance, - const_iterator p_it) const + const_iterator start_it) const { - // use first point of the polyline instead of p - distance += curve_segment_length(start_point(),p,CGAL::POSITIVE, points_.begin(), p_it); - - // If polyline is a loop, ensure that distance is given from start_point() - if ( is_loop() ) + Point_3 p = start_pt; + if(start_it == first_segment_source()) { - if ( distance < FT(0) ) { distance += length(); } - else if ( distance > length() ) { distance -= length(); } + if(is_loop() && distance < 0.) + { + p = *start_it; + start_it = previous_segment_source(start_it); + distance += last_segment_length(); + } } + while(distance < 0.) + { + // Move to previous point + distance += this->distance(p, *start_it); + p = *start_it; + if(start_it != first_segment_source()) + start_it = previous_segment_source(start_it); + } + + // use first point of the polyline instead of p + distance += this->distance(*start_it, p); + +// // If polyline is a loop, ensure that distance is given from start_point() +// if ( is_loop() ) +// { +// if ( distance < FT(0) ) { distance += length(); } +// else if ( distance > length() ) { distance -= length(); } +// } + CGAL_assertion( distance >= FT(0) ); CGAL_assertion( distance <= length() ); // Initialize iterators - const_iterator pit = points_.begin(); + const_iterator pit = start_it; const_iterator previous = pit++; // Iterate to find which segment contains the point we want to construct @@ -332,13 +372,14 @@ public: return res; } - bool are_ordered_along(const Point_3& p, const Point_3& q) const + bool are_ordered_along(const Point_3& p, const Point_3& q, + const_iterator pit, const_iterator qit) const { CGAL_precondition(!is_loop()); - // Locate p & q on polyline - const_iterator pit = locate(p); - const_iterator qit = locate(q,true); +// // Locate p & q on polyline +// const_iterator pit = locate(p); +// const_iterator qit = locate(q,true); // Points are not located on the same segment if ( pit != qit ) { return (pit <= qit); } @@ -360,6 +401,12 @@ private: return (points_.end() - 2); } + FT last_segment_length() const + { + CGAL_precondition(is_valid()); + return distance( *(points_.end() - 2), *(points_.end() - 1) ); + } + /// returns an iterator on the starting point of the segment of the /// polyline which contains p /// if end_point_first is true, then --end is returned instead of begin @@ -774,7 +821,9 @@ public: /// implements `MeshDomainWithFeatures_3::distance_sign()`. CGAL::Sign distance_sign(const Point_3& p, const Point_3& q, - const Curve_index& index) const; + const Curve_index& index, + Polyline_const_iterator pit, + Polyline_const_iterator qit) const; /// implements `MeshDomainWithFeatures_3::is_loop()`. bool is_loop(const Curve_index& index) const; @@ -844,6 +893,8 @@ public: #endif // CGAL_NO_DEPRECATED_CODE FT signed_geodesic_distance(const Point_3& p, const Point_3& q, + Polyline_const_iterator pit, + Polyline_const_iterator qit, const Curve_index& curve_index) const; template @@ -1255,6 +1306,8 @@ template typename Mesh_domain_with_polyline_features_3::FT Mesh_domain_with_polyline_features_3:: signed_geodesic_distance(const Point_3& p, const Point_3& q, + Polyline_const_iterator pit, + Polyline_const_iterator qit, const Curve_index& curve_index) const { // Get corresponding polyline @@ -1262,7 +1315,7 @@ signed_geodesic_distance(const Point_3& p, const Point_3& q, CGAL_assertion(eit != edges_.end()); // Compute geodesic_distance - return eit->second.signed_geodesic_distance(p,q); + return eit->second.signed_geodesic_distance(p, q, pit, qit); } @@ -1548,7 +1601,9 @@ template CGAL::Sign Mesh_domain_with_polyline_features_3:: distance_sign(const Point_3& p, const Point_3& q, - const Curve_index& index) const + const Curve_index& index, + Polyline_const_iterator pit, + Polyline_const_iterator qit) const { typename Edges::const_iterator eit = edges_.find(index); CGAL_assertion(eit != edges_.end()); @@ -1556,7 +1611,7 @@ distance_sign(const Point_3& p, const Point_3& q, if ( p == q ) return CGAL::ZERO; - else if ( eit->second.are_ordered_along(p,q) ) + else if ( eit->second.are_ordered_along(p,q,pit,qit) ) return CGAL::POSITIVE; else return CGAL::NEGATIVE; From 389bb1139085bb239e07ed76373dfc914d72dea6 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 23 Sep 2025 13:06:59 +0200 Subject: [PATCH 06/19] fix point_at() for the loop case --- .../Mesh_domain_with_polyline_features_3.h | 104 +++++++++++------- 1 file changed, 64 insertions(+), 40 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 61974b4468a..5acc1002529 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -292,6 +292,23 @@ public: return it; } + const_iterator next_segment_source(const_iterator it) const + { + CGAL_assertion(it != points_.end()); + if(it == last_segment_source()) + { + if(is_loop()) + return first_segment_source(); + else + return points_.end(); + } + else + { + ++it; + return it; + } + } + /// returns a point at geodesic distance `distance` from `start_pt` along the /// polyline. The polyline is oriented from starting point to end point. /// The distance could be negative. @@ -299,65 +316,53 @@ public: FT distance, const_iterator start_it) const { - Point_3 p = start_pt; - if(start_it == first_segment_source()) + distance += curve_segment_length(*start_it, start_pt, CGAL::POSITIVE); + + // If polyline is a loop, ensure that distance is given from start_it + if(is_loop()) { - if(is_loop() && distance < 0.) - { - p = *start_it; - start_it = previous_segment_source(start_it); - distance += last_segment_length(); - } + std::cout << "is loop" << std::endl; + if(distance < FT(0)) { distance += length(); } + else if(distance > length()) { distance -= length(); } } - while(distance < 0.) - { - // Move to previous point - distance += this->distance(p, *start_it); - p = *start_it; - if(start_it != first_segment_source()) - start_it = previous_segment_source(start_it); - } + CGAL_assertion(distance >= FT(0)); + CGAL_assertion(distance <= length()); - // use first point of the polyline instead of p - distance += this->distance(*start_it, p); - -// // If polyline is a loop, ensure that distance is given from start_point() -// if ( is_loop() ) -// { -// if ( distance < FT(0) ) { distance += length(); } -// else if ( distance > length() ) { distance -= length(); } -// } - - CGAL_assertion( distance >= FT(0) ); - CGAL_assertion( distance <= length() ); + if(distance == FT(0) || distance == length()) + return {*start_it, start_it}; // Initialize iterators - const_iterator pit = start_it; + const_iterator pit = start_it; // start at start_it, and walk forward const_iterator previous = pit++; // Iterate to find which segment contains the point we want to construct - FT segment_length = this->distance(*previous,*pit); - while ( distance > segment_length ) + FT segment_length = this->distance(*previous, *pit); + while(distance > segment_length) { distance -= segment_length; // Increment iterators and update length - ++previous; - ++pit; + previous = next_segment_source(previous); + pit = next_segment_source(pit); - if (pit == points_.end()) - return {*previous, last_segment_source()}; - - segment_length = this->distance(*previous,*pit); + if(pit == points_.end()) + { + CGAL_assertion(distance < this->distance(*previous, end_point())); + break; + } + else + segment_length = this->distance(*previous, *pit); } // return point at distance from current segment source - typedef typename Kernel::Vector_3 Vector_3; - Vector_3 v (*previous, *pit); + using Vector_3 = typename Kernel::Vector_3; + auto vector = Kernel().construct_vector_3_object(); + Vector_3 v = (pit != points_.end()) ? vector(*previous, *pit) + : vector(*previous, end_point()); return {(*previous) + (distance / CGAL::sqrt(v.squared_length())) * v, - previous}; + previous}; } const_iterator locate_corner(const Point_3& p) const @@ -1018,6 +1023,25 @@ public: #endif } // build_curves_aabb_tree() + void dump_curve(const Curve_index& index, const std::string& prefix) const + { + std::string filename(prefix); + filename += std::to_string(index) + ".polylines.txt"; + + std::ofstream os(filename); + typename Edges::const_iterator eit = edges_.find(index); + if(eit == edges_.end()) { + os << "No curve with index " << index << std::endl; + return; + } + const Polyline& polyline = eit->second; + os << polyline.points_.size(); + for(const auto& p : polyline.points_) + os << " " << p; + os << std::endl; + os.close(); + } + /// @endcond }; // class Mesh_domain_with_polyline_features_3 From 2e15f13dc45a77d99e0f0a4f0c6efdaa43a9351e Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 26 Sep 2025 10:16:16 +0200 Subject: [PATCH 07/19] remove cout --- Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 5acc1002529..938ec3d14fe 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -321,7 +321,6 @@ public: // If polyline is a loop, ensure that distance is given from start_it if(is_loop()) { - std::cout << "is loop" << std::endl; if(distance < FT(0)) { distance += length(); } else if(distance > length()) { distance -= length(); } } From 2d2773e4ed5d549c4a86b46372461806f16608e7 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 26 Sep 2025 10:17:56 +0200 Subject: [PATCH 08/19] fix set_polyline_iterator() in insert_balls_on_edges() the iterator we need is the one of p, not the one of `other_point`! + constify what can be + make loops more readable --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 58 ++++++++++--------- 1 file changed, 32 insertions(+), 26 deletions(-) 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 bf292c68236..a58b1b52767 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 @@ -621,7 +621,9 @@ Protect_edges_sizing_field:: insert_corners() { // Iterate on domain corners - typedef std::vector< std::pair > Initial_corners; + using Initial_corner = std::pair; + using Initial_corners = std::vector; + Initial_corners corners; domain_.get_corners(std::back_inserter(corners)); @@ -632,20 +634,16 @@ insert_corners() #endif Dt dt; - for ( typename Initial_corners::iterator it = corners.begin(), - end = corners.end() ; it != end ; ++it ) + for ( const auto& [_,p] : corners ) { if(forced_stop()) break; - const Bare_point& p = it->second; dt.insert(p); } - for ( typename Initial_corners::iterator cit = corners.begin(), - end = corners.end() ; cit != end ; ++cit ) + for ( const auto& [corner_index, p] : corners ) { if(forced_stop()) break; - const Bare_point& p = cit->second; - Index p_index = domain_.index_from_corner_index(cit->first); + Index p_index = domain_.index_from_corner_index(corner_index); #if CGAL_MESH_3_PROTECTION_DEBUG & 1 std::cerr << "\n** treat corner #" << CGAL::IO::oformat(p_index) << std::endl; @@ -694,7 +692,7 @@ insert_corners() // As C3t3::add_to_complex modifies the 'in_dimension' of the vertex, // we need to backup and re-set the 'is_special' marker after. const bool special_ball = is_special(v); - c3t3_.add_to_complex(v,cit->first); + c3t3_.add_to_complex(v, corner_index); if(special_ball) { set_special(v); } @@ -1051,6 +1049,7 @@ insert_balls_on_edges() #endif const Bare_point& p = ft.point_s_.first; const Index& p_index = ft.point_s_.second; + const Polyline_iterator& p_polyline_iter = ft.polyline_begin_; Vertex_handle vp,vq; if ( ! domain_.is_loop(curve_index) ) @@ -1078,7 +1077,7 @@ insert_balls_on_edges() FT curve_length = domain_.curve_length(curve_index); - auto [other_point, polyline_iter] = + auto [other_point, _] = domain_.construct_point_on_curve(p, curve_index, curve_length / 2, @@ -1091,7 +1090,7 @@ insert_balls_on_edges() p_index, Vertex_handle(), CGAL::Emptyset_iterator()).first; - set_polyline_iterator(vp, polyline_iter); + set_polyline_iterator(vp, p_polyline_iter); } // No 'else' because in that case 'is_vertex(..)' already filled @@ -1305,7 +1304,7 @@ insert_balls(const Vertex_handle& vp, } } // nonlinear_growth_of_balls - FT r = (sq - sp) / FT(n+1); + const FT r = (sq - sp) / FT(n+1); #if CGAL_MESH_3_PROTECTION_DEBUG & 1 std::cerr << " n=" << n @@ -1314,9 +1313,9 @@ insert_balls(const Vertex_handle& vp, // Adjust size of steps, D = covered distance - FT D = sp*FT(n+1) + FT((n+1)*(n+2)) / FT(2) * r ; + const FT D = sp*FT(n+1) + FT((n+1)*(n+2)) / FT(2) * r ; - FT dleft_frac = d / D; + const FT dleft_frac = d / D; // Initialize step sizes FT step_size = sp + r; @@ -1359,21 +1358,24 @@ insert_balls(const Vertex_handle& vp, #endif } + // Index and dimension + const int dim = 1; // new_point is on edge + const Index index = domain_.index_from_curve_index(curve_index); + // Launch balls - for ( int i = 1 ; i <= n ; ++i ) + Polyline_iterator p_loc = locate_in_polyline(vp, curve_index); + Bare_point prev_pt = p; + FT dist_to_prev = pt_dist; + + for(int i = 1; i <= n; ++i) { // New point position - const auto [bp, polyline_iter] = - domain_.construct_point_on_curve(p, curve_index, pt_dist, locate_in_polyline(vp, curve_index)); - const Bare_point& new_point = bp; + const auto [new_point, polyline_iter] = + domain_.construct_point_on_curve(prev_pt, curve_index, dist_to_prev, p_loc); // Weight (use as size the min between norm_step_size and linear interpolation) - FT current_size = (std::min)(norm_step_size, sp + CGAL::abs(pt_dist)/d*(sq-sp)); - FT point_weight = current_size * current_size; - - // Index and dimension - Index index = domain_.index_from_curve_index(curve_index); - int dim = 1; // new_point is on edge + const FT current_size = (std::min)(norm_step_size, sp + CGAL::abs(pt_dist)/d*(sq-sp)); + const FT point_weight = current_size * current_size; // Insert point into c3t3 std::pair pair = @@ -1392,8 +1394,12 @@ insert_balls(const Vertex_handle& vp, step_size += r; norm_step_size = dleft_frac * step_size; - // Increment distance - pt_dist += d_signF * norm_step_size; + // Update distance + dist_to_prev = d_signF* norm_step_size; + pt_dist += dist_to_prev; + + prev_pt = new_point; + p_loc = polyline_iter; } // Insert last edge into c3t3 From 26130f161d3004b0d14f3b13f63cbcac6d7eea00 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 26 Sep 2025 15:19:18 +0200 Subject: [PATCH 09/19] fix API for Sizing_field_with_aabb_tree, and deal with the case where distance < 0 in point_at(), because we don't want to start from start_point all over --- .../Mesh_domain_with_polyline_features_3.h | 65 +++++++++++++++++-- 1 file changed, 60 insertions(+), 5 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 938ec3d14fe..baa51e964d6 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -309,6 +309,14 @@ public: } } + /// returns a point at geodesic distance `distance` from p along the + /// polyline. The polyline is oriented from starting point to end point. + /// The distance could be negative. + Point_3 point_at(const Point_3& start_pt, FT distance) const + { + return point_at(start_pt, distance, locate(start_pt)).first; + } + /// returns a point at geodesic distance `distance` from `start_pt` along the /// polyline. The polyline is oriented from starting point to end point. /// The distance could be negative. @@ -316,7 +324,8 @@ public: FT distance, const_iterator start_it) const { - distance += curve_segment_length(*start_it, start_pt, CGAL::POSITIVE); + distance += curve_segment_length(*start_it, start_pt, CGAL::POSITIVE, + start_it, start_it); // If polyline is a loop, ensure that distance is given from start_it if(is_loop()) @@ -324,13 +333,19 @@ public: if(distance < FT(0)) { distance += length(); } else if(distance > length()) { distance -= length(); } } + else if(distance < FT(0)) // If polyline is not a loop and distance is < 0, go backward + { + Point_3 new_start = start_pt; + while(distance < FT(0)) + { + start_it = previous_segment_source(start_it); + distance += this->distance(new_start, *start_it); + } + } CGAL_assertion(distance >= FT(0)); CGAL_assertion(distance <= length()); - if(distance == FT(0) || distance == length()) - return {*start_it, start_it}; - // Initialize iterators const_iterator pit = start_it; // start at start_it, and walk forward const_iterator previous = pit++; @@ -348,7 +363,7 @@ public: if(pit == points_.end()) { CGAL_assertion(distance < this->distance(*previous, end_point())); - break; + break; // return {*previous, previous} } else segment_length = this->distance(*previous, *pit); @@ -817,6 +832,13 @@ public: const Curve_index& curve_index, FT distance, Polyline_const_iterator starting_point_it) const; + + /// implements `MeshDomainWithFeatures_3::construct_point_on_curve()`. + Point_3 + construct_point_on_curve(const Point_3& starting_point, + const Curve_index& curve_index, + FT distance) const; + /// implements `MeshDomainWithFeatures_3::distance_sign_along_loop()`. CGAL::Sign distance_sign_along_loop(const Point_3& p, const Point_3& q, @@ -896,6 +918,9 @@ public: } #endif // CGAL_NO_DEPRECATED_CODE + FT signed_geodesic_distance(const Point_3& p, const Point_3& q, + const Curve_index& curve_index) const; + FT signed_geodesic_distance(const Point_3& p, const Point_3& q, Polyline_const_iterator pit, Polyline_const_iterator qit, @@ -1174,6 +1199,21 @@ construct_point_on_curve(const Point_3& starting_point, return eit->second.point_at(starting_point, distance, starting_point_it); } +template +typename Mesh_domain_with_polyline_features_3::Point_3 +Mesh_domain_with_polyline_features_3:: +construct_point_on_curve(const Point_3& starting_point, + const Curve_index& curve_index, + FT distance) const +{ + // Get corresponding polyline + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + + // Return point at geodesic_distance distance from starting_point + return eit->second.point_at(starting_point, distance); +} + /// @cond DEVELOPERS template @@ -1325,6 +1365,21 @@ add_features_and_incidences(InputIterator first, InputIterator end, } /// @cond DEVELOPERS +template +typename Mesh_domain_with_polyline_features_3::FT +Mesh_domain_with_polyline_features_3:: +signed_geodesic_distance(const Point_3& p, + const Point_3& q, + const Curve_index& curve_index) const +{ + // Get corresponding polyline + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + + // Compute geodesic_distance + return eit->second.signed_geodesic_distance(p, q); +} + template typename Mesh_domain_with_polyline_features_3::FT Mesh_domain_with_polyline_features_3:: From 2c98528a97a0f7eb385e8e49e714bc10e045f44b Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Mon, 29 Sep 2025 11:32:53 +0200 Subject: [PATCH 10/19] avoid redefinition of Polyline --- Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index baa51e964d6..9544313d071 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -968,8 +968,6 @@ private: private: typedef std::map Corners; - - typedef Mesh_3::internal::Polyline Polyline; typedef std::map Edges; typedef std::map Edges_incidences; typedef std::map > Corners_tmp_incidences; From 2bde2295f4f9188e40e00516efc2e08f62a5a031 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 3 Oct 2025 13:55:13 +0200 Subject: [PATCH 11/19] fix location of start_point in construct_point_on_curve() and add tons of assertions --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 1 + .../Mesh_domain_with_polyline_features_3.h | 47 ++++++++++++++----- 2 files changed, 37 insertions(+), 11 deletions(-) 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 a58b1b52767..d76d8c5854f 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 @@ -1364,6 +1364,7 @@ insert_balls(const Vertex_handle& vp, // Launch balls Polyline_iterator p_loc = locate_in_polyline(vp, curve_index); + CGAL_assertion(p_loc == domain_.locate_point(curve_index, p)); Bare_point prev_pt = p; FT dist_to_prev = pt_dist; diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 9544313d071..bf14dcb9b51 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -171,6 +171,8 @@ public: const_iterator q_it) const { CGAL_assertion(orientation != CGAL::ZERO); + CGAL_assertion(p_it == locate(p)); + CGAL_assertion(q_it == locate(q)); if(p_it == q_it) { const CGAL::Comparison_result cmp = compare_distance(*p_it,p,q); @@ -243,6 +245,9 @@ public: FT signed_geodesic_distance(const Point_3& p, const Point_3& q, const_iterator pit, const_iterator qit) const { + CGAL_assertion(pit == locate(p)); + CGAL_assertion(qit == locate(q, false)); + // If p and q are in the same segment of the polyline if ( pit == qit ) { @@ -324,8 +329,14 @@ public: FT distance, const_iterator start_it) const { - distance += curve_segment_length(*start_it, start_pt, CGAL::POSITIVE, - start_it, start_it); + CGAL_assertion(start_it == locate(start_pt)); + + const Point_3& start_it_pt = *start_it; + const_iterator start_it_locate_pt + = (start_it == points_.begin()) ? start_it : std::prev(start_it); + + distance += curve_segment_length(start_it_pt, start_pt, CGAL::POSITIVE, + start_it_locate_pt, start_it); // If polyline is a loop, ensure that distance is given from start_it if(is_loop()) @@ -340,6 +351,7 @@ public: { start_it = previous_segment_source(start_it); distance += this->distance(new_start, *start_it); + new_start = *start_it; } } @@ -387,18 +399,24 @@ public: else if(p == end_point()) res = last_segment_source(); + CGAL_assertion(res == locate(p)); CGAL_assertion(res != points_.end()); return res; } + const_iterator locate_point(const Point_3& p) const + { + return locate(p); + } + bool are_ordered_along(const Point_3& p, const Point_3& q, const_iterator pit, const_iterator qit) const { CGAL_precondition(!is_loop()); -// // Locate p & q on polyline -// const_iterator pit = locate(p); -// const_iterator qit = locate(q,true); + // Locate p & q on polyline + CGAL_assertion(pit == locate(p)); + CGAL_assertion(qit == locate(q, true)); // Points are not located on the same segment if ( pit != qit ) { return (pit <= qit); } @@ -420,12 +438,6 @@ private: return (points_.end() - 2); } - FT last_segment_length() const - { - CGAL_precondition(is_valid()); - return distance( *(points_.end() - 2), *(points_.end() - 1) ); - } - /// returns an iterator on the starting point of the segment of the /// polyline which contains p /// if end_point_first is true, then --end is returned instead of begin @@ -866,6 +878,8 @@ public: Polyline_const_iterator locate_corner(const Curve_index& curve_index, const Point_3& p) const; + Polyline_const_iterator locate_point(const Curve_index& curve_index, const Point_3& p) const; + /** * Returns the index to be stored in a vertex lying on the surface identified * by `index`. @@ -1761,6 +1775,17 @@ locate_corner(const Curve_index& curve_index, return eit->second.locate_corner(p); } +template +typename Mesh_domain_with_polyline_features_3::Polyline_const_iterator +Mesh_domain_with_polyline_features_3:: +locate_point(const Curve_index& curve_index, + const Point_3& p) const +{ + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + return eit->second.locate_point(p); +} + } //namespace CGAL #endif // CGAL_MESH_DOMAIN_WITH_POLYLINE_FEATURES_3_H From 0a6df0877a5c5f9069b8985fbd98a003eae974ab Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 13 Nov 2025 15:38:20 +0100 Subject: [PATCH 12/19] wip using locate_corner() --- Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h index 74346f5ce94..2471f6028b1 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -228,7 +228,8 @@ public: //fill incidences of corners with curves d_ptr->corners_incident_curves.resize(d_ptr->corners.size()); - for(const typename Corners_indices::value_type& pair : d_ptr->corners_indices) { + for(const typename Corners_indices::value_type& pair : d_ptr->corners_indices) + { d_ptr->dt.insert(pair.first); // Fill `corners_incident_curves[corner_id]` @@ -236,15 +237,19 @@ public: d_ptr->domain.get_corner_incident_curves(pair.second, std::inserter(incident_curves, incident_curves.end())); - // For each incident loops, insert a point on the loop, as far as + // For each incident loop, insert a point on the loop, as far as // possible. - for(Curve_index curve_index : incident_curves) { + for(Curve_index curve_index : incident_curves) + { if(domain.is_loop(curve_index)) { FT curve_length = d_ptr->domain.curve_length(curve_index); + auto loc = + d_ptr->domain.locate_corner(curve_index, pair.first); Point_3 other_point = d_ptr->domain.construct_point_on_curve(pair.first, curve_index, - curve_length / 2); + curve_length / 2, + loc); d_ptr->dt.insert(other_point); } } From 789d4c4f9d9262a10721aa678599ce29f5c77adc Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 13 Nov 2025 15:38:20 +0100 Subject: [PATCH 13/19] wip using locate_corner() --- Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h index 74346f5ce94..2471f6028b1 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -228,7 +228,8 @@ public: //fill incidences of corners with curves d_ptr->corners_incident_curves.resize(d_ptr->corners.size()); - for(const typename Corners_indices::value_type& pair : d_ptr->corners_indices) { + for(const typename Corners_indices::value_type& pair : d_ptr->corners_indices) + { d_ptr->dt.insert(pair.first); // Fill `corners_incident_curves[corner_id]` @@ -236,15 +237,19 @@ public: d_ptr->domain.get_corner_incident_curves(pair.second, std::inserter(incident_curves, incident_curves.end())); - // For each incident loops, insert a point on the loop, as far as + // For each incident loop, insert a point on the loop, as far as // possible. - for(Curve_index curve_index : incident_curves) { + for(Curve_index curve_index : incident_curves) + { if(domain.is_loop(curve_index)) { FT curve_length = d_ptr->domain.curve_length(curve_index); + auto loc = + d_ptr->domain.locate_corner(curve_index, pair.first); Point_3 other_point = d_ptr->domain.construct_point_on_curve(pair.first, curve_index, - curve_length / 2); + curve_length / 2, + loc); d_ptr->dt.insert(other_point); } } From b44e1325653a167722248e74c253f4b4b6ed4600 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 14 Nov 2025 10:34:02 +0100 Subject: [PATCH 14/19] fix compilation for return type --- Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h index 2471f6028b1..db291ec470f 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -245,7 +245,7 @@ public: FT curve_length = d_ptr->domain.curve_length(curve_index); auto loc = d_ptr->domain.locate_corner(curve_index, pair.first); - Point_3 other_point = + auto [other_point, _] = d_ptr->domain.construct_point_on_curve(pair.first, curve_index, curve_length / 2, From 5c5c48ba1676ad8bb0d1de65baa31e926fbcd1eb Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 14 Nov 2025 10:46:12 +0100 Subject: [PATCH 15/19] improve readability --- Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h index db291ec470f..f21ed1ea621 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -228,13 +228,13 @@ public: //fill incidences of corners with curves d_ptr->corners_incident_curves.resize(d_ptr->corners.size()); - for(const typename Corners_indices::value_type& pair : d_ptr->corners_indices) + for(const auto& [corner_pt, corner_index] : d_ptr->corners_indices) { - d_ptr->dt.insert(pair.first); + d_ptr->dt.insert(corner_pt); // Fill `corners_incident_curves[corner_id]` - Curves_ids& incident_curves = d_ptr->corners_incident_curves[pair.second]; - d_ptr->domain.get_corner_incident_curves(pair.second, + Curves_ids& incident_curves = d_ptr->corners_incident_curves[corner_index]; + d_ptr->domain.get_corner_incident_curves(corner_index, std::inserter(incident_curves, incident_curves.end())); // For each incident loop, insert a point on the loop, as far as @@ -244,9 +244,9 @@ public: if(domain.is_loop(curve_index)) { FT curve_length = d_ptr->domain.curve_length(curve_index); auto loc = - d_ptr->domain.locate_corner(curve_index, pair.first); + d_ptr->domain.locate_corner(curve_index, corner_pt); auto [other_point, _] = - d_ptr->domain.construct_point_on_curve(pair.first, + d_ptr->domain.construct_point_on_curve(corner_pt, curve_index, curve_length / 2, loc); From e74e5de28a5374cd42e55e22dd12924442678f51 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Mon, 17 Nov 2025 12:32:00 +0100 Subject: [PATCH 16/19] move map to Mesh_domain_with_polyline_features_3 and use it in Protect_edges_sizing_field --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 80 ++++++++--------- .../Mesh_domain_with_polyline_features_3.h | 85 +++++++++++++++++-- .../CGAL/Sizing_field_with_aabb_tree.h | 7 +- 3 files changed, 121 insertions(+), 51 deletions(-) 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 d76d8c5854f..2a2cd85960b 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 @@ -494,24 +494,7 @@ private: return use_edge_distance_; } - Polyline_iterator locate_in_polyline(Vertex_handle vh, const Curve_index& index) const - { - CGAL_assertion(vh->in_dimension() < 2); - if(vh->in_dimension() == 0) //corner - { - auto cp = c3t3_.triangulation().geom_traits().construct_point_3_object(); - return domain_.locate_corner(index, cp(vh->point())); - } - else - return vertex_to_polyline_iterator_.at(vh); - } - - void set_polyline_iterator(Vertex_handle vh, Polyline_iterator it) - { - CGAL_assertion(vh->in_dimension() == 1); - vertex_to_polyline_iterator_[vh] = it; - } private: C3T3& c3t3_; @@ -526,7 +509,6 @@ private: int refine_balls_iteration_nb; bool nonlinear_growth_of_balls; const std::size_t maximal_number_of_vertices_; - std::unordered_map vertex_to_polyline_iterator_; Mesh_error_code* const error_code_; #ifndef CGAL_NO_ATOMIC /// Pointer to the atomic Boolean that can stop the process @@ -1090,7 +1072,7 @@ insert_balls_on_edges() p_index, Vertex_handle(), CGAL::Emptyset_iterator()).first; - set_polyline_iterator(vp, p_polyline_iter); + domain_.set_polyline_iterator(p, p_polyline_iter); } // No 'else' because in that case 'is_vertex(..)' already filled @@ -1264,11 +1246,12 @@ insert_balls(const Vertex_handle& vp, curve_index, d_sign) << ")\n"; #endif - const auto [bp, polyline_iter] = - domain_.construct_point_on_curve(cp(vp_wp), + const Bare_point p = cp(vp_wp); + const auto [bp, polyline_iter] = //[Bare_point, Polyline_const_iterator] + domain_.construct_point_on_curve(p, curve_index, d_signF * d / 2, - locate_in_polyline(vp, curve_index)); + domain_.locate_in_polyline(p, vp->in_dimension(), curve_index)); const Bare_point new_point = bp; const int dim = 1; // new_point is on edge const Index index = domain_.index_from_curve_index(curve_index); @@ -1285,7 +1268,7 @@ insert_balls(const Vertex_handle& vp, Vertex_handle(), out); const Vertex_handle new_vertex = pair.first; - set_polyline_iterator(new_vertex, polyline_iter); + domain_.set_polyline_iterator(new_point, polyline_iter); if(forced_stop()) return out; out = pair.second; @@ -1363,7 +1346,7 @@ insert_balls(const Vertex_handle& vp, const Index index = domain_.index_from_curve_index(curve_index); // Launch balls - Polyline_iterator p_loc = locate_in_polyline(vp, curve_index); + Polyline_iterator p_loc = domain_.locate_in_polyline(p, vp->in_dimension(), curve_index); CGAL_assertion(p_loc == domain_.locate_point(curve_index, p)); Bare_point prev_pt = p; FT dist_to_prev = pt_dist; @@ -1383,7 +1366,7 @@ insert_balls(const Vertex_handle& vp, smart_insert_point(new_point, point_weight, dim, index, prev, out); Vertex_handle new_vertex = pair.first; out = pair.second; - set_polyline_iterator(new_vertex, polyline_iter); + domain_.set_polyline_iterator(new_point, polyline_iter); // Add edge to c3t3 if(!c3t3_.is_in_complex(prev, new_vertex)) { @@ -1600,8 +1583,8 @@ approx_is_too_large(const Edge& e, const bool is_edge_in_complex) const const Bare_point& pb = vb->point().point(); const Curve_index curve_index = c3t3_.curve_index(e); - Polyline_iterator pa_it = locate_in_polyline(va, curve_index); - Polyline_iterator pb_it = locate_in_polyline(vb, curve_index); + Polyline_iterator pa_it = domain_.locate_in_polyline(pa, va->in_dimension(), curve_index); + Polyline_iterator pb_it = domain_.locate_in_polyline(pb, vb->in_dimension(), curve_index); // Construct the geodesic middle point const FT signed_geodesic_distance @@ -1875,14 +1858,14 @@ curve_segment_length(const Vertex_handle v1, v2_valid_curve_index = (domain_.curve_index(v2->index()) == curve_index); } - const Weighted_point& v1_wp = c3t3_.triangulation().point(v1); - const Weighted_point& v2_wp = c3t3_.triangulation().point(v2); + const Bare_point p1 = cp(c3t3_.triangulation().point(v1)); + const Bare_point p2 = cp(c3t3_.triangulation().point(v2)); FT arc_length = (v1_valid_curve_index && v2_valid_curve_index) - ? domain_.curve_segment_length(cp(v1_wp), - cp(v2_wp), - locate_in_polyline(v1, curve_index), - locate_in_polyline(v2, curve_index), + ? domain_.curve_segment_length(p1, + p2, + domain_.locate_in_polyline(p1, v1->in_dimension(), curve_index), + domain_.locate_in_polyline(p2, v2->in_dimension(), curve_index), curve_index, orientation) : compute_distance(v1, v2); //curve polyline may not be consistent @@ -1934,12 +1917,15 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2, const Weighted_point& v1_wp = c3t3_.triangulation().point(v1); const Weighted_point& v2_wp = c3t3_.triangulation().point(v2); + const Bare_point p1 = cp(v1_wp); + const Bare_point p2 = cp(v2_wp); + const bool cov = domain_.is_curve_segment_covered(curve_index, - orientation, - cp(v1_wp), cp(v2_wp), - cw(v1_wp), cw(v2_wp), - locate_in_polyline(v1, curve_index), - locate_in_polyline(v2, curve_index)); + orientation, + p1, p2, + cw(v1_wp), cw(v2_wp), + domain_.locate_in_polyline(p1, v1->in_dimension(), curve_index), + domain_.locate_in_polyline(p2, v2->in_dimension(), curve_index)); #if CGAL_MESH_3_PROTECTION_DEBUG & 1 if(cov) { @@ -1973,20 +1959,28 @@ orientation_of_walk(const Vertex_handle& start, const Weighted_point& start_wp = c3t3_.triangulation().point(start); const Weighted_point& next_wp = c3t3_.triangulation().point(next); + const Bare_point start_p = cp(start_wp); + const Bare_point next_p = cp(next_wp); if(domain_.is_loop(curve_index)) { // if the curve is a cycle, the direction is the direction passing // through the next vertex, and the next-next vertex Vertex_handle next_along_curve = next_vertex_along_curve(next,start,curve_index); const Weighted_point& next_along_curve_wp = c3t3_.triangulation().point(next_along_curve); + const Bare_point next_along_curve_p = cp(next_along_curve_wp); return domain_.distance_sign_along_loop( - cp(start_wp), cp(next_wp), cp(next_along_curve_wp), curve_index); - } else { + start_p, next_p, next_along_curve_p, curve_index, + domain_.locate_in_polyline(start_p, start->in_dimension(), curve_index), + domain_.locate_in_polyline(next_p, next->in_dimension(), curve_index), + domain_.locate_in_polyline(next_along_curve_p, next_along_curve->in_dimension(), curve_index)); + } + else + { // otherwise, the sign is just the sign of the geodesic distance - return domain_.distance_sign(cp(start_wp), cp(next_wp), curve_index, - locate_in_polyline(start, curve_index), - locate_in_polyline(next, curve_index)); + return domain_.distance_sign(start_p, next_p, curve_index, + domain_.locate_in_polyline(start_p, start->in_dimension(), curve_index), + domain_.locate_in_polyline(next_p, next->in_dimension(), curve_index)); } } diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index bf14dcb9b51..416569741dd 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -36,6 +36,7 @@ #include #include #include +#include #include #include @@ -218,6 +219,7 @@ public: /// returns the length of the polyline FT length() const { + length_timer.start(); if(length_ < 0.) { FT result(0); @@ -229,6 +231,7 @@ public: } length_ = result; } + length_timer.stop(); return length_; } @@ -445,20 +448,27 @@ private: const_iterator locate(const Point_3& p, bool end_point_first=false) const { CGAL_precondition(is_valid()); + old_locate_timer.start(); // First look if p is one of the points of the polyline const_iterator result = std::find(points_.begin(), points_.end(), p); if ( result != points_.end() ) { - if ( result != points_.begin() ) - { return --result; } + if ( result != points_.begin() ) { + old_locate_timer.stop(); + return --result; + } else { // Treat loops - if ( end_point_first && p == end_point() ) - { return last_segment_source(); } - else - { return result; } + if ( end_point_first && p == end_point() ) { + old_locate_timer.stop(); + return last_segment_source(); + } + else { + old_locate_timer.stop(); + return result; + } } } @@ -515,7 +525,7 @@ private: previous = it; } // end the while loop on the vertices of the polyline - + old_locate_timer.stop(); if(result == points_.begin()) { return (end_point_first && !nearest_is_a_segment) ? last_segment_source() : points_.begin(); } else { @@ -857,6 +867,15 @@ public: const Point_3& r, const Curve_index& index) const; + CGAL::Sign distance_sign_along_loop(const Point_3& p, + const Point_3& q, + const Point_3& r, + const Curve_index& index, + Polyline_const_iterator pit, + Polyline_const_iterator qit, + Polyline_const_iterator rit) const; + + /// implements `MeshDomainWithFeatures_3::distance_sign()`. CGAL::Sign distance_sign(const Point_3& p, const Point_3& q, const Curve_index& index, @@ -1012,6 +1031,7 @@ public: private: mutable std::shared_ptr curves_aabb_tree_ptr_; mutable bool curves_aabb_tree_is_built; + mutable std::unordered_map vertex_to_polyline_iterator_; public: const Corners_incidences& corners_incidences_map() const @@ -1059,6 +1079,29 @@ public: #endif } // build_curves_aabb_tree() + Polyline_const_iterator locate_in_polyline(const Point_3& p, + const int dim, + const Curve_index& index) const + { + ++locate_counter; + + CGAL_assertion(dim < 2); + + Polyline_const_iterator it; + if(dim == 0) // corner + it = locate_corner(index, p); + else + it = vertex_to_polyline_iterator_.at(p); + + //locate_timer.stop(); + return it; + } + + void set_polyline_iterator(const Point_3& p, Polyline_const_iterator it) const + { + vertex_to_polyline_iterator_[p] = it; + } + void dump_curve(const Curve_index& index, const std::string& prefix) const { std::string filename(prefix); @@ -1733,6 +1776,34 @@ distance_sign_along_loop(const Point_3& p, else { return CGAL::NEGATIVE; } } +template +CGAL::Sign Mesh_domain_with_polyline_features_3:: +distance_sign_along_loop(const Point_3& p, + const Point_3& q, + const Point_3& r, + const Curve_index& index, + Polyline_const_iterator pit, + Polyline_const_iterator qit, + Polyline_const_iterator rit) const +{ + CGAL_assertion(p != q); + CGAL_assertion(p != r); + CGAL_assertion(r != q); + + // Find edge + typename Edges::const_iterator eit = edges_.find(index); + CGAL_assertion(eit != edges_.end()); + CGAL_assertion(eit->second.is_loop()); + + FT pq = eit->second.curve_segment_length(p,q,CGAL::POSITIVE,pit,qit); + FT pr = eit->second.curve_segment_length(p,r,CGAL::POSITIVE,pit,rit); + + // Compare pq and pr + if ( pq <= pr ) { return CGAL::POSITIVE; } else { + return CGAL::NEGATIVE; + } +} + template bool Mesh_domain_with_polyline_features_3:: diff --git a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h index f21ed1ea621..cbbbdd5f4cf 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -575,7 +575,12 @@ public: if (const Point_3* pp = std::get_if(&*int_res)) { FT new_sqd = CGAL::squared_distance(p, *pp); - FT dist = CGAL::abs(d_ptr->domain.signed_geodesic_distance(p, *pp, curve_id)); + auto p_polyline_const_it = ppid.second.second; + auto pp_polyline_const_it = prim.id().second; + FT dist = CGAL::abs(d_ptr->domain.signed_geodesic_distance(p, *pp, + p_polyline_const_it, + pp_polyline_const_it, + curve_id)); #ifdef CGAL_MESH_3_PROTECTION_HIGH_VERBOSITY std::cerr << "Intersection point : Point_3(" << *pp << ") "; From 6ec5daa6b99c02009a4cef854933f913a23ace34 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Mon, 17 Nov 2025 12:46:47 +0100 Subject: [PATCH 17/19] cleaning --- .../Mesh_domain_with_polyline_features_3.h | 35 +++++-------------- 1 file changed, 9 insertions(+), 26 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 416569741dd..8dd91482c48 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -219,7 +219,6 @@ public: /// returns the length of the polyline FT length() const { - length_timer.start(); if(length_ < 0.) { FT result(0); @@ -231,7 +230,6 @@ public: } length_ = result; } - length_timer.stop(); return length_; } @@ -448,27 +446,20 @@ private: const_iterator locate(const Point_3& p, bool end_point_first=false) const { CGAL_precondition(is_valid()); - old_locate_timer.start(); // First look if p is one of the points of the polyline const_iterator result = std::find(points_.begin(), points_.end(), p); if ( result != points_.end() ) { - if ( result != points_.begin() ) { - old_locate_timer.stop(); - return --result; - } + if ( result != points_.begin() ) + { return --result; } else { // Treat loops - if ( end_point_first && p == end_point() ) { - old_locate_timer.stop(); - return last_segment_source(); - } - else { - old_locate_timer.stop(); - return result; - } + if ( end_point_first && p == end_point() ) + { return last_segment_source(); } + else + { return result; } } } @@ -525,7 +516,6 @@ private: previous = it; } // end the while loop on the vertices of the polyline - old_locate_timer.stop(); if(result == points_.begin()) { return (end_point_first && !nearest_is_a_segment) ? last_segment_source() : points_.begin(); } else { @@ -533,16 +523,12 @@ private: } } - // FT squared_distance(const Point_3& p, const Point_3& q) const - // { - // typename Kernel::Compute_squared_distance_3 sq_distance = - // Kernel().compute_squared_distance_3_object(); - // return sq_distance(p,q); - // } FT distance(const Point_3& p, const Point_3& q) const { - return CGAL::sqrt(squared_distance(p, q)); + typename Kernel::Compute_squared_distance_3 sq_distance = + Kernel().compute_squared_distance_3_object(); + return CGAL::sqrt(sq_distance(p, q)); } Angle angle(const Point_3& p, @@ -1083,8 +1069,6 @@ public: const int dim, const Curve_index& index) const { - ++locate_counter; - CGAL_assertion(dim < 2); Polyline_const_iterator it; @@ -1093,7 +1077,6 @@ public: else it = vertex_to_polyline_iterator_.at(p); - //locate_timer.stop(); return it; } From b804ec9394d3b059121ea7c1233f44ac884aedae Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 18 Nov 2025 10:10:54 +0100 Subject: [PATCH 18/19] duplicate functions to keep original API valid --- .../Mesh_domain_with_polyline_features_3.h | 57 ++++++++++++++++++- .../Protect_edges_sizing_field.h | 32 +++++++---- 2 files changed, 76 insertions(+), 13 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 8dd91482c48..348e765198a 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -157,6 +157,14 @@ public: return cover_pred(s1, s2, *c2_it, c2); } + bool is_curve_segment_covered(CGAL::Orientation orientation, + const Point_3& c1, const Point_3& c2, + const FT sq_r1, const FT sq_r2) const + { + return is_curve_segment_covered(orientation, + c1, c2, sq_r1, sq_r2, locate(c1), locate(c2)); + } + FT curve_segment_length(const Point_3& p, const Point_3 q, CGAL::Orientation orientation) const { @@ -426,6 +434,11 @@ public: return ( compare_distance(*pit,p,q) != CGAL::LARGER ); } + bool are_ordered_along(const Point_3& p, const Point_3& q) const + { + return are_ordered_along(p, q, locate(p), locate(q, true)); + } + private: const_iterator first_segment_source() const { @@ -863,7 +876,12 @@ public: /// implements `MeshDomainWithFeatures_3::distance_sign()`. - CGAL::Sign distance_sign(const Point_3& p, const Point_3& q, + CGAL::Sign distance_sign(const Point_3& p, + const Point_3& q, + const Curve_index& index) const; + + CGAL::Sign distance_sign(const Point_3& p, + const Point_3& q, const Curve_index& index, Polyline_const_iterator pit, Polyline_const_iterator qit) const; @@ -872,6 +890,11 @@ public: bool is_loop(const Curve_index& index) const; /// implements `MeshDomainWithFeatures_3::is_curve_segment_covered()`. + bool is_curve_segment_covered(const Curve_index& index, + CGAL::Orientation orientation, + const Point_3& c1, const Point_3& c2, + const FT sq_r1, const FT sq_r2) const; + bool is_curve_segment_covered(const Curve_index& index, CGAL::Orientation orientation, const Point_3& c1, const Point_3& c2, @@ -1733,6 +1756,24 @@ distance_sign(const Point_3& p, const Point_3& q, return CGAL::NEGATIVE; } +template +CGAL::Sign +Mesh_domain_with_polyline_features_3:: +distance_sign(const Point_3& p, + const Point_3& q, + const Curve_index& curve_index) const +{ + typename Edges::const_iterator eit = edges_.find(curve_index); + CGAL_assertion(eit != edges_.end()); + CGAL_precondition(!eit->second.is_loop()); + + if(p == q) + return CGAL::ZERO; + else if(eit->second.are_ordered_along(p, q)) + return CGAL::POSITIVE; + else + return CGAL::NEGATIVE; +} template CGAL::Sign @@ -1818,6 +1859,20 @@ is_curve_segment_covered(const Curve_index& index, c1_it, c2_it); } +template +bool +Mesh_domain_with_polyline_features_3:: +is_curve_segment_covered(const Curve_index& index, + CGAL::Orientation orientation, + const Point_3& c1, const Point_3& c2, + const FT sq_r1, const FT sq_r2) const +{ + typename Edges::const_iterator eit = edges_.find(index); + CGAL_assertion(eit != edges_.end()); + + return eit->second.is_curve_segment_covered(orientation, c1, c2, sq_r1, sq_r2); +} + template typename Mesh_domain_with_polyline_features_3::Polyline_const_iterator Mesh_domain_with_polyline_features_3:: diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h index fd1e0b9e104..33f1ddc2c3e 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h @@ -98,6 +98,8 @@ public: typedef typename MeshDomain::Corner_index Corner_index; typedef typename MeshDomain::Index Index; + using Polyline_iterator = typename CGAL::Mesh_3::internal::Polyline::const_iterator; + private: typedef typename CGAL::Kernel_traits::Kernel Kernel; @@ -1884,35 +1886,40 @@ Protect_edges_sizing_field:: insert_balls_on_edges() { // Get features - typedef std::tuple, - std::pair > Feature_tuple; - typedef std::vector Input_features; + // Get features + struct Feature_tuple + { + Curve_index curve_index_; + Polyline_iterator polyline_begin_; + std::pair point_s_; + std::pair point_t_; + }; + typedef std::vector Input_features; Input_features input_features; domain_.get_curves(std::back_inserter(input_features)); // Iterate on edges - for(typename Input_features::iterator fit = input_features.begin(), - end = input_features.end() ; fit != end ; ++fit) + for(const Feature_tuple& ft : input_features) { - const Curve_index& curve_index = std::get<0>(*fit); + const Curve_index& curve_index = ft.curve_index_; if(! is_treated(curve_index)) { #if CGAL_MESH_3_PROTECTION_DEBUG & 1 std::cerr << "** treat curve #" << curve_index << std::endl; std::cerr << "is it a loop? " << domain_.is_loop(curve_index) << std::endl; #endif - const Bare_point& p = std::get<1>(*fit).first; - const Bare_point& q = std::get<2>(*fit).first; - - const Index& p_index = std::get<1>(*fit).second; - const Index& q_index = std::get<2>(*fit).second; + const Bare_point& p = ft.point_s_.first; + const Index& p_index = ft.point_s_.second; + const Polyline_iterator& p_polyline_iter = ft.polyline_begin_; Vertex_handle vp,vq; if(! domain_.is_loop(curve_index)) { vp = get_vertex_corner_from_point(p, p_index); + + const Bare_point& q = ft.point_t_.first; + const Index& q_index = ft.point_t_.second; vq = get_vertex_corner_from_point(q, q_index); } else @@ -1943,6 +1950,7 @@ insert_balls_on_edges() p_index, curve_index, CGAL::Emptyset_iterator()).first; + domain_.set_polyline_iterator(p, p_polyline_iter); } // No 'else' because in that case 'is_vertex(..)' already filled // the variable 'vp'. From 5b574bfaae42a763506cfbaf440e7d813733f0b5 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 18 Nov 2025 12:08:25 +0100 Subject: [PATCH 19/19] move CGAL::Mesh_3::internal::Polyline class to its own header --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 2 +- .../include/CGAL/Mesh_3/internal/Polyline.h | 561 ++++++++++++++++++ .../Mesh_domain_with_polyline_features_3.h | 525 +--------------- .../Protect_edges_sizing_field.h | 1 + 4 files changed, 564 insertions(+), 525 deletions(-) create mode 100644 Mesh_3/include/CGAL/Mesh_3/internal/Polyline.h 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 2a2cd85960b..c160948ffa8 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 @@ -46,7 +46,7 @@ #include #include #include -#include +#include #include diff --git a/Mesh_3/include/CGAL/Mesh_3/internal/Polyline.h b/Mesh_3/include/CGAL/Mesh_3/internal/Polyline.h new file mode 100644 index 00000000000..73589c34ba2 --- /dev/null +++ b/Mesh_3/include/CGAL/Mesh_3/internal/Polyline.h @@ -0,0 +1,561 @@ +// Copyright (c) 2009-2010 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Stéphane Tayeb, Laurent Rineau +// +//****************************************************************************** +// File Description : +// +//****************************************************************************** + +#ifndef CGAL_MESH_3_INTERNAL_POLYLINE_H +#define CGAL_MESH_3_INTERNAL_POLYLINE_H + +#include + +#include + +#include + +namespace CGAL { + +/// @cond DEVELOPERS +namespace Mesh_3 { +namespace internal { + +template +class Polyline +{ + typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Segment_3 Segment_3; + typedef typename Kernel::FT FT; + + typedef std::vector Data; + +public: + typedef typename Data::const_iterator const_iterator; + typedef std::pair Point_and_location; + + Polyline() {} + ~Polyline() {} + + /// adds a point at the end of the polyline + void add_point(const Point_3& p) + { + if( points_.empty() || p != end_point() ) { + points_.push_back(p); + } + } + + /// returns the starting point of the polyline + const Point_3& start_point() const + { + CGAL_assertion( ! points_.empty() ); + return points_.front(); + } + + /// returns the ending point of the polyline + const Point_3& end_point() const + { + CGAL_assertion( ! points_.empty() ); + return points_.back(); + } + + /// returns `true` if the polyline is not degenerate + bool is_valid() const + { + return points_.size() > 1; + } + + /// returns `true` if polyline is a loop + bool is_loop() const + { + return start_point() == end_point(); + } + + const_iterator next(const_iterator it, Orientation orientation) const { + if(orientation == POSITIVE) { + CGAL_assertion(it != (points_.end() - 1)); + if(it == (points_.end() - 2)) { + CGAL_assertion(is_loop()); + it = points_.begin(); + } else { + ++it; + } + } else { + CGAL_assertion(orientation == NEGATIVE); + CGAL_assertion(it != points_.begin()); + if(it == (points_.begin() + 1)) { + CGAL_assertion(is_loop()); + it = points_.end() - 1; + } else { + --it; + } + } + return it; + } + + bool is_curve_segment_covered(CGAL::Orientation orientation, + const Point_3& c1, const Point_3& c2, + const FT sq_r1, const FT sq_r2, + const const_iterator cc1_it, + const const_iterator cc2_it) const + { + CGAL_assertion(orientation != CGAL::ZERO); + typename Kernel::Has_on_bounded_side_3 cover_pred = + Kernel().has_on_bounded_side_3_object(); + + typedef typename Kernel::Sphere_3 Sphere_3; + const Sphere_3 s1(c1, sq_r1); + const Sphere_3 s2(c2, sq_r2); + + const_iterator c1_it = cc1_it; + const_iterator c2_it = cc2_it; + + if(orientation == CGAL::NEGATIVE) { + ++c1_it; + ++c2_it; + CGAL_assertion(c1_it != points_.end()); + CGAL_assertion(c2_it != points_.end()); + } + + if(c1_it == c2_it) return cover_pred(s1, s2, c1, c2); + const_iterator next_it = this->next(c1_it, orientation); + + if(!cover_pred(s1, s2, c1, *next_it)) return false; + + for(const_iterator it = next_it; it != c2_it; /* in body */) { + next_it = this->next(it, orientation); + if(!cover_pred(s1, s2, *it, *next_it)) return false; + it = next_it; + } // end loop ]c1_it, c2_it[ + + return cover_pred(s1, s2, *c2_it, c2); + } + + bool is_curve_segment_covered(CGAL::Orientation orientation, + const Point_3& c1, const Point_3& c2, + const FT sq_r1, const FT sq_r2) const + { + return is_curve_segment_covered(orientation, + c1, c2, sq_r1, sq_r2, locate(c1), locate(c2)); + } + + FT curve_segment_length(const Point_3& p, const Point_3 q, + CGAL::Orientation orientation) const + { + CGAL_assertion(orientation != CGAL::ZERO); + const_iterator p_it = locate(p); + const_iterator q_it = locate(q); + return curve_segment_length(p, q, orientation, p_it, q_it); + } + + FT curve_segment_length(const Point_3& p, const Point_3 q, + CGAL::Orientation orientation, + const_iterator p_it, + const_iterator q_it) const + { + CGAL_assertion(orientation != CGAL::ZERO); + CGAL_assertion(p_it == locate(p)); + CGAL_assertion(q_it == locate(q)); + + if(p_it == q_it) { + const CGAL::Comparison_result cmp = compare_distance(*p_it,p,q); + if( (cmp != LARGER && orientation == POSITIVE) || + (cmp != SMALLER && orientation == NEGATIVE) ) + { + // If the orientation of `p` and `q` on the segment is compatible + // with `orientation`, then return the distance between the two + // points. + return distance(p, q); + } + } + + if(orientation == CGAL::NEGATIVE) { + ++p_it; + ++q_it; + CGAL_assertion(p_it != points_.end()); + CGAL_assertion(q_it != points_.end()); + } + + const_iterator next_it = this->next(p_it, orientation); + FT result = distance(p, *next_it); + for(const_iterator it = next_it; it != q_it; /* in body */) { + next_it = this->next(it, orientation); + result += distance(*it, *next_it); + it = next_it; + } // end loop ]p_it, q_it[ + result += distance(*q_it, q); + return result; + } + + + /// returns the angle at the first point. + /// \pre The polyline must be a loop. + Angle angle_at_first_point() const { + CGAL_precondition(is_loop()); + const Point_3& first = points_.front(); + const Point_3& next_p = points_[1]; + const Point_3& prev = points_[points_.size() - 2]; + return angle(prev, first, next_p); + } + + /// returns the length of the polyline + FT length() const + { + if(length_ < 0.) + { + FT result(0); + const_iterator it = points_.begin(); + const_iterator previous = it++; + + for(const_iterator end = points_.end(); it != end; ++it, ++previous) { + result += distance(*previous, *it); + } + length_ = result; + } + return length_; + } + + /// returns the signed geodesic distance between `p` and `q`. + FT signed_geodesic_distance(const Point_3& p, const Point_3& q) const + { + // Locate p & q on polyline + const_iterator pit = locate(p); + const_iterator qit = locate(q,false); + + return signed_geodesic_distance(p, q, pit, qit); + } + + FT signed_geodesic_distance(const Point_3& p, const Point_3& q, + const_iterator pit, const_iterator qit) const + { + CGAL_assertion(pit == locate(p)); + CGAL_assertion(qit == locate(q, false)); + + // If p and q are in the same segment of the polyline + if ( pit == qit ) + { + FT result = distance(p,q); + + // Find the closest point to *pit + if ( compare_distance(*pit,p,q) != CGAL::LARGER ) + { return result; } + else + { return -result; } + } + if(is_loop()) + { + FT positive_distance, negative_distance; + if(pit <= qit) + { + positive_distance = curve_segment_length(p, q, CGAL::POSITIVE, pit, qit); + negative_distance = length() - positive_distance; + } + else + { + negative_distance = curve_segment_length(q, p, CGAL::POSITIVE, qit, pit); + positive_distance = length() - negative_distance; + } + return (positive_distance < negative_distance) ? positive_distance : (-negative_distance); + } + else + { + return (pit <= qit) + ? curve_segment_length(p, q, CGAL::POSITIVE, pit, qit) + : ( - curve_segment_length(p, q, CGAL::NEGATIVE, pit, qit) ); + } + } + + const_iterator previous_segment_source(const_iterator it) const + { + CGAL_assertion(it != points_.end()); + if(it == first_segment_source()) + { + CGAL_assertion(is_loop()); + it = last_segment_source(); + } + else + { + --it; + } + return it; + } + + const_iterator next_segment_source(const_iterator it) const + { + CGAL_assertion(it != points_.end()); + if(it == last_segment_source()) + { + if(is_loop()) + return first_segment_source(); + else + return points_.end(); + } + else + { + ++it; + return it; + } + } + + /// returns a point at geodesic distance `distance` from p along the + /// polyline. The polyline is oriented from starting point to end point. + /// The distance could be negative. + Point_3 point_at(const Point_3& start_pt, FT distance) const + { + return point_at(start_pt, distance, locate(start_pt)).first; + } + + /// returns a point at geodesic distance `distance` from `start_pt` along the + /// polyline. The polyline is oriented from starting point to end point. + /// The distance could be negative. + Point_and_location point_at(const Point_3& start_pt, + FT distance, + const_iterator start_it) const + { + CGAL_assertion(start_it == locate(start_pt)); + + const Point_3& start_it_pt = *start_it; + const_iterator start_it_locate_pt + = (start_it == points_.begin()) ? start_it : std::prev(start_it); + + distance += curve_segment_length(start_it_pt, start_pt, CGAL::POSITIVE, + start_it_locate_pt, start_it); + + // If polyline is a loop, ensure that distance is given from start_it + if(is_loop()) + { + if(distance < FT(0)) { distance += length(); } + else if(distance > length()) { distance -= length(); } + } + else if(distance < FT(0)) // If polyline is not a loop and distance is < 0, go backward + { + Point_3 new_start = start_pt; + while(distance < FT(0)) + { + start_it = previous_segment_source(start_it); + distance += this->distance(new_start, *start_it); + new_start = *start_it; + } + } + + CGAL_assertion(distance >= FT(0)); + CGAL_assertion(distance <= length()); + + // Initialize iterators + const_iterator pit = start_it; // start at start_it, and walk forward + const_iterator previous = pit++; + + // Iterate to find which segment contains the point we want to construct + FT segment_length = this->distance(*previous, *pit); + while(distance > segment_length) + { + distance -= segment_length; + + // Increment iterators and update length + previous = next_segment_source(previous); + pit = next_segment_source(pit); + + if(pit == points_.end()) + { + CGAL_assertion(distance < this->distance(*previous, end_point())); + break; // return {*previous, previous} + } + else + segment_length = this->distance(*previous, *pit); + } + + // return point at distance from current segment source + using Vector_3 = typename Kernel::Vector_3; + auto vector = Kernel().construct_vector_3_object(); + Vector_3 v = (pit != points_.end()) ? vector(*previous, *pit) + : vector(*previous, end_point()); + + return {(*previous) + (distance / CGAL::sqrt(v.squared_length())) * v, + previous}; + } + + const_iterator locate_corner(const Point_3& p) const + { + const_iterator res = points_.end(); + if(p == start_point()) + res = points_.begin(); + else if(p == end_point()) + res = last_segment_source(); + + CGAL_assertion(res == locate(p)); + CGAL_assertion(res != points_.end()); + return res; + } + + const_iterator locate_point(const Point_3& p) const + { + return locate(p); + } + + bool are_ordered_along(const Point_3& p, const Point_3& q, + const_iterator pit, const_iterator qit) const + { + CGAL_precondition(!is_loop()); + + // Locate p & q on polyline + CGAL_assertion(pit == locate(p)); + CGAL_assertion(qit == locate(q, true)); + + // Points are not located on the same segment + if ( pit != qit ) { return (pit <= qit); } + + // pit == qit, then we have to sort p&q along (pit,pit+1) + return ( compare_distance(*pit,p,q) != CGAL::LARGER ); + } + + bool are_ordered_along(const Point_3& p, const Point_3& q) const + { + return are_ordered_along(p, q, locate(p), locate(q, true)); + } + +private: + const_iterator first_segment_source() const + { + CGAL_precondition(is_valid()); + return points_.begin(); + } + + const_iterator last_segment_source() const + { + CGAL_precondition(is_valid()); + return (points_.end() - 2); + } + + /// returns an iterator on the starting point of the segment of the + /// polyline which contains p + /// if end_point_first is true, then --end is returned instead of begin + /// if p is the starting point of a loop. + const_iterator locate(const Point_3& p, bool end_point_first=false) const + { + CGAL_precondition(is_valid()); + + // First look if p is one of the points of the polyline + const_iterator result = std::find(points_.begin(), points_.end(), p); + if ( result != points_.end() ) + { + if ( result != points_.begin() ) + { return --result; } + else + { + // Treat loops + if ( end_point_first && p == end_point() ) + { return last_segment_source(); } + else + { return result; } + } + } + + CGAL_assertion(result == points_.end()); + + // Get result by projecting p on the polyline + const_iterator it = points_.begin(); + const_iterator previous = it; + Segment_3 nearest_segment; + const_iterator nearest_vertex = it; + result = nearest_vertex; + bool nearest_is_a_segment = false; + + while ( ++it != points_.end() ) + { + Segment_3 seg (*previous, *it); + + if(nearest_is_a_segment) + { + if(compare_distance(p, *it, nearest_segment) == CGAL::SMALLER) + { + nearest_vertex = it; + nearest_is_a_segment = false; + result = it; + if (possibly(angle(*previous, *it, p) == CGAL::ACUTE) && + compare_distance(p, seg, *nearest_vertex) == CGAL::SMALLER) + { + nearest_segment = seg; + nearest_is_a_segment = true; + result = previous; + } + } + else if(compare_distance(p, seg, nearest_segment) == CGAL::SMALLER) + { + nearest_segment = seg; + result = previous; + } + } + else { + if(compare_distance(p, *it, *nearest_vertex) == CGAL::SMALLER) + { + nearest_vertex = it; + result = it; + } + if ((nearest_vertex != it || + possibly(angle(*previous, *it, p) == CGAL::ACUTE)) && + compare_distance(p, seg, *nearest_vertex) == CGAL::SMALLER) + { + nearest_segment = seg; + nearest_is_a_segment = true; + result = previous; + } + } + previous = it; + } // end the while loop on the vertices of the polyline + + if(result == points_.begin()) { + return (end_point_first && !nearest_is_a_segment) ? last_segment_source() : points_.begin(); + } else { + return result; + } + } + + + FT distance(const Point_3& p, const Point_3& q) const + { + typename Kernel::Compute_squared_distance_3 sq_distance = + Kernel().compute_squared_distance_3_object(); + return CGAL::sqrt(sq_distance(p, q)); + } + + Angle angle(const Point_3& p, + const Point_3& angle_vertex_point, + const Point_3& q) const + { + typename Kernel::Angle_3 compute_angle = Kernel().angle_3_object(); + return compute_angle(p,angle_vertex_point,q); + } + + template + CGAL::Sign compare_distance(const Point_3& p, + const T1& obj1, + const T2& obj2) const + { + typename Kernel::Compare_distance_3 compare_distance = + Kernel().compare_distance_3_object(); + return compare_distance(p,obj1,obj2); + } + +public: + Data points_; + +private: + mutable FT length_ = -1.; + +}; // end class Polyline + + +} // end namespace internal +} // end namespace Mesh_3 +} // end namespace CGAL + +#endif // CGAL_MESH_3_INTERNAL_POLYLINE_H diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index 348e765198a..4d595db09f8 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -47,530 +48,6 @@ namespace CGAL { namespace Mesh_3 { namespace internal { -template -class Polyline -{ - typedef typename Kernel::Point_3 Point_3; - typedef typename Kernel::Segment_3 Segment_3; - typedef typename Kernel::FT FT; - - typedef std::vector Data; - -public: - typedef typename Data::const_iterator const_iterator; - typedef std::pair Point_and_location; - - Polyline() {} - ~Polyline() {} - - /// adds a point at the end of the polyline - void add_point(const Point_3& p) - { - if( points_.empty() || p != end_point() ) { - points_.push_back(p); - } - } - - /// returns the starting point of the polyline - const Point_3& start_point() const - { - CGAL_assertion( ! points_.empty() ); - return points_.front(); - } - - /// returns the ending point of the polyline - const Point_3& end_point() const - { - CGAL_assertion( ! points_.empty() ); - return points_.back(); - } - - /// returns `true` if the polyline is not degenerate - bool is_valid() const - { - return points_.size() > 1; - } - - /// returns `true` if polyline is a loop - bool is_loop() const - { - return start_point() == end_point(); - } - - const_iterator next(const_iterator it, Orientation orientation) const { - if(orientation == POSITIVE) { - CGAL_assertion(it != (points_.end() - 1)); - if(it == (points_.end() - 2)) { - CGAL_assertion(is_loop()); - it = points_.begin(); - } else { - ++it; - } - } else { - CGAL_assertion(orientation == NEGATIVE); - CGAL_assertion(it != points_.begin()); - if(it == (points_.begin() + 1)) { - CGAL_assertion(is_loop()); - it = points_.end() - 1; - } else { - --it; - } - } - return it; - } - - bool is_curve_segment_covered(CGAL::Orientation orientation, - const Point_3& c1, const Point_3& c2, - const FT sq_r1, const FT sq_r2, - const const_iterator cc1_it, - const const_iterator cc2_it) const - { - CGAL_assertion(orientation != CGAL::ZERO); - typename Kernel::Has_on_bounded_side_3 cover_pred = - Kernel().has_on_bounded_side_3_object(); - - typedef typename Kernel::Sphere_3 Sphere_3; - const Sphere_3 s1(c1, sq_r1); - const Sphere_3 s2(c2, sq_r2); - - const_iterator c1_it = cc1_it; - const_iterator c2_it = cc2_it; - - if(orientation == CGAL::NEGATIVE) { - ++c1_it; - ++c2_it; - CGAL_assertion(c1_it != points_.end()); - CGAL_assertion(c2_it != points_.end()); - } - - if(c1_it == c2_it) return cover_pred(s1, s2, c1, c2); - const_iterator next_it = this->next(c1_it, orientation); - - if(!cover_pred(s1, s2, c1, *next_it)) return false; - - for(const_iterator it = next_it; it != c2_it; /* in body */) { - next_it = this->next(it, orientation); - if(!cover_pred(s1, s2, *it, *next_it)) return false; - it = next_it; - } // end loop ]c1_it, c2_it[ - - return cover_pred(s1, s2, *c2_it, c2); - } - - bool is_curve_segment_covered(CGAL::Orientation orientation, - const Point_3& c1, const Point_3& c2, - const FT sq_r1, const FT sq_r2) const - { - return is_curve_segment_covered(orientation, - c1, c2, sq_r1, sq_r2, locate(c1), locate(c2)); - } - - FT curve_segment_length(const Point_3& p, const Point_3 q, - CGAL::Orientation orientation) const - { - CGAL_assertion(orientation != CGAL::ZERO); - const_iterator p_it = locate(p); - const_iterator q_it = locate(q); - return curve_segment_length(p, q, orientation, p_it, q_it); - } - - FT curve_segment_length(const Point_3& p, const Point_3 q, - CGAL::Orientation orientation, - const_iterator p_it, - const_iterator q_it) const - { - CGAL_assertion(orientation != CGAL::ZERO); - CGAL_assertion(p_it == locate(p)); - CGAL_assertion(q_it == locate(q)); - - if(p_it == q_it) { - const CGAL::Comparison_result cmp = compare_distance(*p_it,p,q); - if( (cmp != LARGER && orientation == POSITIVE) || - (cmp != SMALLER && orientation == NEGATIVE) ) - { - // If the orientation of `p` and `q` on the segment is compatible - // with `orientation`, then return the distance between the two - // points. - return distance(p, q); - } - } - - if(orientation == CGAL::NEGATIVE) { - ++p_it; - ++q_it; - CGAL_assertion(p_it != points_.end()); - CGAL_assertion(q_it != points_.end()); - } - - const_iterator next_it = this->next(p_it, orientation); - FT result = distance(p, *next_it); - for(const_iterator it = next_it; it != q_it; /* in body */) { - next_it = this->next(it, orientation); - result += distance(*it, *next_it); - it = next_it; - } // end loop ]p_it, q_it[ - result += distance(*q_it, q); - return result; - } - - - /// returns the angle at the first point. - /// \pre The polyline must be a loop. - Angle angle_at_first_point() const { - CGAL_precondition(is_loop()); - const Point_3& first = points_.front(); - const Point_3& next_p = points_[1]; - const Point_3& prev = points_[points_.size() - 2]; - return angle(prev, first, next_p); - } - - /// returns the length of the polyline - FT length() const - { - if(length_ < 0.) - { - FT result(0); - const_iterator it = points_.begin(); - const_iterator previous = it++; - - for(const_iterator end = points_.end(); it != end; ++it, ++previous) { - result += distance(*previous, *it); - } - length_ = result; - } - return length_; - } - - /// returns the signed geodesic distance between `p` and `q`. - FT signed_geodesic_distance(const Point_3& p, const Point_3& q) const - { - // Locate p & q on polyline - const_iterator pit = locate(p); - const_iterator qit = locate(q,false); - - return signed_geodesic_distance(p, q, pit, qit); - } - - FT signed_geodesic_distance(const Point_3& p, const Point_3& q, - const_iterator pit, const_iterator qit) const - { - CGAL_assertion(pit == locate(p)); - CGAL_assertion(qit == locate(q, false)); - - // If p and q are in the same segment of the polyline - if ( pit == qit ) - { - FT result = distance(p,q); - - // Find the closest point to *pit - if ( compare_distance(*pit,p,q) != CGAL::LARGER ) - { return result; } - else - { return -result; } - } - if(is_loop()) - { - FT positive_distance, negative_distance; - if(pit <= qit) - { - positive_distance = curve_segment_length(p, q, CGAL::POSITIVE, pit, qit); - negative_distance = length() - positive_distance; - } - else - { - negative_distance = curve_segment_length(q, p, CGAL::POSITIVE, qit, pit); - positive_distance = length() - negative_distance; - } - return (positive_distance < negative_distance) ? positive_distance : (-negative_distance); - } - else - { - return (pit <= qit) - ? curve_segment_length(p, q, CGAL::POSITIVE, pit, qit) - : ( - curve_segment_length(p, q, CGAL::NEGATIVE, pit, qit) ); - } - } - - const_iterator previous_segment_source(const_iterator it) const - { - CGAL_assertion(it != points_.end()); - if(it == first_segment_source()) - { - CGAL_assertion(is_loop()); - it = last_segment_source(); - } - else - { - --it; - } - return it; - } - - const_iterator next_segment_source(const_iterator it) const - { - CGAL_assertion(it != points_.end()); - if(it == last_segment_source()) - { - if(is_loop()) - return first_segment_source(); - else - return points_.end(); - } - else - { - ++it; - return it; - } - } - - /// returns a point at geodesic distance `distance` from p along the - /// polyline. The polyline is oriented from starting point to end point. - /// The distance could be negative. - Point_3 point_at(const Point_3& start_pt, FT distance) const - { - return point_at(start_pt, distance, locate(start_pt)).first; - } - - /// returns a point at geodesic distance `distance` from `start_pt` along the - /// polyline. The polyline is oriented from starting point to end point. - /// The distance could be negative. - Point_and_location point_at(const Point_3& start_pt, - FT distance, - const_iterator start_it) const - { - CGAL_assertion(start_it == locate(start_pt)); - - const Point_3& start_it_pt = *start_it; - const_iterator start_it_locate_pt - = (start_it == points_.begin()) ? start_it : std::prev(start_it); - - distance += curve_segment_length(start_it_pt, start_pt, CGAL::POSITIVE, - start_it_locate_pt, start_it); - - // If polyline is a loop, ensure that distance is given from start_it - if(is_loop()) - { - if(distance < FT(0)) { distance += length(); } - else if(distance > length()) { distance -= length(); } - } - else if(distance < FT(0)) // If polyline is not a loop and distance is < 0, go backward - { - Point_3 new_start = start_pt; - while(distance < FT(0)) - { - start_it = previous_segment_source(start_it); - distance += this->distance(new_start, *start_it); - new_start = *start_it; - } - } - - CGAL_assertion(distance >= FT(0)); - CGAL_assertion(distance <= length()); - - // Initialize iterators - const_iterator pit = start_it; // start at start_it, and walk forward - const_iterator previous = pit++; - - // Iterate to find which segment contains the point we want to construct - FT segment_length = this->distance(*previous, *pit); - while(distance > segment_length) - { - distance -= segment_length; - - // Increment iterators and update length - previous = next_segment_source(previous); - pit = next_segment_source(pit); - - if(pit == points_.end()) - { - CGAL_assertion(distance < this->distance(*previous, end_point())); - break; // return {*previous, previous} - } - else - segment_length = this->distance(*previous, *pit); - } - - // return point at distance from current segment source - using Vector_3 = typename Kernel::Vector_3; - auto vector = Kernel().construct_vector_3_object(); - Vector_3 v = (pit != points_.end()) ? vector(*previous, *pit) - : vector(*previous, end_point()); - - return {(*previous) + (distance / CGAL::sqrt(v.squared_length())) * v, - previous}; - } - - const_iterator locate_corner(const Point_3& p) const - { - const_iterator res = points_.end(); - if(p == start_point()) - res = points_.begin(); - else if(p == end_point()) - res = last_segment_source(); - - CGAL_assertion(res == locate(p)); - CGAL_assertion(res != points_.end()); - return res; - } - - const_iterator locate_point(const Point_3& p) const - { - return locate(p); - } - - bool are_ordered_along(const Point_3& p, const Point_3& q, - const_iterator pit, const_iterator qit) const - { - CGAL_precondition(!is_loop()); - - // Locate p & q on polyline - CGAL_assertion(pit == locate(p)); - CGAL_assertion(qit == locate(q, true)); - - // Points are not located on the same segment - if ( pit != qit ) { return (pit <= qit); } - - // pit == qit, then we have to sort p&q along (pit,pit+1) - return ( compare_distance(*pit,p,q) != CGAL::LARGER ); - } - - bool are_ordered_along(const Point_3& p, const Point_3& q) const - { - return are_ordered_along(p, q, locate(p), locate(q, true)); - } - -private: - const_iterator first_segment_source() const - { - CGAL_precondition(is_valid()); - return points_.begin(); - } - - const_iterator last_segment_source() const - { - CGAL_precondition(is_valid()); - return (points_.end() - 2); - } - - /// returns an iterator on the starting point of the segment of the - /// polyline which contains p - /// if end_point_first is true, then --end is returned instead of begin - /// if p is the starting point of a loop. - const_iterator locate(const Point_3& p, bool end_point_first=false) const - { - CGAL_precondition(is_valid()); - - // First look if p is one of the points of the polyline - const_iterator result = std::find(points_.begin(), points_.end(), p); - if ( result != points_.end() ) - { - if ( result != points_.begin() ) - { return --result; } - else - { - // Treat loops - if ( end_point_first && p == end_point() ) - { return last_segment_source(); } - else - { return result; } - } - } - - CGAL_assertion(result == points_.end()); - - // Get result by projecting p on the polyline - const_iterator it = points_.begin(); - const_iterator previous = it; - Segment_3 nearest_segment; - const_iterator nearest_vertex = it; - result = nearest_vertex; - bool nearest_is_a_segment = false; - - while ( ++it != points_.end() ) - { - Segment_3 seg (*previous, *it); - - if(nearest_is_a_segment) - { - if(compare_distance(p, *it, nearest_segment) == CGAL::SMALLER) - { - nearest_vertex = it; - nearest_is_a_segment = false; - result = it; - if (possibly(angle(*previous, *it, p) == CGAL::ACUTE) && - compare_distance(p, seg, *nearest_vertex) == CGAL::SMALLER) - { - nearest_segment = seg; - nearest_is_a_segment = true; - result = previous; - } - } - else if(compare_distance(p, seg, nearest_segment) == CGAL::SMALLER) - { - nearest_segment = seg; - result = previous; - } - } - else { - if(compare_distance(p, *it, *nearest_vertex) == CGAL::SMALLER) - { - nearest_vertex = it; - result = it; - } - if ((nearest_vertex != it || - possibly(angle(*previous, *it, p) == CGAL::ACUTE)) && - compare_distance(p, seg, *nearest_vertex) == CGAL::SMALLER) - { - nearest_segment = seg; - nearest_is_a_segment = true; - result = previous; - } - } - previous = it; - } // end the while loop on the vertices of the polyline - - if(result == points_.begin()) { - return (end_point_first && !nearest_is_a_segment) ? last_segment_source() : points_.begin(); - } else { - return result; - } - } - - - FT distance(const Point_3& p, const Point_3& q) const - { - typename Kernel::Compute_squared_distance_3 sq_distance = - Kernel().compute_squared_distance_3_object(); - return CGAL::sqrt(sq_distance(p, q)); - } - - Angle angle(const Point_3& p, - const Point_3& angle_vertex_point, - const Point_3& q) const - { - typename Kernel::Angle_3 compute_angle = Kernel().angle_3_object(); - return compute_angle(p,angle_vertex_point,q); - } - - template - CGAL::Sign compare_distance(const Point_3& p, - const T1& obj1, - const T2& obj2) const - { - typename Kernel::Compare_distance_3 compare_distance = - Kernel().compare_distance_3_object(); - return compare_distance(p,obj1,obj2); - } - -public: - Data points_; - -private: - mutable FT length_ = -1.; - -}; // end class Polyline - - template struct Mesh_domain_segment_of_curve_primitive{ typedef typename std::iterator_traits::value_type Map_value_type; diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h index 33f1ddc2c3e..5bfa8821168 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h @@ -45,6 +45,7 @@ #include #include #include +#include #include