1 // Copyright (c) 2009 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h $
7 // $Id: Triangulation_helpers.h 53d4c9b 2019-10-28T11:29:08+01:00 Mael Rouxel-Labbé
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Stephane Tayeb
12 //
13 //******************************************************************************
14 // File Description :
15 //******************************************************************************
16
17 #ifndef CGAL_MESH_3_TRIANGULATION_HELPERS_H
18 #define CGAL_MESH_3_TRIANGULATION_HELPERS_H
19
20 #include <CGAL/license/Mesh_3.h>
21
22 #include <CGAL/enum.h>
23 #include <CGAL/internal/Has_nested_type_Bare_point.h>
24 #include <CGAL/Time_stamper.h>
25
26 #include <boost/mpl/if.hpp>
27 #include <boost/mpl/identity.hpp>
28 #include <boost/type_traits/is_same.hpp>
29 #include <boost/unordered_set.hpp>
30 #include <boost/utility/enable_if.hpp>
31
32 #include <algorithm>
33 #include <iostream>
34 #include <iterator>
35 #include <limits>
36 #include <utility>
37 #include <vector>
38
39 namespace CGAL {
40
41 namespace Mesh_3 {
42
43 template<typename Tr>
44 class Triangulation_helpers
45 {
46 typedef typename Tr::Geom_traits Gt;
47
48 typedef typename Gt::FT FT;
49 typedef typename Gt::Vector_3 Vector_3;
50
51 // If `Tr` is not a triangulation that has defined Bare_point,
52 // use Point_3 as defined in the traits class.
53 typedef typename boost::mpl::eval_if_c<
54 CGAL::internal::Has_nested_type_Bare_point<Tr>::value,
55 typename CGAL::internal::Bare_point_type<Tr>,
56 boost::mpl::identity<typename Gt::Point_3>
57 >::type Bare_point;
58
59 // 'Point' is either a bare point or a weighted point, depending on the triangulation.
60 // Since 'Triangulation_helpers' can be templated by an unweighted triangulation,
61 // this is one of the rare cases where we use 'Point'.
62 typedef typename Tr::Point Point;
63
64 typedef typename Tr::Vertex Vertex;
65 typedef typename Tr::Vertex_handle Vertex_handle;
66 typedef typename Tr::Cell Cell;
67 typedef typename Tr::Cell_handle Cell_handle;
68 typedef typename Tr::Cell_iterator Cell_iterator;
69
70 typedef std::vector<Cell_handle> Cell_vector;
71
72 /**
73 * A functor to reset visited flag of each facet of a cell
74 */
75 struct Reset_facet_visited
76 {
operatorReset_facet_visited77 void operator()(const Cell_handle& c) const
78 {
79 for ( int i=0; i<4 ; ++i )
80 c->reset_visited(i);
81 }
82 };
83
84 /**
85 * A functor to get the point of a vertex vh, but replacing
86 * it by m_p when vh == m_vh
87 */
88 struct Point_getter
89 {
90 /// When the requested will be about vh, the returned point will be p
91 /// instead of vh->point()
Point_getterPoint_getter92 Point_getter(const Vertex_handle &vh, const Point&p)
93 : m_vh(vh), m_p(p)
94 {}
95
operatorPoint_getter96 const Point& operator()(const Vertex_handle &vh) const
97 {
98 return (vh == m_vh ? m_p : vh->point());
99 }
100
101 private:
102 const Vertex_handle m_vh;
103 const Point &m_p;
104 };
105
106 public:
107 /// Constructor / Destructor
Triangulation_helpers()108 Triangulation_helpers() {}
~Triangulation_helpers()109 ~Triangulation_helpers() {}
110
111 /**
112 * Returns true if moving \c v to \c p makes no topological
113 * change in \c tr
114 */
115 bool no_topological_change(Tr& tr,
116 const Vertex_handle v,
117 const Vector_3& move,
118 const Point& p,
119 Cell_vector& cells_tos) const;
120 bool no_topological_change__without_set_point(
121 const Tr& tr,
122 const Vertex_handle v,
123 const Point& p,
124 Cell_vector& cells_tos) const;
125
126 bool no_topological_change(Tr& tr,
127 const Vertex_handle v,
128 const Vector_3& move,
129 const Point& p) const;
130 bool no_topological_change__without_set_point(
131 const Tr& tr,
132 const Vertex_handle v,
133 const Point& p) const;
134
135 bool inside_protecting_balls(const Tr& tr,
136 const Vertex_handle v,
137 const Bare_point& p) const;
138
139 /**
140 * Returns the squared distance from \c vh to its closest vertex
141 *
142 * \pre `vh` is not the infinite vertex
143 */
144 template<typename Tag> // Two versions to distinguish using 'Has_visited_for_vertex_extractor'
145 FT get_sq_distance_to_closest_vertex(const Tr& tr,
146 const Vertex_handle& vh,
147 const Cell_vector& incident_cells,
148 typename boost::enable_if_c<Tag::value>::type* = nullptr) const;
149
150 // @todo are the two versions really worth it, I can't tell the difference from a time POV...
151 template<typename Tag>
152 FT get_sq_distance_to_closest_vertex(const Tr& tr,
153 const Vertex_handle& vh,
154 const Cell_vector& incident_cells,
155 typename boost::disable_if_c<Tag::value>::type* = nullptr) const;
156
157 private:
158 /**
159 * Returns true if \c v is well_oriented on each cell of \c cell_tos
160 */
161 // For sequential version
162 bool well_oriented(const Tr& tr,
163 const Cell_vector& cell_tos) const;
164 // For parallel version
165 bool well_oriented(const Tr& tr,
166 const Cell_vector& cell_tos,
167 const Point_getter& pg) const;
168 };
169
170 template<typename Tr>
171 bool
172 Triangulation_helpers<Tr>::
no_topological_change(Tr & tr,const Vertex_handle v0,const Vector_3 & move,const Point & p,Cell_vector & cells_tos)173 no_topological_change(Tr& tr,
174 const Vertex_handle v0,
175 const Vector_3& move,
176 const Point& p,
177 Cell_vector& cells_tos) const
178 {
179 if(tr.point(v0) == p)
180 return true;
181
182 // For periodic triangulations, calling this function actually takes longer than
183 // just removing and re-inserting the point directly because the periodic mesh
184 // triangulation's side_of_power_sphere() function is somewhat brute-forcy:
185 // it has to check for 27 copies of the query point to determine correctly the side.
186 // Thus, we simply return 'false' if the triangulation is a periodic triangulation.
187 //
188 // Note that the function was nevertheless adapted to work with periodic triangulation
189 // so this hack can be disabled if one day 'side_of_power_sphere()' is improved.
190 if(boost::is_same<typename Tr::Periodic_tag, Tag_true>::value)
191 return false;
192
193 typename Gt::Construct_opposite_vector_3 cov =
194 tr.geom_traits().construct_opposite_vector_3_object();
195
196 bool np = true;
197 const Point fp = tr.point(v0);
198
199 // move the point
200 tr.set_point(v0, move, p);
201
202 if(!well_oriented(tr, cells_tos))
203 {
204 // Reset (restore) v0
205 tr.set_point(v0, cov(move), fp);
206 return false;
207 }
208
209 // Reset visited tags of facets
210 std::for_each(cells_tos.begin(), cells_tos.end(), Reset_facet_visited());
211
212 for ( typename Cell_vector::iterator cit = cells_tos.begin() ;
213 cit != cells_tos.end() ;
214 ++cit )
215 {
216 Cell_handle c = *cit;
217 for(int j=0; j<4; j++)
218 {
219 // Treat each facet only once
220 if(c->is_facet_visited(j))
221 continue;
222
223 // Set facet and its mirrored version as visited
224 Cell_handle cj = c->neighbor(j);
225 int mj = tr.mirror_index(c, j);
226 c->set_facet_visited(j);
227 cj->set_facet_visited(mj);
228
229 if(tr.is_infinite(c->vertex(j)))
230 {
231 if(tr.side_of_power_sphere(c, tr.point(cj->vertex(mj)), false)
232 != CGAL::ON_UNBOUNDED_SIDE)
233 {
234 np = false;
235 break;
236 }
237 }
238 else
239 {
240 if(tr.side_of_power_sphere(cj, tr.point(c->vertex(j)), false)
241 != CGAL::ON_UNBOUNDED_SIDE)
242 {
243 np = false;
244 break;
245 }
246 }
247 }
248 }
249
250 // Reset (restore) v0
251 tr.set_point(v0, cov(move), fp);
252
253 return np;
254 }
255
256 template<typename Tr>
257 bool
258 Triangulation_helpers<Tr>::
no_topological_change__without_set_point(const Tr & tr,const Vertex_handle v0,const Point & p,Cell_vector & cells_tos)259 no_topological_change__without_set_point(
260 const Tr& tr,
261 const Vertex_handle v0,
262 const Point& p,
263 Cell_vector& cells_tos) const
264 {
265 bool np = true;
266
267 Point_getter pg(v0, p);
268
269 if(!well_oriented(tr, cells_tos, pg))
270 {
271 return false;
272 }
273
274 // Reset visited tags of facets
275 std::for_each(cells_tos.begin(), cells_tos.end(), Reset_facet_visited());
276
277 for ( typename Cell_vector::iterator cit = cells_tos.begin() ;
278 cit != cells_tos.end() ;
279 ++cit )
280 {
281 Cell_handle c = *cit;
282 for(int j=0; j<4; j++)
283 {
284 // Treat each facet only once
285 if(c->is_facet_visited(j))
286 continue;
287
288 // Set facet and it's mirror's one visited
289 Cell_handle cj = c->neighbor(j);
290 int mj = tr.mirror_index(c, j);
291 c->set_facet_visited(j);
292 cj->set_facet_visited(mj);
293
294 Vertex_handle v1 = c->vertex(j);
295 typedef typename Tr::Triangulation_data_structure TDS;
296 typedef typename TDS::Cell_range Cell_range;
297 typedef typename TDS::Vertex_range Vertex_range;
298 if(tr.is_infinite(v1))
299 {
300 // Build a copy of c, and replace V0 by a temporary vertex (position "p")
301 typename Cell_handle::value_type c_copy (*c);
302 Cell_range::Time_stamper_impl::initialize_time_stamp(&c_copy);
303 int i_v0;
304 typename Vertex_handle::value_type v;
305 Vertex_range::Time_stamper_impl::initialize_time_stamp(&v);
306 if (c_copy.has_vertex(v0, i_v0))
307 {
308 v.set_point(p);
309 c_copy.set_vertex(i_v0, Vertex_range::s_iterator_to(v));
310 }
311
312 Cell_handle c_copy_h = Cell_range::s_iterator_to(c_copy);
313 if(tr.side_of_power_sphere(c_copy_h, pg(cj->vertex(mj)), false)
314 != CGAL::ON_UNBOUNDED_SIDE)
315 {
316 np = false;
317 break;
318 }
319 }
320 else
321 {
322 // Build a copy of cj, and replace V0 by a temporary vertex (position "p")
323 typename Cell_handle::value_type cj_copy (*cj);
324 Cell_range::Time_stamper_impl::initialize_time_stamp(&cj_copy);
325 int i_v0;
326 typename Vertex_handle::value_type v;
327 Vertex_range::Time_stamper_impl::initialize_time_stamp(&v);
328 if (cj_copy.has_vertex(v0, i_v0))
329 {
330 v.set_point(p);
331 cj_copy.set_vertex(i_v0, Vertex_range::s_iterator_to(v));
332 }
333
334 Cell_handle cj_copy_h = Cell_range::s_iterator_to(cj_copy);
335 if(tr.side_of_power_sphere(cj_copy_h, pg(v1), false)
336 != CGAL::ON_UNBOUNDED_SIDE)
337 {
338 np = false;
339 break;
340 }
341 }
342 }
343 }
344
345 return np;
346 }
347
348
349 template<typename Tr>
350 bool
351 Triangulation_helpers<Tr>::
no_topological_change(Tr & tr,const Vertex_handle v0,const Vector_3 & move,const Point & p)352 no_topological_change(Tr& tr,
353 const Vertex_handle v0,
354 const Vector_3& move,
355 const Point& p) const
356 {
357 Cell_vector cells_tos;
358 cells_tos.reserve(64);
359 tr.incident_cells(v0, std::back_inserter(cells_tos));
360 return no_topological_change(tr, v0, move, p, cells_tos);
361 }
362
363 template<typename Tr>
364 bool
365 Triangulation_helpers<Tr>::
no_topological_change__without_set_point(const Tr & tr,const Vertex_handle v0,const Point & p)366 no_topological_change__without_set_point(
367 const Tr& tr,
368 const Vertex_handle v0,
369 const Point& p) const
370 {
371 Cell_vector cells_tos;
372 cells_tos.reserve(64);
373 tr.incident_cells(v0, std::back_inserter(cells_tos));
374 return no_topological_change__without_set_point(tr, v0, p, cells_tos);
375 }
376
377
378 template<typename Tr>
379 bool
380 Triangulation_helpers<Tr>::
inside_protecting_balls(const Tr & tr,const Vertex_handle v,const Bare_point & p)381 inside_protecting_balls(const Tr& tr,
382 const Vertex_handle v,
383 const Bare_point& p) const
384 {
385 typename Gt::Compare_weighted_squared_radius_3 cwsr =
386 tr.geom_traits().compare_weighted_squared_radius_3_object();
387
388 Vertex_handle nv = tr.nearest_power_vertex(p, v->cell());
389 const Point& nvwp = tr.point(nv);
390 if(cwsr(nvwp, FT(0)) == CGAL::SMALLER)
391 {
392 typename Tr::Geom_traits::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
393 const Point& nvwp = tr.point(nv);
394 // 'true' if the distance between 'p' and 'nv' is smaller or equal than the weight of 'nv'
395 return (cwsr(nvwp , - tr.min_squared_distance(p, cp(nvwp))) != CGAL::LARGER);
396 }
397
398 return false;
399 }
400
401 /// Return the squared distance from vh to its closest vertex
402 /// if `Has_visited_for_vertex_extractor` is `true`
403 template<typename Tr>
404 template<typename Tag>
405 typename Triangulation_helpers<Tr>::FT
406 Triangulation_helpers<Tr>::
get_sq_distance_to_closest_vertex(const Tr & tr,const Vertex_handle & vh,const Cell_vector & incident_cells,typename boost::enable_if_c<Tag::value>::type *)407 get_sq_distance_to_closest_vertex(const Tr& tr,
408 const Vertex_handle& vh,
409 const Cell_vector& incident_cells,
410 typename boost::enable_if_c<Tag::value>::type*) const
411 {
412 CGAL_precondition(!tr.is_infinite(vh));
413
414 typedef std::vector<Vertex_handle> Vertex_container;
415
416 // There is no need to use tr.min_squared_distance() here because we are computing
417 // distances between 'v' and a neighbor within their common cell, which means
418 // that even if we are using a periodic triangulation, the distance is correctly computed.
419 typename Gt::Compute_squared_distance_3 csqd = tr.geom_traits().compute_squared_distance_3_object();
420 typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
421
422 Vertex_container treated_vertices;
423 FT min_sq_dist = std::numeric_limits<FT>::infinity();
424
425 for(typename Cell_vector::const_iterator cit = incident_cells.begin();
426 cit != incident_cells.end(); ++cit)
427 {
428 const Cell_handle c = (*cit);
429 const int k = (*cit)->index(vh);
430 const Point& wpvh = tr.point(c, k);
431
432 // For each vertex of the cell
433 for(int i=1; i<4; ++i)
434 {
435 const int n = (k+i)&3;
436 const Vertex_handle& vn = c->vertex(n);
437
438 if(vn == Vertex_handle() ||
439 tr.is_infinite(vn) ||
440 vn->visited_for_vertex_extractor)
441 continue;
442
443 vn->visited_for_vertex_extractor = true;
444 treated_vertices.push_back(vn);
445
446 const Point& wpvn = tr.point(c, n);
447 const FT sq_d = csqd(cp(wpvh), cp(wpvn));
448
449 if(sq_d < min_sq_dist)
450 min_sq_dist = sq_d;
451 }
452 }
453
454 for(std::size_t i=0; i < treated_vertices.size(); ++i)
455 treated_vertices[i]->visited_for_vertex_extractor = false;
456
457 return min_sq_dist;
458 }
459
460 /// Return the squared distance from vh to its closest vertex
461 /// if `Has_visited_for_vertex_extractor` is `false`
462 template<typename Tr>
463 template<typename Tag>
464 typename Triangulation_helpers<Tr>::FT
465 Triangulation_helpers<Tr>::
get_sq_distance_to_closest_vertex(const Tr & tr,const Vertex_handle & vh,const Cell_vector & incident_cells,typename boost::disable_if_c<Tag::value>::type *)466 get_sq_distance_to_closest_vertex(const Tr& tr,
467 const Vertex_handle& vh,
468 const Cell_vector& incident_cells,
469 typename boost::disable_if_c<Tag::value>::type*) const
470 {
471 CGAL_precondition(!tr.is_infinite(vh));
472
473 typedef CGAL::Hash_handles_with_or_without_timestamps Hash_fct;
474 typedef boost::unordered_set<Vertex_handle, Hash_fct> Vertex_container;
475 typedef typename Vertex_container::iterator VC_it;
476
477 // There is no need to use tr.min_squared_distance() here because we are computing
478 // distances between 'v' and a neighbor within their common cell, which means
479 // that even if we are using a periodic triangulation, the distance is correctly computed.
480 typename Gt::Compute_squared_distance_3 csqd = tr.geom_traits().compute_squared_distance_3_object();
481 typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
482
483 Vertex_container treated_vertices;
484 FT min_sq_dist = std::numeric_limits<FT>::infinity();
485
486 for(typename Cell_vector::const_iterator cit = incident_cells.begin();
487 cit != incident_cells.end(); ++cit)
488 {
489 const Cell_handle c = (*cit);
490 const int k = (*cit)->index(vh);
491 const Point& wpvh = tr.point(c, k);
492
493 // For each vertex of the cell
494 for(int i=1; i<4; ++i)
495 {
496 const int n = (k+i)&3;
497 const Vertex_handle& vn = c->vertex(n);
498
499 if(vn == Vertex_handle() ||
500 tr.is_infinite(vn))
501 continue;
502
503 std::pair<VC_it, bool> is_insert_successful = treated_vertices.insert(vn);
504 if(! is_insert_successful.second) // vertex has already been treated
505 continue;
506
507 const Point& wpvn = tr.point(c, n);
508 const FT sq_d = csqd(cp(wpvh), cp(wpvn));
509
510 if(sq_d < min_sq_dist)
511 min_sq_dist = sq_d;
512 }
513 }
514
515 return min_sq_dist;
516 }
517
518
519 /// This function well_oriented is called by no_topological_change after the
520 /// position of the vertex has been (tentatively) modified.
521 template<typename Tr>
522 bool
523 Triangulation_helpers<Tr>::
well_oriented(const Tr & tr,const Cell_vector & cells_tos)524 well_oriented(const Tr& tr,
525 const Cell_vector& cells_tos) const
526 {
527 typedef typename Tr::Geom_traits Gt;
528 typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object();
529 typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
530
531 typename Cell_vector::const_iterator it = cells_tos.begin();
532 for( ; it != cells_tos.end() ; ++it)
533 {
534 Cell_handle c = *it;
535 if( tr.is_infinite(c) )
536 {
537 int iv = c->index(tr.infinite_vertex());
538 Cell_handle cj = c->neighbor(iv);
539 int mj = tr.mirror_index(c, iv);
540
541 const Point& mjwp = tr.point(cj, mj);
542 const Point& fwp1 = tr.point(c, (iv+1)&3);
543 const Point& fwp2 = tr.point(c, (iv+2)&3);
544 const Point& fwp3 = tr.point(c, (iv+3)&3);
545
546 if(orientation(cp(mjwp), cp(fwp1), cp(fwp2), cp(fwp3)) != CGAL::NEGATIVE)
547 return false;
548 }
549 else
550 {
551 const Point& cwp0 = tr.point(c, 0);
552 const Point& cwp1 = tr.point(c, 1);
553 const Point& cwp2 = tr.point(c, 2);
554 const Point& cwp3 = tr.point(c, 3);
555
556 if(orientation(cp(cwp0), cp(cwp1), cp(cwp2), cp(cwp3)) != CGAL::POSITIVE)
557 return false;
558 }
559 }
560 return true;
561 }
562
563 /// Another version for the parallel version
564 /// Here, the set_point is not done before, but we use a Point_getter instance
565 /// to get the point of a vertex.
566 template<typename Tr>
567 bool
568 Triangulation_helpers<Tr>::
well_oriented(const Tr & tr,const Cell_vector & cells_tos,const Point_getter & pg)569 well_oriented(const Tr& tr,
570 const Cell_vector& cells_tos,
571 const Point_getter& pg) const
572 {
573 typedef typename Tr::Geom_traits Gt;
574 typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object();
575 typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
576
577 typename Cell_vector::const_iterator it = cells_tos.begin();
578 for( ; it != cells_tos.end() ; ++it)
579 {
580 Cell_handle c = *it;
581 if( tr.is_infinite(c) )
582 {
583 int iv = c->index(tr.infinite_vertex());
584 Cell_handle cj = c->neighbor(iv);
585 int mj = tr.mirror_index(c, iv);
586 if(orientation(cp(pg(cj->vertex(mj))),
587 cp(pg(c->vertex((iv+1)&3))),
588 cp(pg(c->vertex((iv+2)&3))),
589 cp(pg(c->vertex((iv+3)&3)))) != CGAL::NEGATIVE)
590 return false;
591 }
592 else if(orientation(cp(pg(c->vertex(0))),
593 cp(pg(c->vertex(1))),
594 cp(pg(c->vertex(2))),
595 cp(pg(c->vertex(3)))) != CGAL::POSITIVE)
596 return false;
597 }
598 return true;
599 }
600
601
602
603 } // end namespace Mesh_3
604
605 } //namespace CGAL
606
607 #endif // CGAL_MESH_3_TRIANGULATION_HELPERS_H
608