songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge 5b574bfaae into dbcace69e6

This commit is contained in:
Jane Tournois 2025-12-02 12:16:55 +01:00 committed by GitHub
parent dbcace69e6 5b574bfaae
commit 71dcddb981
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
5 changed files with 949 additions and 487 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
View File

@ -46,6 +46,7 @@
#include <CGAL/iterator.h>
#include <CGAL/number_utils.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Mesh_3/internal/Polyline.h>
#include <CGAL/boost/iterator/transform_iterator.hpp>
@ -142,6 +143,8 @@ public:
using Distance_Function =
typename CGAL::Default::Get<DistanceFunction, NoDistanceFunction>::type;
using Polyline_iterator = typename CGAL::Mesh_3::internal::Polyline<GT>::const_iterator;
private:
typedef typename CGAL::Kernel_traits<MeshDomain>::Kernel Kernel;
typedef Delaunay_triangulation_3<Kernel> Dt;
@ -491,6 +494,8 @@ private:
return use_edge_distance_;
}
private:
C3T3& c3t3_;
const MeshDomain& domain_;
@ -598,7 +603,9 @@ Protect_edges_sizing_field<C3T3, MD, Sf, Df>::
insert_corners()
{
// Iterate on domain corners
typedef std::vector< std::pair<Corner_index, Bare_point> > Initial_corners;
using Initial_corner = std::pair<Corner_index, Bare_point>;
using Initial_corners = std::vector<Initial_corner>;
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<Bare_point, Index> point_s_;
std::pair<Bare_point, Index> 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<Vertex_handle, ErasedVeOutIt> 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));
}
}

561
Mesh_3/include/CGAL/Mesh_3/internal/Polyline.h Normal file
View File

@ -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 <CGAL/license/Mesh_3.h>
#include <CGAL/iterator.h>
#include <vector>
namespace CGAL {
/// @cond DEVELOPERS
namespace Mesh_3 {
namespace internal {
template <typename Kernel>
class Polyline
{
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Segment_3 Segment_3;
typedef typename Kernel::FT FT;
typedef std::vector<Point_3> Data;
public:
typedef typename Data::const_iterator const_iterator;
typedef std::pair<Point_3, const_iterator> 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 <typename T1, typename T2>
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

658
Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h
View File

@ -30,12 +30,14 @@
#include <CGAL/Real_timer.h>
#include <CGAL/property_map.h>
#include <CGAL/SMDS_3/internal/indices_management.h>
#include <CGAL/Mesh_3/internal/Polyline.h>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <type_traits>
#include <unordered_map>
#include <variant>
#include <memory>
@ -46,411 +48,6 @@ namespace CGAL {
namespace Mesh_3 {
namespace internal {
template <typename Kernel>
class Polyline
{
typedef typename Kernel::Point_3 Point_3;
typedef typename Kernel::Segment_3 Segment_3;
typedef typename Kernel::FT FT;
typedef std::vector<Point_3> 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 <typename T1, typename T2>
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 <typename GT, typename MapIterator>
struct Mesh_domain_segment_of_curve_primitive{
typedef typename std::iterator_traits<MapIterator>::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<GT>;
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 <typename Surf_p_index, typename IncidenceMap>
void reindex_patches(const std::vector<Surf_p_index>& map, IncidenceMap& incidence_map);
@ -837,8 +487,6 @@ private:
private:
typedef std::map<Point_3,Corner_index> Corners;
typedef Mesh_3::internal::Polyline<GT> Polyline;
typedef std::map<Curve_index, Polyline> Edges;
typedef std::map<Curve_index, Surface_patch_index_set > Edges_incidences;
typedef std::map<Corner_index, std::set<Curve_index> > Corners_tmp_incidences;
@ -851,6 +499,8 @@ private:
typedef CGAL::AABB_traits_3<GT,
Curves_primitives> 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> curves_aabb_tree_ptr_;
mutable bool curves_aabb_tree_is_built;
mutable std::unordered_map<Point_3, Polyline_const_iterator> 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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::FT
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::FT
@ -1014,6 +721,22 @@ curve_length(const Curve_index& curve_index) const
}
template <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::Point_and_location
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::Point_3
Mesh_domain_with_polyline_features_3<MD_>::
@ -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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::FT
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::FT
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
CGAL::Sign
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
CGAL::Sign
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
CGAL::Sign
@ -1517,6 +1277,34 @@ distance_sign_along_loop(const Point_3& p,
else { return CGAL::NEGATIVE; }
}
template <class MD_>
CGAL::Sign Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
bool
Mesh_domain_with_polyline_features_3<MD_>::
@ -1535,13 +1323,53 @@ Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
bool
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::Polyline_const_iterator
Mesh_domain_with_polyline_features_3<MD_>::
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 <class MD_>
typename Mesh_domain_with_polyline_features_3<MD_>::Polyline_const_iterator
Mesh_domain_with_polyline_features_3<MD_>::
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

30
Mesh_3/include/CGAL/Sizing_field_with_aabb_tree.h
View File

@ -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<Point_3>(&*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 << ") ";

33
Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h
View File

@ -45,6 +45,7 @@
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Time_stamper.h>
#include <CGAL/Mesh_3/internal/Polyline.h>
#include <CGAL/boost/iterator/transform_iterator.hpp>
@ -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<Gt>::const_iterator;
private:
typedef typename CGAL::Kernel_traits<MeshDomain>::Kernel Kernel;
@ -1884,35 +1887,40 @@ Protect_edges_sizing_field<C3T3, MD, Sf>::
insert_balls_on_edges()
{
// Get features
typedef std::tuple<Curve_index,
std::pair<Bare_point,Index>,
std::pair<Bare_point,Index> > Feature_tuple;
typedef std::vector<Feature_tuple> Input_features;
// Get features
struct Feature_tuple
{
Curve_index curve_index_;
Polyline_iterator polyline_begin_;
std::pair<Bare_point, Index> point_s_;
std::pair<Bare_point, Index> point_t_;
};
typedef std::vector<Feature_tuple> 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'.