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..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,6 +46,7 @@ #include #include #include +#include #include @@ -142,6 +143,8 @@ public: using Distance_Function = typename CGAL::Default::Get::type; + using Polyline_iterator = typename CGAL::Mesh_3::internal::Polyline::const_iterator; + private: typedef typename CGAL::Kernel_traits::Kernel Kernel; typedef Delaunay_triangulation_3 Dt; @@ -491,6 +494,8 @@ private: return use_edge_distance_; } + + private: C3T3& c3t3_; const MeshDomain& domain_; @@ -598,7 +603,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)); @@ -609,20 +616,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; @@ -671,7 +674,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); } @@ -1009,6 +1012,7 @@ insert_balls_on_edges() struct Feature_tuple { Curve_index curve_index_; + Polyline_iterator polyline_begin_; std::pair point_s_; std::pair point_t_; }; @@ -1026,16 +1030,17 @@ 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; + 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); - 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 { @@ -1054,10 +1059,11 @@ insert_balls_on_edges() FT curve_length = domain_.curve_length(curve_index); - Bare_point other_point = + auto [other_point, _] = 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, @@ -1066,7 +1072,9 @@ insert_balls_on_edges() p_index, Vertex_handle(), 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'. vq = vp; @@ -1238,10 +1246,13 @@ insert_balls(const Vertex_handle& vp, curve_index, d_sign) << ")\n"; #endif - const Bare_point new_point = - 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); + d_signF * d / 2, + 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); const FT point_weight = CGAL::square(size_(new_point, dim, index)); @@ -1256,8 +1267,10 @@ insert_balls(const Vertex_handle& vp, index, Vertex_handle(), out); - if(forced_stop()) return out; const Vertex_handle new_vertex = pair.first; + domain_.set_polyline_iterator(new_point, polyline_iter); + + if(forced_stop()) return out; out = pair.second; const FT sn = get_radius(new_vertex); if(sp <= sn) { @@ -1274,7 +1287,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 @@ -1283,9 +1296,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; @@ -1328,26 +1341,32 @@ 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 = 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; + + for(int i = 1; i <= n; ++i) { // New point position - Bare_point new_point = - domain_.construct_point_on_curve(p, curve_index, pt_dist); + 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 = smart_insert_point(new_point, point_weight, dim, index, prev, out); Vertex_handle new_vertex = pair.first; out = pair.second; + domain_.set_polyline_iterator(new_point, polyline_iter); // Add edge to c3t3 if(!c3t3_.is_in_complex(prev, new_vertex)) { @@ -1359,8 +1378,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 @@ -1553,15 +1576,28 @@ 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(); + + const Curve_index curve_index = c3t3_.curve_index(e); + 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 Curve_index curve_index = c3t3_.curve_index(e); - const FT signed_geodesic_distance = domain_.signed_geodesic_distance(pa, pb, curve_index); - const Bare_point geodesic_middle = (signed_geodesic_distance >= FT(0)) - ? domain_.construct_point_on_curve(pa, curve_index, signed_geodesic_distance / 2) - : domain_.construct_point_on_curve(pb, curve_index, -signed_geodesic_distance / 2); + const 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, + pa_it) + : domain_.construct_point_on_curve(pb, + curve_index, + -signed_geodesic_distance / 2, + pb_it); const Bare_point edge_middle = CGAL::midpoint(pa, pb); const FT squared_evaluated_distance = CGAL::squared_distance(edge_middle, geodesic_middle); @@ -1822,12 +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), + ? 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 @@ -1879,10 +1917,16 @@ 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)); + 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) { std::cerr << " But the curve is locally covered\n"; @@ -1915,18 +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); + 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_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 d31ac841756..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,12 +30,14 @@ #include #include #include +#include #include #include #include #include #include +#include #include #include @@ -46,411 +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; - - 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 - { - 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 = locate(c1); - const_iterator c2_it = locate(c2); - - 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); - } - - 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); - - 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 - { - //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 ) - { - result += distance(*previous,*it); - } - - return result; - } - - /// 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); - - // 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()) { - 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 { - return (pit <= qit) - ? curve_segment_length(p, q, CGAL::POSITIVE, pit, qit) - : ( - curve_segment_length(p, q, CGAL::NEGATIVE, pit, qit) ); - } - } - - - /// 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 - { - // use first point of the polyline instead of p - distance += curve_segment_length(start_point(),p,CGAL::POSITIVE); - - // 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 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; - ++pit; - - if (pit == points_.end()) - return *previous; - - 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); - - return (*previous) + (distance / CGAL::sqrt(v.squared_length())) * v; - } - - bool are_ordered_along(const Point_3& p, const Point_3& q) const - { - CGAL_precondition(!is_loop()); - - // 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); } - - // pit == qit, then we have to sort p&q along (pit,pit+1) - return ( compare_distance(*pit,p,q) != CGAL::LARGER ); - } - -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 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)); - } - - 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_; -}; // end class Polyline - - template struct Mesh_domain_segment_of_curve_primitive{ typedef typename std::iterator_traits::value_type Map_value_type; @@ -567,6 +164,10 @@ public: typedef GT R; typedef typename MD::Point_3 Point_3; + 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 /// @{ @@ -713,24 +314,55 @@ 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; + /// implements `MeshDomainWithFeatures_3::construct_point_on_curve()`. + Point_and_location + construct_point_on_curve(const Point_3& starting_point, + 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, 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, + 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; + /// implements `MeshDomainWithFeatures_3::is_loop()`. bool is_loop(const Curve_index& index) const; @@ -740,6 +372,19 @@ public: 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, + const FT sq_r1, const FT sq_r2, + 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; + + 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`. @@ -795,6 +440,11 @@ public: 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, + const Curve_index& curve_index) const; + template void reindex_patches(const std::vector& map, IncidenceMap& incidence_map); @@ -837,8 +487,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; @@ -851,6 +499,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_; @@ -867,6 +517,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 @@ -914,6 +565,45 @@ public: #endif } // build_curves_aabb_tree() + Polyline_const_iterator locate_in_polyline(const Point_3& p, + const int dim, + const Curve_index& index) const + { + 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); + + 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); + 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 @@ -962,6 +652,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)}; } @@ -1000,6 +691,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 @@ -1014,6 +721,22 @@ curve_length(const Curve_index& curve_index) const } +template +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, + 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, starting_point_it); +} + template typename Mesh_domain_with_polyline_features_3::Point_3 Mesh_domain_with_polyline_features_3:: @@ -1026,7 +749,7 @@ construct_point_on_curve(const Point_3& starting_point, 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); } @@ -1183,7 +906,8 @@ add_features_and_incidences(InputIterator first, InputIterator end, 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, +signed_geodesic_distance(const Point_3& p, + const Point_3& q, const Curve_index& curve_index) const { // Get corresponding polyline @@ -1191,7 +915,23 @@ 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); +} + +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 + 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, pit, qit); } @@ -1477,7 +1217,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()); @@ -1485,12 +1227,30 @@ 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; } +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 @@ -1517,6 +1277,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:: @@ -1535,13 +1323,53 @@ 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); +} + +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:: +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); +} + +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 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..cbbbdd5f4cf 100644 --- a/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h +++ b/Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h @@ -228,23 +228,28 @@ 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) { - d_ptr->dt.insert(pair.first); + for(const auto& [corner_pt, corner_index] : d_ptr->corners_indices) + { + 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 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); - Point_3 other_point = - d_ptr->domain.construct_point_on_curve(pair.first, + auto loc = + d_ptr->domain.locate_corner(curve_index, corner_pt); + auto [other_point, _] = + d_ptr->domain.construct_point_on_curve(corner_pt, curve_index, - curve_length / 2); + curve_length / 2, + loc); d_ptr->dt.insert(other_point); } } @@ -570,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 << ") "; 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..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 @@ -98,6 +99,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 +1887,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 +1951,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'.