1 // Copyright (c) 1997-2013 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/Periodic_2_triangulation_2/include/CGAL/Periodic_2_triangulation_2.h $
7 // $Id: Periodic_2_triangulation_2.h 4e519a3 2021-05-05T13:15:37+02:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 // Author(s) : Nico Kruithof <Nico@nghk.nl>
11
12 #ifndef CGAL_PERIODIC_2_TRIANGULATION_2_H
13 #define CGAL_PERIODIC_2_TRIANGULATION_2_H
14
15 #include <CGAL/license/Periodic_2_triangulation_2.h>
16
17
18 #include <CGAL/basic.h>
19
20 #include <list>
21 #include <vector>
22 #include <map>
23 #include <algorithm>
24 #include <utility>
25 #include <iostream>
26
27 #include <CGAL/iterator.h>
28 #include <CGAL/Iterator_project.h>
29 #include <CGAL/function_objects.h>
30
31 #include <CGAL/triangulation_assertions.h>
32 #include <CGAL/Triangulation_utils_2.h>
33
34 #include <CGAL/Triangulation_data_structure_2.h>
35 #include <CGAL/Periodic_2_triangulation_face_base_2.h>
36 #include <CGAL/Periodic_2_triangulation_vertex_base_2.h>
37 #include <CGAL/Periodic_2_triangulation_iterators_2.h>
38 #include <CGAL/spatial_sort.h>
39
40 #include <boost/random/linear_congruential.hpp>
41 #include <boost/random/uniform_smallint.hpp>
42 #include <boost/random/variate_generator.hpp>
43
44 #include <CGAL/utility.h>
45 #include <array>
46
47 namespace CGAL
48 {
49
50 /// Periodic triangulation class.
51 /// Its main functionality is:
52 /// - Insertion of points
53 /// - Deletion of points
54 /// - Point location
55 template < class Gt,
56 class Tds = Triangulation_data_structure_2 <
57 Periodic_2_triangulation_vertex_base_2<Gt>,
58 Periodic_2_triangulation_face_base_2<Gt> > >
59 class Periodic_2_triangulation_2
60 : public Triangulation_cw_ccw_2
61 {
62 typedef Periodic_2_triangulation_2<Gt, Tds> Self;
63
64 public:
65 // Public types of Periodic_2_triangulation_2
66 /// The triangulation data structure type
67 typedef Tds Triangulation_data_structure;
68 /// The traits class
69 typedef Gt Geom_traits;
70
71 /// The periodic offset type
72 typedef typename Gt::Periodic_2_offset_2 Offset;
73 /// The iso rectangle type
74 typedef typename Gt::Iso_rectangle_2 Iso_rectangle;
75 /// Integer tuple to store the number of sheets in each direction of space.
76 typedef std::array<int, 2> Covering_sheets;
77
78 /// The point type
79 typedef typename Gt::Point_2 Point;
80 /// The vector type
81 typedef typename Gt::Segment_2 Segment;
82 /// The segment type
83 typedef typename Gt::Vector_2 Vector;
84 /// The triangle type
85 typedef typename Gt::Triangle_2 Triangle;
86
87 /// Represents a point-offset pair. The point in the pair lies in the original domain.
88 typedef std::pair<Point, Offset> Periodic_point;
89 /// A pair of periodic points representing a segment in the periodic domain.
90 typedef std::array<std::pair<Point, Offset>, 2> Periodic_segment;
91 /// A triple of periodic points representing a triangle in the periodic domain.
92 typedef std::array<std::pair<Point, Offset>, 3> Periodic_triangle;
93
94 /// The vertex type
95 typedef typename Tds::Vertex Vertex;
96 /// The face type
97 typedef typename Tds::Face Face;
98 /// The edge type
99 typedef typename Tds::Edge Edge;
100
101 /// Size type (an unsigned integral type)
102 typedef typename Tds::size_type size_type;
103 /// Difference type (a signed integral type)
104 typedef typename Tds::difference_type difference_type;
105
106 /// Handle to a vertex
107 typedef typename Tds::Vertex_handle Vertex_handle;
108 /// Handle to a face
109 typedef typename Tds::Face_handle Face_handle;
110
111 /// Iterator over the faces
112 typedef typename Tds::Face_iterator Face_iterator;
113 /// Iterator over the edges
114 typedef typename Tds::Edge_iterator Edge_iterator;
115 /// Iterator over the vertices
116 typedef typename Tds::Vertex_iterator Vertex_iterator;
117 /// Iterator over the vertices whose corresponding points lie in the
118 /// original domain, i.e. for each set of periodic copies the
119 /// Unique_vertex_iterator iterates over exactly one representative.
120 typedef Periodic_2_triangulation_unique_vertex_iterator_2<Self>
121 Unique_vertex_iterator;
122
123 /// \name For compatibility with the Triangulation_2 class
124 // \{
125 typedef Face_iterator Finite_faces_iterator;
126 typedef Edge_iterator Finite_edges_iterator;
127 typedef Vertex_iterator Finite_vertices_iterator;
128 typedef Face_iterator All_faces_iterator;
129 // \}
130
131 /// Circulator over all faces incident to a vertex
132 typedef typename Tds::Face_circulator Face_circulator;
133 /// Circulator over all edges incident to a vertex
134 typedef typename Tds::Edge_circulator Edge_circulator;
135 /// Circulator over all vertices incident to a vertex
136 typedef typename Tds::Vertex_circulator Vertex_circulator;
137
138 /// \name Periodic iterator types
139 //\{
140 /// Iterator over all periodic triangles
141 typedef Periodic_2_triangulation_triangle_iterator_2<Self>
142 Periodic_triangle_iterator;
143 /// Iterator over all periodic segments
144 typedef Periodic_2_triangulation_segment_iterator_2<Self>
145 Periodic_segment_iterator;
146 /// Iterator over all periodic points
147 typedef Periodic_2_triangulation_point_iterator_2<Self>
148 Periodic_point_iterator;
149 //\}
150
151 /// \name Enumeration types
152 //\{
153
154 /// Type determining how to iterate over the stored simplices in the triangulation
155 enum Iterator_type
156 {
157 STORED = 0,
158 UNIQUE, // 1
159 STORED_COVER_DOMAIN, // 2
160 UNIQUE_COVER_DOMAIN // 3
161 };//3
162
163 /// Return type of a point location query
164 enum Locate_type
165 {
166 /// The query point lies on a vertex
167 VERTEX = 0,
168 /// The query point lies on an edge
169 EDGE,
170 /// The query point lies on a face
171 FACE,
172 /// The query point lies outside the affine hull of the triangulation,
173 /// which is the case when the triangulation is empty.
174 EMPTY,
175 OUTSIDE_CONVEX_HULL, // unused, for compatibility with Alpha_shape_2
176 OUTSIDE_AFFINE_HULL // unused, for compatibility with Alpha_shape_2
177 };
178
179 /// Returns false, no infinite simplices in the periodic triangulation
180 template<class T>
181 bool is_infinite(const T&, int = 0) const { return false; }
182
183 //\}
184
185 // Auxiliary iterators for convenience
186 // do not use default template argument to please VC++
187 /// Functor that returns the point given a vertex
188 typedef Project_point<Vertex> Proj_point;
189
190 /// \name STL types
191 // \{
192 /// value_type similar to stl containers
193 typedef Point value_type; // to have a back_inserter
194 /// const_reference similar to stl containers
195 typedef const value_type& const_reference;
196 /// reference similar to stl containers
197 typedef value_type& reference;
198 // \}
199
200 /// Tag to distinguish regular triangulations from others;
201 typedef Tag_false Weighted_tag;
202
203 /// Tag to distinguish periodic triangulations from others
204 typedef Tag_true Periodic_tag;
205
206 protected:
207 // Protected types of Periodic_2_triangulation_2
208 typedef typename Gt::Orientation_2 Orientation_2;
209 typedef typename Gt::Compare_x_2 Compare_x;
210 typedef typename Gt::Compare_y_2 Compare_y;
211
212
213 typedef typename Gt::FT FT;
214 typedef std::pair<Vertex_handle, Offset> Virtual_vertex;
215 typedef std::map<Vertex_handle, Virtual_vertex> Virtual_vertex_map;
216 typedef typename Virtual_vertex_map::const_iterator Virtual_vertex_map_it;
217
218 /// Vector is contains virtual copies with offset off:
219 /// virtual copy with offset off is stored at position: i=3*off[0]+off[1]-1
220 typedef std::map<Vertex_handle, std::vector<Vertex_handle> >
221 Virtual_vertex_reverse_map;
222 typedef typename Virtual_vertex_reverse_map::const_iterator
223 Virtual_vertex_reverse_map_it;
224 typedef std::map<Vertex_handle, std::list<Vertex_handle> > Too_long_edges_map;
225 typedef typename Too_long_edges_map::const_iterator Too_long_edges_map_it;
226
227 /// \name Functors
228 // \{
229 /// Functor for symbolically perturbing points
230 class Perturbation_order
231 {
232 const Self *t;
233
234 public:
235 // Perturbation_order, public interface
Perturbation_order(const Self * tr)236 Perturbation_order(const Self *tr) :
237 t(tr)
238 {
239 }
240
operator()241 bool operator()(const Point *p, const Point *q) const
242 {
243 return t->compare_xy(*p, *q) == SMALLER;
244 }
operator()245 bool operator()(const Periodic_point *p, const Periodic_point *q) const
246 {
247 return t->compare_xy(p->first, q->first, p->second, q->second) == SMALLER;
248 }
249 };
250 // \}
251 friend class Perturbation_order;
252
253 private:
254 // Copy constructor helpers
255 class Finder;
256 void copy_multiple_covering(const Periodic_2_triangulation_2 & tr);
257
258 public:
259 // Public functions of Periodic_2_triangulation_2
260 /// \name Constructors
261 //\{
262 /// Constructor
263 Periodic_2_triangulation_2(
264 const Iso_rectangle &domain = Iso_rectangle(0, 0, 1, 1),
265 const Geom_traits &geom_traits = Geom_traits());
266
267 /// Copy constructor
268 Periodic_2_triangulation_2(const Periodic_2_triangulation_2<Gt, Tds> &tr);
269 /// Assignment
270 Periodic_2_triangulation_2 &operator=(const Periodic_2_triangulation_2 &tr);
271
272 /// Copy the triangulation
273 void copy_triangulation(const Periodic_2_triangulation_2 &tr);
274 /// Swap two triangulations
275 void swap(Periodic_2_triangulation_2 &tr);
276 /// Clear the triangulation
277 void clear();
278
279 /// Serialize the triangulation to an output stream
280 std::ostream& save(std::ostream& os) const;
281
282 /// Deserialize the triangulation from an input stream
283 std::istream& load(std::istream& is);
284
285 //\}
286
287 /// \name Access functions
288 //\{
289
290 /// Returns the geometric traits used for the predicates and constructions.
geom_traits()291 const Geom_traits& geom_traits() const
292 {
293 return _gt;
294 }
295 /// Returns the data structure storing the triangulation.
tds()296 const Triangulation_data_structure & tds() const
297 {
298 return _tds;
299 }
300 /// Returns the data structure storing the triangulation.
tds()301 Triangulation_data_structure & tds()
302 {
303 return _tds;
304 }
305 /// Returns the domain of the 1-sheeted cover.
domain()306 const Iso_rectangle & domain() const
307 {
308 return _domain;
309 }
310 /// Returns the number of copies of the 1-sheeted cover stored in each of
311 /// the principal directions.
number_of_sheets()312 Covering_sheets number_of_sheets() const
313 {
314 return _cover;
315 }
316 /// Returns the dimension of the triangulation.
dimension()317 int dimension() const
318 {
319 return _tds.dimension() == 2 ? 2 : 0;
320 }
321
322 //\}
323
324 /// \name Number of simplices
325 //\{
326 /// Returns whether the triangulation is empty.
empty()327 bool empty() const
328 {
329 return _tds.dimension() < 2;
330 }
331
332 /// Returns the number of vertices. Counts all vertices that are
333 /// representatives of the same point in the 1-cover as one vertex.
number_of_vertices()334 size_type number_of_vertices() const
335 {
336 if (is_1_cover())
337 return _tds.number_of_vertices();
338 else
339 return _tds.number_of_vertices() / 9;
340 }
341
342 /// Returns the number of edges. Counts all edges that are
343 /// representatives of the same segment in the 1-cover as one edge.
number_of_edges()344 size_type number_of_edges() const
345 {
346 if (is_1_cover())
347 return _tds.number_of_edges();
348 else
349 return _tds.number_of_edges() / 9;
350 }
351 /// Returns the number of faces. Counts all faces that are
352 /// representatives of the same triangle in the 1-cover as one face.
number_of_faces()353 size_type number_of_faces() const
354 {
355 if (is_1_cover())
356 return _tds.number_of_faces();
357 else
358 return _tds.number_of_faces() / 9;
359 }
360 /// Returns the number of vertices stored in the data structure.
number_of_stored_vertices()361 size_type number_of_stored_vertices() const
362 {
363 return _tds.number_of_vertices();
364 }
365 /// Returns the number of edges stored in the data structure.
number_of_stored_edges()366 size_type number_of_stored_edges() const
367 {
368 return _tds.number_of_edges();
369 }
370 /// Returns the number of faces stored in the data structure.
number_of_stored_faces()371 size_type number_of_stored_faces() const
372 {
373 return _tds.number_of_faces();
374 }
375 //\}
376
377
378 /// \name Methods regarding the covering
379 /// \{
380
381 /// Checks whether the triangulation is a valid simplicial complex in the one cover.
382 bool is_triangulation_in_1_sheet() const;
383
384 /// Convert a 9 sheeted cover (used for sparse triangulations) to a single sheeted cover.
385 /// \pre !is_1_cover();
386 void convert_to_1_sheeted_covering();
387 /// Convert a single sheeted cover (used for dense triangulations) to a 9 sheeted cover.
388 /// \pre is_1_cover();
389 void convert_to_9_sheeted_covering();
390 // \}
391
392 /// \name Geometric access functions
393 // \{
394
395 /// Returns the periodic point given by vertex v. If t is
396 /// represented in the 1-sheeted covering space, the offset is
397 /// always zero. Otherwise v can correspond to a periodic copy
398 /// outside domain of an input point.
periodic_point(Vertex_handle v)399 Periodic_point periodic_point(Vertex_handle v) const
400 {
401 return Periodic_point(v->point(), get_offset(v));
402 }
403
404 /// If t is represented in the 1-sheeted covering space, this
405 /// function returns the periodic point given by the i-th vertex of
406 /// face f, that is the point in the original domain and the offset
407 /// of the vertex in f. If t is represented in the 9-sheeted
408 /// covering space, this offset is possibly added to another offset
409 /// determining the periodic copy.
410 /// \pre i == {0,1,2}
periodic_point(Face_handle f,int i)411 Periodic_point periodic_point(Face_handle f, int i) const
412 {
413 return Periodic_point(f->vertex(i)->point(), get_offset(f, i));
414 }
415
416 /// Returns the periodic segment formed by the two point-offset
417 /// pairs corresponding to the two vertices of edge (f,i).
418 /// \pre i == {0,1,2}
periodic_segment(Face_handle f,int i)419 Periodic_segment periodic_segment(Face_handle f, int i) const
420 {
421 CGAL_triangulation_precondition( number_of_vertices() != 0 );
422 CGAL_triangulation_precondition( i >= 0 && i <= 2);
423 return make_array(periodic_point(f, ccw(i)),
424 periodic_point(f, cw(i)));
425 }
426
427 /// Same as the previous method for edge e.
periodic_segment(const Edge & e)428 Periodic_segment periodic_segment(const Edge &e) const
429 {
430 return periodic_segment(e.first, e.second);
431 }
432
433 /// Returns the periodic triangle formed by the three point-offset
434 /// pairs corresponding to the three vertices of facet f.
periodic_triangle(Face_handle f)435 Periodic_triangle periodic_triangle(Face_handle f) const
436 {
437 return make_array(periodic_point(f, 0),
438 periodic_point(f, 1),
439 periodic_point(f, 2));
440 }
441
442 /// Converts the Periodic_point pp (point-offset pair) to the corresponding
443 /// Point in \f$R^2\f$.
point(const Periodic_point & pp)444 Point point(const Periodic_point & pp) const
445 {
446 return construct_point(pp.first, pp.second);
447 }
point(Vertex_handle v)448 Point point(Vertex_handle v) const
449 {
450 return point(periodic_point(v));
451 }
point(Face_handle fh,int i)452 Point point(Face_handle fh, int i) const
453 {
454 return point(periodic_point(fh, i));
455 }
456
457 /// Converts the Periodic_segment ps to a Segment in \f$R^2\f$.
segment(const Periodic_segment & ps)458 Segment segment(const Periodic_segment &ps) const
459 {
460 return construct_segment(ps[0].first, ps[1].first, ps[0].second, ps[1].second);
461 }
462 /// Converts the Periodic_triangle pt to a Triagle in \f$R^2\f$.
triangle(const Periodic_triangle & pt)463 Triangle triangle(const Periodic_triangle &pt) const
464 {
465 Triangle triang = construct_triangle(pt[0].first, pt[1].first, pt[2].first,
466 pt[0].second, pt[1].second, pt[2].second);
467 return triang;
468 }
469
470 /// Constructs the segment associated with the edge (f,i), respects the offset
segment(Face_handle f,int i)471 Segment segment(Face_handle f, int i) const
472 {
473 return segment(periodic_segment(f, i));
474 }
475 /// Constructs the segment associated with the edge e, respects the offset
segment(const Edge & e)476 Segment segment(const Edge& e) const
477 {
478 return segment(periodic_segment(e));
479 }
480 /// Constructs the segment associated with the edge ec, respects the offset
segment(const Edge_circulator & ec)481 Segment segment(const Edge_circulator& ec) const
482 {
483 return segment(periodic_segment(ec->first, ec->second));
484 }
485 /// Constructs the segment associated with the edge ei, respects the offset
segment(const Edge_iterator & ei)486 Segment segment(const Edge_iterator& ei) const
487 {
488 return segment(periodic_segment(ei->first, ei->second));
489 }
490
491 /// Constructs the triangle associated with the face f, respects the offset
triangle(Face_handle f)492 Triangle triangle(Face_handle f) const
493 {
494 return triangle(periodic_triangle(f));
495 }
496 //\}
497
move_in_domain(const Point & p)498 Point move_in_domain(const Point &p)
499 {
500 typename Gt::FT x = p.x();
501 typename Gt::FT y = p.y();
502
503 while (x < _domain.xmin()) x += _domain.xmax() - _domain.xmin();
504 while (x >= _domain.xmax()) x -= _domain.xmax() - _domain.xmin();
505
506 while (y < _domain.ymin()) y += _domain.ymax() - _domain.ymin();
507 while (y >= _domain.ymax()) y -= _domain.ymax() - _domain.ymin();
508
509 return Point(x, y);
510 }
511
512 /// \name Queries on simplices
513 // \{
is_edge(Vertex_handle va,Vertex_handle vb)514 bool is_edge(Vertex_handle va, Vertex_handle vb) const
515 {
516 return _tds.is_edge(va, vb);
517 }
is_edge(Vertex_handle va,Vertex_handle vb,Face_handle & fr,int & i)518 bool is_edge(Vertex_handle va, Vertex_handle vb, Face_handle& fr, int & i) const
519 {
520 return _tds.is_edge(va, vb, fr, i);
521 }
is_face(Vertex_handle v1,Vertex_handle v2,Vertex_handle v3)522 bool is_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3) const
523 {
524 return _tds.is_face(v1, v2, v3);
525 }
is_face(Vertex_handle v1,Vertex_handle v2,Vertex_handle v3,Face_handle & fr)526 bool is_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3,
527 Face_handle &fr) const
528 {
529 return _tds.is_face(v1, v2, v3, fr);
530 }
531 // \}
532
533 /// \name Queries
534 // \{
535
536 /// Wrapper function for locate if only the requested point is given.
537 Face_handle locate(const Point &p, Face_handle start = Face_handle()) const
538 {
539 Locate_type lt;
540 int li;
541 return locate(p, Offset(), lt, li, start);
542 }
543
544 /// Wrapper function for locate if the offset is omitted.
545 Face_handle locate(const Point& p, Locate_type& lt, int& li,
546 Face_handle start = Face_handle()) const
547 {
548 return locate(p, Offset(), lt, li, start);
549 }
550
551 /// Returns the oriented side of the point p with respect to the
552 /// triangle defined by the face f
oriented_side(Face_handle f,const Point & p)553 Oriented_side oriented_side(Face_handle f, const Point& p) const
554 {
555 return oriented_side(f, p, Offset());
556 }
557
558 // \}
559
560
561 /// \name Traversal of the Triangulation
562 // \{
563 /// Iterator over all stored vertices. Starts at an arbitrary vertex.
564 /// Returns vertices_end() if t.number_of_vertices()=0.
vertices_begin()565 Vertex_iterator vertices_begin() const
566 {
567 return _tds.vertices_begin();
568 }
569 /// Past the end Vertex_iterator.
vertices_end()570 Vertex_iterator vertices_end() const
571 {
572 return _tds.vertices_end();
573 }
574 /// Iterator over all stored edges. Starts at an arbitrary edge.
575 /// Returns edges_end() if t.number_of_vertices()=0.
edges_begin()576 Edge_iterator edges_begin() const
577 {
578 return _tds.edges_begin();
579 }
580 /// Past the end Edge_iterator.
edges_end()581 Edge_iterator edges_end() const
582 {
583 return _tds.edges_end();
584 }
585 /// Iterator over all stored faces. Starts at an arbitrary face.
586 /// Returns faces_end() if t.number_of_vertices()=0.
faces_begin()587 Face_iterator faces_begin() const
588 {
589 return _tds.faces_begin();
590 }
591 /// Past the end Face_iterator.
faces_end()592 Face_iterator faces_end() const
593 {
594 return _tds.faces_end();
595 }
596
597 /// Iterator over all stored vertices. Starts at an arbitrary vertex.
598 /// Returns vertices_end() if t.number_of_vertices()=0.
finite_vertices_begin()599 Vertex_iterator finite_vertices_begin() const
600 {
601 return _tds.vertices_begin();
602 }
603 /// Past the end Vertex_iterator.
finite_vertices_end()604 Vertex_iterator finite_vertices_end() const
605 {
606 return _tds.vertices_end();
607 }
608 /// Iterator over all stored edges. Starts at an arbitrary edge.
609 /// Returns edges_end() if t.number_of_vertices()=0.
finite_edges_begin()610 Edge_iterator finite_edges_begin() const
611 {
612 return _tds.edges_begin();
613 }
614 /// Past the end Edge_iterator.
finite_edges_end()615 Edge_iterator finite_edges_end() const
616 {
617 return _tds.edges_end();
618 }
619 /// Iterator over all stored faces. Starts at an arbitrary face.
620 /// Returns faces_end() if t.number_of_vertices()=0.
finite_faces_begin()621 Face_iterator finite_faces_begin() const
622 {
623 return _tds.faces_begin();
624 }
625 /// Past the end Face_iterator.
finite_faces_end()626 Face_iterator finite_faces_end() const
627 {
628 return _tds.faces_end();
629 }
630
631 /// Iterator over all stored vertices. Starts at an arbitrary vertex.
632 /// Returns vertices_end() if t.number_of_vertices()=0.
all_vertices_begin()633 Vertex_iterator all_vertices_begin() const
634 {
635 return _tds.vertices_begin();
636 }
637 /// Past the end Vertex_iterator.
all_vertices_end()638 Vertex_iterator all_vertices_end() const
639 {
640 return _tds.vertices_end();
641 }
642 /// Iterator over all stored edges. Starts at an arbitrary edge.
643 /// Returns edges_end() if t.number_of_vertices()=0.
all_edges_begin()644 Edge_iterator all_edges_begin() const
645 {
646 return _tds.edges_begin();
647 }
648 /// Past the end Edge_iterator.
all_edges_end()649 Edge_iterator all_edges_end() const
650 {
651 return _tds.edges_end();
652 }
653 /// Iterator over all stored faces. Starts at an arbitrary face.
654 /// Returns faces_end() if t.number_of_vertices()=0.
all_faces_begin()655 Face_iterator all_faces_begin() const
656 {
657 return _tds.faces_begin();
658 }
659 /// Past the end Face_iterator.
all_faces_end()660 Face_iterator all_faces_end() const
661 {
662 return _tds.faces_end();
663 }
664
665 /// begin iterator over the non-virtual vertices
unique_vertices_begin()666 Unique_vertex_iterator unique_vertices_begin() const
667 {
668 return CGAL::filter_iterator(vertices_end(),
669 Periodic_2_triangulation_2_internal::Domain_tester<Self>(this),
670 vertices_begin());
671 }
672 /// past-the-end iterator over the non-virtual vertices
unique_vertices_end()673 Unique_vertex_iterator unique_vertices_end() const
674 {
675 return CGAL::filter_iterator(vertices_end(),
676 Periodic_2_triangulation_2_internal::Domain_tester<Self>(this));
677 }
678
679 // \}
680 /// \name Geometric iterators
681 //\{
682
683 /// Start iterator over the points
684 Periodic_point_iterator periodic_points_begin(Iterator_type it = STORED) const
685 {
686 return Periodic_point_iterator(this, it);
687 }
688 /// Past-the-end iterator over the points
689 Periodic_point_iterator periodic_points_end(Iterator_type it = STORED) const
690 {
691 return Periodic_point_iterator(this, 1, it);
692 }
693
694 /// Start iterator over the segments
695 Periodic_segment_iterator periodic_segments_begin(Iterator_type it = STORED) const
696 {
697 return Periodic_segment_iterator(this, it);
698 }
699 /// Past-the-end iterator over the segments
700 Periodic_segment_iterator periodic_segments_end(Iterator_type it = STORED) const
701 {
702 return Periodic_segment_iterator(this, 1, it);
703 }
704
705 /// Start iterator over the triangles
706 Periodic_triangle_iterator periodic_triangles_begin(Iterator_type it = STORED) const
707 {
708 return Periodic_triangle_iterator(this, it);
709 }
710 /// Past-the-end iterator over the triangles
711 Periodic_triangle_iterator periodic_triangles_end(Iterator_type it = STORED) const
712 {
713 return Periodic_triangle_iterator(this, 1, it);
714 }
715 //\}
716
717 /// \name Incident simplices
718 // \{
719
720 Face_circulator incident_faces(Vertex_handle v, Face_handle f = Face_handle()) const
721 {
722 return _tds.incident_faces(v, f);
723 }
724 Edge_circulator incident_edges(Vertex_handle v, Face_handle f = Face_handle()) const
725 {
726 return _tds.incident_edges(v, f);
727 }
728 /* Vertex_circulator incident_vertices(Vertex_handle v, Face_handle f = */
729 /* Face_handle()) const */
730 /* { */
731 /* bool DEPRECATED_USE_ADJACENT_VERTICES; */
732 /* return adjacent_vertices(v, f); */
733 /* } */
734 Vertex_circulator adjacent_vertices(Vertex_handle v, Face_handle f =
735 Face_handle()) const
736 {
737 return _tds.incident_vertices(v, f);
738 }
739 // \}
740
741 /// \name Traversal between adjacent faces
742 // \{
743
mirror_vertex(Face_handle f,int i)744 Vertex_handle mirror_vertex(Face_handle f, int i) const
745 {
746 return _tds.mirror_vertex(f, i);
747 }
mirror_index(Face_handle f,int i)748 int mirror_index(Face_handle f, int i) const
749 {
750 return _tds.mirror_index(f, i);
751 }
752 //\}
753
754
755 /// \name Modifiers
756 // \{
757
758 /// Flips the edge and all periodic copies
759 void flip(Face_handle f, int i);
760
761 /// Inserts a point in the triangulation
762 /// \param p the point to be inserted
763 /// \param start the start face for point location
764 /// \return The new vertex handle or an existing Vertex_handle if p was inserted before
765 Vertex_handle insert(const Point &p, Face_handle start = Face_handle());
766
767 /// Inserts a point in the triangulation
768 /// \pre The point has been located in the triangulation
769 Vertex_handle insert(const Point& p, Locate_type lt, Face_handle loc, int li);
770
771 /// Insert a point in the triangulation
push_back(const Point & p)772 Vertex_handle push_back(const Point& p)
773 {
774 return insert(p);
775 }
776
777 // \}
778
779 /// \name Advanced modifiers
780 //\{
781
782 /// Insert the first vertex in the triangulation and creates the 9-cover.
783 Vertex_handle insert_first(const Point& p);
784 /// Inserts p in the face f and sets the offsets of the newly created faces
785 /// Insert periodic copies in all periodic copies of the domain
786 Vertex_handle insert_in_face(const Point& p, Face_handle f);
787 /// Inserts (p,o) in the edge (f,i) and sets the offsets of the newly created faces
788 /// Insert periodic copies in all periodic copies of the domain
789 Vertex_handle insert_in_edge(const Point& p, Face_handle f, int i);
790
791 /// Remove a degree 3 vertex from a 2D triangulation
792 void remove_degree_3(Vertex_handle v);
793
794 bool remove_degree_init(Vertex_handle v, const Offset &v_o,
795 std::vector<Face_handle> &f,
796 std::vector<Vertex_handle> &w, std::vector<Offset> &offset_w,
797 std::vector<int> &i, int &d, int &maxd, bool &simplicity_criterion);
798
799 /// Remove a vertex from a 2D triangulation with number_of_vertices() == 1
800 void remove_first(Vertex_handle v);
801
802 /// Changes the domain. Note that this function calls clear(), i.e.,
803 /// it erases the existing triangulation.
set_domain(const Iso_rectangle & domain)804 void set_domain(const Iso_rectangle &domain)
805 {
806 clear();
807 _domain = domain;
808 _gt.set_domain(_domain);
809 _edge_length_threshold =
810 FT(0.166) * (_domain.xmax() - _domain.xmin()) * (_domain.xmax() - _domain.xmin());
811 }
812 //\}
813
814
815 /// \name Point location
816
817 /// Do a remembering heuristic walk to locate point (p,o)
818 Face_handle
819 march_locate_2D(Face_handle f, const Point& p, const Offset& o,
820 Locate_type& lt, int& li) const;
821
822 /// Checks whether the result of two point location queries are equivalent.
823 bool
824 compare_walks(const Point& p, Face_handle c1, Face_handle c2,
825 Locate_type& lt1, Locate_type& lt2, int li1, int li2) const;
826
827 /// Testing where the point (p,off) lies w.r.t. the face f
828 Bounded_side side_of_face(const Point &p, const Offset &off, Face_handle f,
829 Locate_type <, int &li) const;
830 /// Testing where the point (p,off) lies w.r.t. the face f
side_of_face(const Point & p,Face_handle f,Locate_type & lt,int & li)831 Bounded_side side_of_face(const Point &p, Face_handle f, Locate_type <,
832 int &li) const
833 {
834 return side_of_face(p, Offset(), f, lt, li);
835 }
836 //\}
837
838 /// \name Predicates and Constructions
839 //\{
840 /// Determines whether the point p lies on the (un-)bounded side of
841 /// the triangle (p0,p1,p2)
842 Bounded_side
843 bounded_side(const Point &p0, const Point &p1, const Point &p2, const Point &p) const;
844
845 /// Determines whether the point (p,o) lies on the (un-)bounded side of
846 /// the triangle ((p0,o0),(p1,o1),(p2,o2))
847 Bounded_side
848 bounded_side(const Point &p0, const Point &p1, const Point &p2, const Point &p,
849 const Offset &o0, const Offset &o1, const Offset &o2, const Offset &o) const;
850
851 /// Determines whether the point q lies strictly between the points p and r
852 /// p,q and r are supposed to be collinear points
853 bool
854 collinear_between(const Point& p, const Point& q, const Point& r) const;
855
856 /// Determines whether the point (q,o_q) lies strictly between the points (p,o_p) and (r,o_r)
857 /// (q,o_q), (p,o_p) and (r,o_r) are supposed to be collinear points
858 bool
859 collinear_between(const Point& p, const Point& q, const Point& r,
860 const Offset& o_p, const Offset& o_q, const Offset& o_r) const;
861
862 /// Compares the x-coordinates of p and q
863 Comparison_result compare_x(const Point& p, const Point& q) const;
864 /// Compares the x-coordinates of (p,o_p) and (q,o_q)
865 Comparison_result compare_x(const Point& p, const Point& q,
866 const Offset &o_p, const Offset &o_q) const;
867
868 /// Compares p and q lexicographically
869 Comparison_result compare_xy(const Point& p, const Point& q,
870 const Offset &o_p, const Offset &o_q) const;
871 /// Compares (p,o_p) and (q,o_q) lexicographically
872 Comparison_result compare_xy(const Point& p, const Point& q) const;
873 /// Compares the x-coordinates of p and q
874 Comparison_result compare_y(const Point& p, const Point& q) const;
875 /// Compares the x-coordinates of (p,o_p) and (q,o_q)
876 Comparison_result compare_y(const Point& p, const Point& q,
877 const Offset &o_p, const Offset &o_q) const;
878 /// Checks for equality of p and q
879 bool xy_equal(const Point& p, const Point& q) const;
880 /// Returns the orientation of p1,p2,p3
881 Orientation orientation(const Point& p1, const Point& p2, const Point& p3) const;
882 /// Returns the orientation of (p1,o1), (p2,o2), (p3,o3)
883 Orientation orientation(const Point& p1, const Point& p2, const Point& p3,
884 const Offset& o1, const Offset& o2, const Offset& o3) const;
885 //\}
886
887 /// \name Miscellaneous
888 //\{
889
890 /// Returns whether the union of the faces f and f->neighbor(i) form
891 /// a convex quadrilateral.
892 bool flippable(Face_handle f, int i);
893
degree(Vertex_handle v)894 size_type degree(Vertex_handle v) const
895 {
896 return _tds.degree(v);
897 }
898
899 /// Checks if the triangulation is valid.
900 bool is_valid(bool verbose = false, int level = 0) const;
901 /// Checks if the face is valid.
902 bool is_valid(Face_handle fh, bool verbose = false, int level = 0) const;
903
904 //\}
905
906
907 /// \name Undocumented functions, needed by the geometric iterators
908 // \{
909 /// [Undoc] Returns whether the stored triangulation covers a 1-cover.
is_1_cover()910 bool is_1_cover() const
911 {
912 return (_cover[0] == 1) && (_cover[1] == 1);
913 }
914
915 /// [Undoc] Combines two offsets, where the first offset is defined by the
916 /// virtual vertex and the second by the face.
combine_offsets(const Offset & o_c,const Offset & o_t)917 Offset combine_offsets(const Offset& o_c, const Offset& o_t) const
918 {
919 Offset o_ct(_cover[0] * o_t.x(), _cover[1] * o_t.y());
920 return o_c + o_ct;
921 }
922 /// [Undoc] Returns the offset of nb==ch->neighbor(i) with respect to ch.
923 /// Get the offset between the origins of the internal offset coordinate
924 /// systems of two neighboring faces with respect from ch to nb.
925 ///
926 /// - Find two corresponding vertices from each face
927 /// - Return the difference of their offsets.
928 ///
get_neighbor_offset(Face_handle fh,int i)929 Offset get_neighbor_offset(Face_handle fh, int i) const
930 {
931 Face_handle nb = fh->neighbor(i);
932 return get_neighbor_offset(fh, i, nb, nb->index(fh));
933 }
934 /// [Undoc] Returns the offset of nb==ch->neighbor(i) with respect to ch.
935 /// Get the offset between the origins of the internal offset coordinate
936 /// systems of two neighboring faces with respect from ch to nb.
937 ///
938 /// - Find two corresponding vertices from each face
939 /// - Return the difference of their offsets.
940 ///
get_neighbor_offset(Face_handle fh,int i,Face_handle nb,int j)941 Offset get_neighbor_offset(Face_handle fh, int i, Face_handle nb, int j) const
942 {
943 // Redundance in the signature
944 CGAL_triangulation_precondition(fh->neighbor(i) == nb);
945 CGAL_triangulation_precondition(nb->neighbor(j) == fh);
946 CGAL_triangulation_precondition(fh->vertex(cw(i)) == nb->vertex(ccw(j)));
947
948 return int_to_off(nb->offset(ccw(j))) - int_to_off(fh->offset(cw(i)));
949 }
950 /// [Undoc] returns the combined offset of the vertex
951 /// (if we are not on the 1-cover) and the offset defined by the face.
get_offset(Face_handle f,int i)952 Offset get_offset(Face_handle f, int i) const
953 {
954 if (is_1_cover())
955 return int_to_off(f->offset(i));
956
957 Virtual_vertex_map_it it = _virtual_vertices.find(f->vertex(i));
958 if (it != _virtual_vertices.end())
959 return combine_offsets(it->second.second, int_to_off(f->offset(i)));
960 else
961 return combine_offsets(Offset(), int_to_off(f->offset(i)));
962 }
963 /// [Undoc] Returns the offset of the vertex if we are not on the 1-cover.
get_offset(Vertex_handle vh)964 Offset get_offset(Vertex_handle vh) const
965 {
966 if (is_1_cover())
967 return Offset();
968 Virtual_vertex_map_it it = _virtual_vertices.find(vh);
969 if (it != _virtual_vertices.end())
970 return it->second.second;
971 else
972 return Offset();
973 }
974
975 /// Converts an offset to a bit pattern where bit1==offx and bit0==offy.
off_to_int(const Offset & off)976 int off_to_int(const Offset & off) const
977 {
978 CGAL_triangulation_assertion( off.x() == 0 || off.x() == 1 );
979 CGAL_triangulation_assertion( off.y() == 0 || off.y() == 1 );
980 int i = ((off.x() & 1) << 1) + (off.y() & 1);
981 return i;
982 }
983 /// Creates an offset from a bit pattern.
int_to_off(int i)984 Offset int_to_off(int i) const
985 {
986 return Offset((i >> 1) & 1, i & 1);
987 }
988
989 // \}
990 // Protected functions of Periodic_2_triangulation_2
991 /// Const accessor to the virtual vertices reverse map,
992 /// used to optimize point location for periodic copies.
virtual_vertices_reverse()993 const Virtual_vertex_reverse_map &virtual_vertices_reverse() const
994 {
995 return _virtual_vertices_reverse;
996 }
997
998 /// [Undoc] Returns the non-virtual copy of the vertex.
get_original_vertex(Vertex_handle vh)999 Vertex_handle get_original_vertex(Vertex_handle vh) const
1000 {
1001 if (is_1_cover())
1002 return vh;
1003 Virtual_vertex_map_it it = _virtual_vertices.find(vh);
1004 if (it != _virtual_vertices.end())
1005 return it->second.first;
1006 else
1007 return vh;
1008 }
1009
1010
1011 /// Tests whether a vertex is a periodic copy of a vertex in the 3-cover.
is_virtual(Vertex_handle v)1012 bool is_virtual(Vertex_handle v)
1013 {
1014 if (is_1_cover())
1015 return false;
1016 return (_virtual_vertices.find(v) != _virtual_vertices.end());
1017 }
1018
periodic_copies(const Vertex_handle v)1019 const std::vector<Vertex_handle>& periodic_copies(const Vertex_handle v) const
1020 {
1021 CGAL_triangulation_precondition(number_of_sheets() != make_array(1, 1) );
1022 CGAL_triangulation_precondition(_virtual_vertices.find(v) == _virtual_vertices.end());
1023 CGAL_triangulation_assertion(_virtual_vertices_reverse.find(v) != _virtual_vertices_reverse.end());
1024 return _virtual_vertices_reverse.find(v)->second;
1025 }
1026
1027 template<class Stream>
draw_triangulation(Stream & os)1028 Stream& draw_triangulation(Stream& os) const
1029 {
1030 Edge_iterator it = edges_begin();
1031 for (; it != edges_end(); ++it)
1032 {
1033 os << segment(it);
1034 }
1035 return os;
1036 }
1037 protected:
1038 std::vector<Vertex_handle> insert_dummy_points();
1039
1040
try_to_convert_to_one_cover()1041 inline void try_to_convert_to_one_cover()
1042 {
1043 // Fall back to 1-cover if the criterion that the longest edge is shorter
1044 // than sqrt(0.166) is fulfilled.
1045 if (_too_long_edge_counter == 0 && !is_1_cover())
1046 {
1047 CGAL_triangulation_expensive_assertion( is_valid() );
1048 this->convert_to_1_sheeted_covering();
1049 CGAL_triangulation_expensive_assertion( is_valid() );
1050 }
1051 }
1052
1053 protected:
1054 // Protected functions of Periodic_2_triangulation_2
1055
1056 /// Inserts a point with an offset in the triangulation
1057 /// \pre The point has been located in the triangulation
1058 Vertex_handle insert(const Point& p, const Offset& o, Locate_type lt,
1059 Face_handle loc, int li, Vertex_handle vh);
1060
1061 /// \name Helper functions for queries
1062 // \{
1063
1064 /// Locates the simplex containing the point represented by p and o.
1065 ///
1066 /// The type of the simplex is stored in lt.
1067 /// The simplex containing the point is returned using lt and li.
1068 /// The Face_handle start is the start point of the heuristic walk.
1069 Face_handle
1070 locate(const Point& p, const Offset &o, Locate_type& lt, int& li,
1071 Face_handle start = Face_handle()) const;
1072 /// Returns the oriented side of the point (p,o) with respect to the
1073 /// triangle defined by the face f
1074 Oriented_side
1075 oriented_side(Face_handle f, const Point& p, const Offset &o) const;
1076 // \}
1077
1078 /// \name Insertion helpers
1079 //\{
1080 /// Insert too long edges in the star of v
1081 void insert_too_long_edges_in_star(Vertex_handle v);
1082
1083 /// Insert too long edge
1084 void insert_too_long_edge(Face_handle f, int i);
1085
1086 /// Remove too long edges in the star of v
1087 void remove_too_long_edges_in_star(Vertex_handle v);
1088
1089 /// Removes an edge if it is too long
1090 void remove_too_long_edge(Face_handle f, int i);
1091
1092 /// Check whether an edge is too long
1093 bool edge_is_too_long(const Point &p1, const Point &p2) const;
1094
1095 /// Flips the edge, no periodic copies are flipped
1096 void flip_single_edge(Face_handle f, int i);
1097
1098 /// Remove a vertex from the virtual copies maps
1099 /// Used when a Delaunay vertex is removed
1100 void remove_from_virtual_copies(Vertex_handle v);
1101 //\}
1102
1103 /// \name Wrapping the traits
1104 //\{
construct_point(const Point & p,const Offset & o)1105 Point construct_point(const Point& p, const Offset &o) const
1106 {
1107 return geom_traits().construct_point_2_object()(p, o);
1108 }
construct_point(const Periodic_point & pp)1109 Point construct_point(const Periodic_point& pp) const
1110 {
1111 return construct_point(pp.first, pp.second);
1112 }
construct_triangle(const Point & p1,const Point & p2,const Point & p3,const Offset & o1,const Offset & o2,const Offset & o3)1113 Triangle construct_triangle(const Point &p1, const Point &p2,
1114 const Point &p3, const Offset &o1, const Offset &o2, const Offset &o3) const
1115 {
1116 return geom_traits().construct_triangle_2_object()(p1, p2, p3, o1, o2, o3);
1117 }
construct_triangle(const Periodic_triangle & tri)1118 Triangle construct_triangle(const Periodic_triangle& tri) const
1119 {
1120 return construct_triangle(tri[0].first, tri[1].first, tri[2].first,
1121 tri[0].second, tri[1].second, tri[2].second);
1122 }
construct_segment(const Point & p1,const Point & p2,const Offset & o1,const Offset & o2)1123 Segment construct_segment(const Point &p1, const Point &p2, const Offset &o1,
1124 const Offset &o2) const
1125 {
1126 return geom_traits().construct_segment_2_object()(p1, p2, o1, o2);
1127 }
construct_segment(const Periodic_segment & seg)1128 Segment construct_segment(const Periodic_segment& seg) const
1129 {
1130 return construct_segment(seg[0].first, seg[1].first, seg[0].second,
1131 seg[1].second);
1132 }
1133 //\}
1134
1135 /// Test whether removing vertex v decreases the dimension of the triangulation.
test_dim_down(Vertex_handle)1136 bool test_dim_down(Vertex_handle /*v*/) const
1137 {
1138 //test the dimensionality of the resulting triangulation
1139 //upon removing of vertex v
1140 return number_of_vertices() == 1;
1141 }
1142 void make_hole(Vertex_handle v, std::list<Edge> & hole);
1143
1144
1145 Face_handle create_face(Face_handle f1, int i1, Face_handle f2, int i2, Face_handle f3, int i3);
1146 Face_handle create_face(Face_handle f1, int i1, Face_handle f2, int i2);
1147 Face_handle create_face(Face_handle f, int i, Vertex_handle v);
1148 Face_handle create_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3);
1149 Face_handle create_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, Face_handle f1, Face_handle f2, Face_handle f3);
1150 Face_handle create_face();
1151 Face_handle create_face(Face_handle); //calls copy constructor of Face
1152 void delete_face(Face_handle f);
1153 void delete_vertex(Vertex_handle v);
1154
1155 // template members
1156
well_oriented(Vertex_handle v)1157 bool well_oriented(Vertex_handle v) const
1158 {
1159 Face_circulator fc = incident_faces(v), done(fc);
1160 do
1161 {
1162 Orientation o;
1163
1164 Vertex_handle v0 = fc->vertex(0);
1165 Vertex_handle v1 = fc->vertex(1);
1166 Vertex_handle v2 = fc->vertex(2);
1167 if (fc->has_zero_offsets())
1168 {
1169 o = orientation(v0->point(), v1->point(), v2->point());
1170 }
1171 else
1172 {
1173 Offset off0 = get_offset(fc, 0);
1174 Offset off1 = get_offset(fc, 1);
1175 Offset off2 = get_offset(fc, 2);
1176 o = orientation(v0->point(), v1->point(), v2->point(), off0, off1, off2);
1177 }
1178 if (o != COUNTERCLOCKWISE) return false;
1179 }
1180 while (++fc != done);
1181 return true;
1182 }
1183
1184 template<class EdgeIt>
star_hole(const Point & p,EdgeIt edge_begin,EdgeIt edge_end)1185 Vertex_handle star_hole(const Point& p, EdgeIt edge_begin, EdgeIt edge_end)
1186 {
1187 std::list<Face_handle> empty_list;
1188 return star_hole(p, edge_begin, edge_end, empty_list.begin(),
1189 empty_list.end());
1190 }
1191
1192 template<class EdgeIt, class FaceIt>
star_hole(const Point & p,EdgeIt edge_begin,EdgeIt edge_end,FaceIt face_begin,FaceIt face_end)1193 Vertex_handle star_hole(const Point& p, EdgeIt edge_begin, EdgeIt edge_end,
1194 FaceIt face_begin, FaceIt face_end)
1195 {
1196 CGAL_assertion(is_1_cover());
1197
1198 Vertex_handle v = _tds.star_hole(edge_begin, edge_end, face_begin, face_end);
1199 v->set_point(p);
1200 return v;
1201 }
1202
1203 /// Periodic functions
1204 //\{
1205
1206 /// These functions give the pair (vertex, offset) that corresponds
1207 /// to the i-th vertex of face f. The vertex returned is not a virtual copy.
1208 void get_vertex(Face_handle f, int i, Vertex_handle &vh, Offset &off) const;
1209 /// These functions give the pair (vertex, offset) that corresponds
1210 /// to the i-th vertex of vertex vh. The vertex returned is not a virtual copy.
1211 void get_vertex(Vertex_handle vh_i, Vertex_handle &vh, Offset &off) const;
1212 /// Returns the face containing the three vertices defined by vh[0], vh1[1] and vh[2].
1213 inline Face_handle get_face(const Vertex_handle* vh) const;
1214 /// Constructs a list of too long edges in the triangulation.
1215 int find_too_long_edges(
1216 std::map<Vertex_handle, std::list<Vertex_handle> >& edges) const;
1217
1218 /// Returns the offset such that (p, o) lies on the bounded side of the face f.
get_location_offset(Face_handle f,const Point & p,const Offset & o)1219 Offset get_location_offset(Face_handle f, const Point &p, const Offset &o) const
1220 {
1221 CGAL_triangulation_precondition( number_of_vertices() != 0 );
1222
1223 if (is_1_cover() && f->has_zero_offsets())
1224 {
1225 // default case:
1226 return Offset();
1227 }
1228 else
1229 {
1230 int cumm_off = f->offset(0) | f->offset(1) | f->offset(2);
1231 // Special case for the periodic space.
1232 // Fetch vertices and respective offsets of c from _virtual_vertices
1233 const Point *pts[3];
1234 Offset off[3];
1235 for (int i = 0; i < 3; i++)
1236 {
1237 pts[i] = &(f->vertex(i)->point());
1238 off[i] = get_offset(f, i);
1239 }
1240
1241 // Main idea seems to just test all possibilities.
1242 for (int i = 0; i < 4; i++)
1243 {
1244 if (((cumm_off | (~i)) & 3) == 3)
1245 {
1246 if (bounded_side(*pts[0], *pts[1], *pts[2], p, off[0], off[1],
1247 off[2], combine_offsets(o, int_to_off(i))) != ON_UNBOUNDED_SIDE)
1248 {
1249 return int_to_off(i);
1250 }
1251 }
1252 }
1253 }
1254 CGAL_assertion(false);
1255 return Offset();
1256 }
1257
1258 /// Assigns the offsets to the vertices of the face f, and makes the offset minimal in each direction.
set_offsets(Face_handle f,int o0,int o1,int o2)1259 void set_offsets(Face_handle f, int o0, int o1, int o2)
1260 {
1261 int off0[2] = { (o0 >> 1) & 1, (o0 & 1) };
1262 int off1[2] = { (o1 >> 1) & 1, (o1 & 1) };
1263 int off2[2] = { (o2 >> 1) & 1, (o2 & 1) };
1264 // Make sure that there is at least one zero offset in every direction
1265 for (int i = 0; i < 2; i++)
1266 {
1267 int min_off = (std::min)((std::min)(off0[i], off1[i]), off2[i]);
1268 if (min_off != 0)
1269 {
1270 off0[i] -= min_off;
1271 off1[i] -= min_off;
1272 off2[i] -= min_off;
1273 }
1274 }
1275 o0 = ((off0[0] & 1) << 1) + (off0[1] & 1);
1276 o1 = ((off1[0] & 1) << 1) + (off1[1] & 1);
1277 o2 = ((off2[0] & 1) << 1) + (off2[1] & 1);
1278 f->set_offsets(o0, o1, o2);
1279 }
1280
1281 /// Assigns the offsets to the vertices of the face f, and makes the offset minimal in each direction.
1282 template<class Offset>
set_offsets(Face_handle f,const Offset & o0,const Offset & o1,const Offset & o2)1283 void set_offsets(Face_handle f, const Offset &o0, const Offset &o1, const Offset &o2)
1284 {
1285 int off0[2] = { o0.x(), o0.y() };
1286 int off1[2] = { o1.x(), o1.y() };
1287 int off2[2] = { o2.x(), o2.y() };
1288 for (int i = 0; i < 2; i++)
1289 {
1290 int min_off = (std::min)((std::min)(off0[i], off1[i]), off2[i]);
1291 if (min_off != 0)
1292 {
1293 off0[i] -= min_off;
1294 off1[i] -= min_off;
1295 off2[i] -= min_off;
1296 }
1297 }
1298
1299 CGAL_triangulation_assertion((std::min)((std::min)(off0[0], off1[0]), off2[0]) == 0);
1300 CGAL_triangulation_assertion((std::min)((std::min)(off0[1], off1[1]), off2[1]) == 0);
1301 CGAL_triangulation_assertion((0 <= off0[0]) && (off0[0] < 2));
1302 CGAL_triangulation_assertion((0 <= off1[0]) && (off1[0] < 2));
1303 CGAL_triangulation_assertion((0 <= off2[0]) && (off2[0] < 2));
1304 CGAL_triangulation_assertion((0 <= off0[1]) && (off0[1] < 2));
1305 CGAL_triangulation_assertion((0 <= off1[1]) && (off1[1] < 2));
1306 CGAL_triangulation_assertion((0 <= off2[1]) && (off2[1] < 2));
1307
1308 int o0i = ((off0[0] & 1) << 1) + (off0[1] & 1);
1309 int o1i = ((off1[0] & 1) << 1) + (off1[1] & 1);
1310 int o2i = ((off2[0] & 1) << 1) + (off2[1] & 1);
1311 f->set_offsets(o0i, o1i, o2i);
1312 }
1313 //\}
1314
1315 /// Checks the too_long_edges bookkeeping
1316 bool is_valid_too_long_edges(bool verbose = false, int level = 0) const;
1317
1318 /** @name Checking helpers */ //@{
1319 /// calls has_self_edges for every face of the triangulation
has_self_edges()1320 bool has_self_edges() const
1321 {
1322 Face_iterator it;
1323 for ( it = all_faces_begin(); it != all_faces_end(); ++it )
1324 if (has_self_edges(it)) return true;
1325 return false;
1326 }
has_self_edges(Face_handle fh)1327 bool has_self_edges(Face_handle fh) const
1328 {
1329 CGAL_triangulation_assertion((fh->vertex(0) != fh->vertex(1)) ||
1330 (fh->offset(0) != fh->offset(1)));
1331 CGAL_triangulation_assertion((fh->vertex(0) != fh->vertex(2)) ||
1332 (fh->offset(0) != fh->offset(2)));
1333 CGAL_triangulation_assertion((fh->vertex(1) != fh->vertex(2)) ||
1334 (fh->offset(1) != fh->offset(2)));
1335 return ((fh->vertex(0) == fh->vertex(1)) ||
1336 (fh->vertex(0) == fh->vertex(2)) ||
1337 (fh->vertex(1) == fh->vertex(2)));
1338 }
1339
1340 //@}
1341
1342
1343 protected:
1344 // Protected data of Periodic_2_triangulation_2
1345 /// \name Triangulation data members
1346 // \{
1347
1348 /// Geometric traits
1349 Gt _gt;
1350 /// Triangulation data structure
1351 Tds _tds;
1352 // \}
1353
1354 private:
1355 /// Inserts (p,o) in the face f and sets the offsets of the newly created faces
1356 /// Doesn't insert periodic copies
1357 Vertex_handle insert_in_face(const Point& p, const Offset &o,
1358 Face_handle f,
1359 Vertex_handle vh);
1360 /// Inserts (p,o) in the edge (f,i) and sets the offsets of the newly created faces
1361 /// Doesn't insert periodic copies
1362 Vertex_handle insert_in_edge(const Point& p, const Offset &o,
1363 Face_handle f, int i,
1364 Vertex_handle vh);
1365
1366 /// Remove a vertex without removing it's possible periodic copies.
1367 /// Helper functions
1368 void remove_degree_3_single_copy(Vertex_handle vh);
1369
1370 // Private data of Periodic_2_triangulation_2
1371 /// \name Periodic members
1372 //\{
1373
1374 /// Determines if we currently compute in 3-cover or 1-cover.
1375 Covering_sheets _cover;
1376
1377 /// The domain
1378 Iso_rectangle _domain;
1379
1380 protected:
1381 // @fixme this covering stuff should really be at the Delaunay level (will need
1382 // to be if P2RT2 is ever introduced...)
1383
1384 /// This threshold should be chosen such that if all edges are shorter,
1385 /// we can be sure that there are no self-edges anymore.
1386 FT _edge_length_threshold;
1387
1388 /// This adjacency list stores all edges that are longer than
1389 /// edge_length_threshold.
1390 Too_long_edges_map _too_long_edges;
1391 /// Number of edges that are too long
1392 size_t _too_long_edge_counter;
1393
1394 private:
1395 /// map of offsets for periodic copies of vertices
1396 Virtual_vertex_map _virtual_vertices;
1397 /// map of a non-virtual vertex to its virtual copies
1398 Virtual_vertex_reverse_map _virtual_vertices_reverse;
1399 //\}
1400 }; // class Periodic_2_triangulation_2
1401
1402 // CONSTRUCTORS
1403 template<class Gt, class Tds>
Periodic_2_triangulation_2(const Iso_rectangle & domain,const Geom_traits & geom_traits)1404 Periodic_2_triangulation_2<Gt, Tds>::Periodic_2_triangulation_2(
1405 const Iso_rectangle & domain, const Geom_traits& geom_traits)
1406 : _gt(geom_traits), _tds()
1407 , _cover(make_array(1, 1))
1408 , _domain(domain)
1409 , _too_long_edge_counter(0)
1410 {
1411 CGAL_triangulation_precondition(_domain.xmax() - _domain.xmin() ==
1412 _domain.ymax() - _domain.ymin());
1413 set_domain(_domain);
1414 }
1415
1416 // copy constructor duplicates vertices and faces
1417 template<class Gt, class Tds>
Periodic_2_triangulation_2(const Periodic_2_triangulation_2 & tr)1418 Periodic_2_triangulation_2<Gt, Tds>::Periodic_2_triangulation_2(const Periodic_2_triangulation_2 &tr)
1419 {
1420 copy_triangulation(tr);
1421 }
1422
1423 //Assignment
1424 template<class Gt, class Tds>
1425 Periodic_2_triangulation_2<Gt, Tds> &
1426 Periodic_2_triangulation_2<Gt, Tds>::operator=(
1427 const Periodic_2_triangulation_2 &tr)
1428 {
1429 copy_triangulation(tr);
1430 return *this;
1431 }
1432
1433 // Helping functions
1434 template < class GT, class Tds >
1435 class Periodic_2_triangulation_2<GT, Tds>::Finder
1436 {
1437 const Self* _t;
1438 const Point & _p;
1439 public:
Finder(const Self * t,const Point & p)1440 Finder(const Self* t, const Point &p) : _t(t), _p(p) {}
operator()1441 bool operator()(const Vertex_handle v)
1442 {
1443 return _t->xy_equal(v->point(), _p);
1444 }
1445 };
1446
1447 template < class GT, class Tds >
1448 inline void
1449 Periodic_2_triangulation_2<GT, Tds>::
copy_multiple_covering(const Periodic_2_triangulation_2<GT,Tds> & tr)1450 copy_multiple_covering(const Periodic_2_triangulation_2<GT, Tds> & tr)
1451 {
1452 // Write the respective offsets in the vertices to make them
1453 // automatically copy with the tds.
1454 for (Vertex_iterator vit = tr.vertices_begin() ;
1455 vit != tr.vertices_end() ; ++vit)
1456 {
1457 vit->set_offset(tr.get_offset(vit));
1458 }
1459 // copy the tds
1460 _tds = tr.tds();
1461 // make a list of all vertices that belong to the original
1462 // domain and initialize the basic structure of
1463 // virtual_vertices_reverse
1464 std::list<Vertex_handle> vlist;
1465 for (Vertex_iterator vit = vertices_begin() ;
1466 vit != vertices_end() ; ++vit)
1467 {
1468 if (vit->offset() == Offset())
1469 {
1470 vlist.push_back(vit);
1471 _virtual_vertices_reverse.insert(
1472 std::make_pair(vit, std::vector<Vertex_handle>(8)));
1473 CGAL_triangulation_assertion(_virtual_vertices_reverse.find(vit)
1474 ->second.size() == 8);
1475 }
1476 }
1477 // Iterate over all vertices that are not in the original domain
1478 // and construct the respective entries to virtual_vertices and
1479 // virtual_vertices_reverse
1480 for (Vertex_iterator vit2 = vertices_begin() ;
1481 vit2 != vertices_end() ; ++vit2)
1482 {
1483 if (vit2->offset() != Offset())
1484 {
1485 //TODO: use some binding, maybe boost instead of the Finder.
1486 typename std::list<Vertex_handle>::iterator vlist_it
1487 = std::find_if(vlist.begin(), vlist.end(),
1488 Finder(this, vit2->point()));
1489 Offset off = vit2->offset();
1490 _virtual_vertices.insert(std::make_pair(vit2,
1491 std::make_pair(*vlist_it, off)));
1492 _virtual_vertices_reverse.find(*vlist_it)
1493 ->second[3 * off[0] + off[1] - 1] = vit2;
1494 CGAL_triangulation_assertion(get_offset(vit2) == off);
1495 }
1496 }
1497 // Cleanup vertex offsets
1498 for (Vertex_iterator vit = vertices_begin() ;
1499 vit != vertices_end() ; ++vit)
1500 vit->clear_offset();
1501 for (Vertex_iterator vit = tr.vertices_begin() ;
1502 vit != tr.vertices_end() ; ++vit)
1503 vit->clear_offset();
1504 // Build up the too_long_edges container
1505 _too_long_edge_counter = 0;
1506 _too_long_edges.clear();
1507 for (Vertex_iterator vit = vertices_begin() ;
1508 vit != vertices_end() ; ++vit)
1509 _too_long_edges[vit] = std::list<Vertex_handle>();
1510 std::pair<Vertex_handle, Vertex_handle> edge_to_add;
1511 Point p1, p2;
1512 int i, j;
1513 for (Edge_iterator eit = edges_begin() ;
1514 eit != edges_end() ; ++eit)
1515 {
1516 if (&*(eit->first->vertex(cw(eit->second)))
1517 < &*(eit->first->vertex(ccw(eit->second))))
1518 {
1519 i = cw(eit->second);
1520 j = ccw(eit->second);
1521 }
1522 else
1523 {
1524 i = ccw(eit->second);
1525 j = cw(eit->second);
1526 }
1527 edge_to_add = std::make_pair(eit->first->vertex(i),
1528 eit->first->vertex(j));
1529 p1 = construct_point(eit->first->vertex(i)->point(),
1530 get_offset(eit->first, i));
1531 p2 = construct_point(eit->first->vertex(j)->point(),
1532 get_offset(eit->first, j));
1533 Vertex_handle v_no = eit->first->vertex(i);
1534 if (squared_distance(p1, p2) > _edge_length_threshold)
1535 {
1536 CGAL_triangulation_assertion(
1537 find(_too_long_edges[v_no].begin(),
1538 _too_long_edges[v_no].end(),
1539 edge_to_add.second) == _too_long_edges[v_no].end());
1540 _too_long_edges[v_no].push_back(edge_to_add.second);
1541 _too_long_edge_counter++;
1542 }
1543 }
1544 }
1545
1546 template<class Gt, class Tds>
copy_triangulation(const Periodic_2_triangulation_2 & tr)1547 void Periodic_2_triangulation_2<Gt, Tds>::copy_triangulation(
1548 const Periodic_2_triangulation_2 &tr)
1549 {
1550 _tds.clear();
1551 _gt = tr._gt;
1552 _cover = tr._cover;
1553 _domain = tr._domain;
1554 _edge_length_threshold = tr._edge_length_threshold;
1555 _too_long_edge_counter = tr._too_long_edge_counter;
1556 if (tr.is_1_cover())
1557 {
1558 _tds = tr.tds();
1559 }
1560 else
1561 {
1562 copy_multiple_covering(tr);
1563 }
1564 CGAL_assertion(_too_long_edge_counter == tr._too_long_edge_counter);
1565 CGAL_triangulation_expensive_postcondition(*this == tr);
1566 }
1567
1568 template<class Gt, class Tds>
swap(Periodic_2_triangulation_2 & tr)1569 void Periodic_2_triangulation_2<Gt, Tds>::swap(Periodic_2_triangulation_2 &tr)
1570 {
1571 _tds.swap(tr._tds);
1572
1573 Geom_traits t = geom_traits();
1574 _gt = tr.geom_traits();
1575 tr._gt = t;
1576
1577 std::swap(tr._cover, _cover);
1578 std::swap(tr._domain, _domain);
1579
1580 std::swap(tr._edge_length_threshold, _edge_length_threshold);
1581 std::swap(tr._too_long_edges, _too_long_edges);
1582 std::swap(tr._too_long_edge_counter, _too_long_edge_counter);
1583
1584 std::swap(tr._virtual_vertices, _virtual_vertices);
1585 std::swap(tr._virtual_vertices_reverse, _virtual_vertices_reverse);
1586 }
1587
1588 template<class Gt, class Tds>
clear()1589 void Periodic_2_triangulation_2<Gt, Tds>::clear()
1590 {
1591 _tds.clear();
1592 _tds.set_dimension(-2);
1593
1594 _too_long_edges.clear();
1595 _too_long_edge_counter = 0;
1596
1597 _virtual_vertices.clear();
1598 _virtual_vertices_reverse.clear();
1599
1600 _cover = make_array(1, 1);
1601 }
1602
1603 template<class Gt, class Tds>
is_valid(Face_handle fh,bool,int)1604 bool Periodic_2_triangulation_2<Gt, Tds>::is_valid(Face_handle fh, bool /*verbose*/, int /*level*/) const
1605 {
1606 bool result = true;
1607
1608 int xmin, xmax, ymin, ymax;
1609 xmin = ymin = 3;
1610 xmax = ymax = 0;
1611 for (int i = 0; i < 3; ++i)
1612 {
1613 Offset o = get_offset(fh, i);
1614 xmin = (std::min)(xmin, o[0]);
1615 xmax = (std::max)(xmax, o[0]);
1616 ymin = (std::min)(ymin, o[1]);
1617 ymax = (std::max)(ymax, o[1]);
1618 }
1619 // Should at most cross 1 border in each direction
1620 result &= (xmax - xmin <= 1);
1621 result &= (ymax - ymin <= 1);
1622 if (!result)
1623 {
1624 std::cerr << "min/max: " << xmin << "," << xmax << " " << ymin << "," << ymax << std::endl;
1625 for (int i = 0; i < 3; ++i)
1626 {
1627 Offset o = get_offset(fh, i);
1628 std::cerr << "Offset: " << o << std::endl;
1629 }
1630 std::cerr << std::endl;
1631 CGAL_triangulation_assertion(false);
1632 }
1633
1634 return result;
1635 }
1636
1637 template<class Gt, class Tds>
is_valid(bool verbose,int level)1638 bool Periodic_2_triangulation_2<Gt, Tds>::is_valid(bool verbose, int level) const
1639 {
1640 bool result = _tds.is_valid(verbose, level);
1641 CGAL_triangulation_assertion(result);
1642
1643 if (dimension() == 2)
1644 {
1645 // Check positive orientation:
1646 const Point *p[3];
1647 Offset off[3];
1648 for (Face_iterator fit = faces_begin(); fit != faces_end(); ++fit)
1649 {
1650 for (int i = 0; i < 3; i++)
1651 {
1652 p[i] = &fit->vertex(i)->point();
1653 off[i] = get_offset(fit, i);
1654 }
1655
1656 if (orientation(*p[0], *p[1], *p[2], off[0], off[1], off[2]) != POSITIVE)
1657 {
1658 if (verbose)
1659 {
1660 std::cerr
1661 << "Periodic_2_triangulation_2: wrong orientation:" << "\n"
1662 << *p[0] << " \t" << off[0] << "\n"
1663 << *p[1] << " \t" << off[1] << "\n"
1664 << *p[2] << " \t" << off[2] << std::endl;
1665 }
1666 result = false;
1667 }
1668 }
1669 }
1670 CGAL_triangulation_assertion(result);
1671
1672 // Check for the right number of simplices
1673 int copies = number_of_sheets()[0] * number_of_sheets()[1];
1674 result &= (number_of_stored_vertices() == copies * number_of_vertices());
1675 result &= (number_of_stored_edges() == copies * number_of_edges());
1676 result &= (number_of_stored_faces() == copies * number_of_faces());
1677 CGAL_triangulation_assertion(result);
1678
1679 // check number of euler characteristic. This cannot be done by the Tds
1680 // which does not know the genus
1681 result &= (number_of_stored_vertices() - number_of_stored_edges()
1682 + number_of_stored_faces() == 0);
1683 CGAL_triangulation_assertion(result);
1684
1685 result &= !has_self_edges();
1686 CGAL_triangulation_assertion(result);
1687
1688 // Edges should not be longer than 1 periodicity
1689 for (Face_iterator fit = faces_begin(); fit != faces_end(); ++fit)
1690 {
1691 result &= is_valid(fit, verbose, level);
1692 }
1693 CGAL_triangulation_assertion(result);
1694
1695 result &= is_1_cover() == _virtual_vertices.empty();
1696 result &= is_1_cover() == _virtual_vertices_reverse.empty();
1697 result &= (_virtual_vertices.size() == (number_of_sheets()[0]
1698 * number_of_sheets()[1] - 1) * _virtual_vertices_reverse.size());
1699 CGAL_triangulation_assertion(result);
1700
1701 for (Virtual_vertex_map_it it = _virtual_vertices.begin(); it
1702 != _virtual_vertices.end(); ++it)
1703 {
1704 const Vertex_handle © = it->first;
1705 const Vertex_handle &orig = it->second.first;
1706 const Offset &off = it->second.second;
1707 size_t index = number_of_sheets()[0] * off[0] + off[1] - 1;
1708 Virtual_vertex_reverse_map_it rev_it = _virtual_vertices_reverse.find(orig);
1709 if (rev_it != _virtual_vertices_reverse.end())
1710 {
1711 if (index < rev_it->second.size())
1712 {
1713 result &= (rev_it->second[index] == copy);
1714 }
1715 else
1716 {
1717 result &= false;
1718 }
1719 }
1720 else
1721 {
1722 result &= false;
1723 }
1724 }
1725 CGAL_triangulation_assertion(result);
1726
1727 for (Virtual_vertex_reverse_map_it it = _virtual_vertices_reverse.begin(); it
1728 != _virtual_vertices_reverse.end(); ++it)
1729 {
1730 const std::vector<Vertex_handle> &copies = it->second;
1731 result &= copies.size() == 8;
1732 for (size_t i = 0; i < copies.size(); ++i)
1733 {
1734 Virtual_vertex_map_it copy_it = _virtual_vertices.find(copies[i]);
1735 if (copy_it != _virtual_vertices.end())
1736 {
1737 result &= copy_it->second.first == it->first;
1738 }
1739 else
1740 {
1741 result &= false;
1742 }
1743 }
1744 }
1745
1746 // Check the too_long_edges administration
1747 result &= is_valid_too_long_edges(verbose, level);
1748
1749 return result;
1750 }
1751
1752 template<class Gt, class Tds>
is_valid_too_long_edges(bool verbose,int)1753 bool Periodic_2_triangulation_2<Gt, Tds>::is_valid_too_long_edges(bool verbose, int /*level*/) const
1754 {
1755 bool result = true;
1756
1757 result &= is_1_cover() == _too_long_edges.empty();
1758 CGAL_triangulation_assertion(result);
1759 size_t too_long_edges = 0;
1760 for (Too_long_edges_map_it it = _too_long_edges.begin(); it
1761 != _too_long_edges.end(); ++it)
1762 {
1763 too_long_edges += it->second.size();
1764 }
1765 CGAL_triangulation_assertion(result);
1766 if (_too_long_edge_counter != too_long_edges)
1767 {
1768 if (verbose) std::cout << "Too long edge counter is incorrect: " << _too_long_edge_counter << " != " << too_long_edges << std::endl;
1769 result = false;
1770 }
1771 CGAL_triangulation_assertion(result);
1772
1773 /// Expensive check whether the right too long edges are in the list
1774 if (is_1_cover())
1775 {
1776 for (Edge_iterator eit = edges_begin(); eit != edges_end(); ++eit)
1777 {
1778 Vertex_handle vh1 = eit->first->vertex(ccw(eit->second));
1779 Vertex_handle vh2 = eit->first->vertex(cw(eit->second));
1780 Point p1 = construct_point(vh1->point(), get_offset(eit->first, ccw(eit->second)));
1781 Point p2 = construct_point(vh2->point(), get_offset(eit->first, cw(eit->second)));
1782 result &= (!edge_is_too_long(p1, p2));
1783 }
1784 CGAL_triangulation_assertion(result);
1785 }
1786 else
1787 {
1788 too_long_edges = 0;
1789 for (Edge_iterator eit = edges_begin(); eit != edges_end(); ++eit)
1790 {
1791 Vertex_handle vh1 = eit->first->vertex(ccw(eit->second));
1792 Vertex_handle vh2 = eit->first->vertex(cw(eit->second));
1793 Point p1 = construct_point(vh1->point(),
1794 get_offset(eit->first, ccw(eit->second)));
1795 Point p2 = construct_point(vh2->point(),
1796 get_offset(eit->first, cw(eit->second)));
1797
1798 if (&*vh2 < &*vh1)
1799 std::swap(vh1, vh2);
1800 CGAL_triangulation_assertion(&*vh1 < &*vh2);
1801
1802 bool too_long = edge_is_too_long(p1, p2);
1803 if (too_long != edge_is_too_long(p2, p1))
1804 {
1805 if (verbose) std::cout << "Long edge criterion not symmetric c(v1,v2) != c(v2,v1)" << std::endl;
1806 result = false;
1807 }
1808 CGAL_triangulation_assertion(result);
1809
1810 Too_long_edges_map_it it = _too_long_edges.find(vh1);
1811 if (it == _too_long_edges.end())
1812 {
1813 if (too_long)
1814 {
1815 if (verbose) std::cout << "1. Too long edge not in the data structure" << std::endl;
1816 result = false;
1817 }
1818 result &= !too_long;
1819 CGAL_triangulation_assertion(result);
1820 }
1821 else
1822 {
1823 typename std::list<Vertex_handle>::const_iterator it2 = find(it->second.begin(), it->second.end(), vh2);
1824 if (too_long)
1825 {
1826 too_long_edges++;
1827 if (it2 == it->second.end())
1828 {
1829 if (verbose) std::cout << "2. Too long edge not in the data structure" << std::endl;
1830 result = false;
1831 }
1832 CGAL_triangulation_assertion(result);
1833 }
1834 else
1835 {
1836 if (it2 != it->second.end())
1837 {
1838 if (verbose) std::cout << "Edge is not too long, but contained in the data structure" << std::endl;
1839 result = false;
1840 }
1841 CGAL_triangulation_assertion(result);
1842 }
1843 }
1844 }
1845
1846 if (_too_long_edge_counter != too_long_edges)
1847 {
1848 if (verbose)
1849 std::cout << "Counts do not match: " << _too_long_edge_counter << " != " << too_long_edges << std::endl;
1850 result = false;
1851 }
1852 CGAL_triangulation_assertion(result);
1853 }
1854
1855 return result;
1856 }
1857
1858 template<class Gt, class Tds>
flippable(Face_handle f,int i)1859 bool Periodic_2_triangulation_2<Gt, Tds>::flippable(Face_handle f, int i)
1860 {
1861 Face_handle nb = f->neighbor(i);
1862 int j = nb->index(f);
1863
1864 const Point *p[4];
1865
1866 p[0] = &f->vertex(i)->point(); // i
1867 p[1] = &nb->vertex(j)->point(); // opposite
1868 p[2] = &f->vertex(ccw(i))->point(); // ccw
1869 p[3] = &f->vertex(cw(i))->point(); // cw
1870
1871 if (is_1_cover() && f->has_zero_offsets() && nb->has_zero_offsets())
1872 {
1873 // if (orientation(*p[0], *p[1], *p[2]) != RIGHT_TURN)
1874 // return false;
1875 // if (orientation(*p[0], *p[1], *p[3]) != LEFT_TURN)
1876 // return false;
1877 if (orientation(*p[0], *p[1], *p[2]) == LEFT_TURN)
1878 return false;
1879 if (orientation(*p[0], *p[1], *p[3]) == RIGHT_TURN)
1880 return false;
1881 }
1882 else
1883 {
1884 Offset off[4];
1885 off[0] = get_offset(f, i);
1886 off[1] = combine_offsets(get_offset(nb, j), get_neighbor_offset(nb, j, f, i));
1887 off[2] = get_offset(f, ccw(i));
1888 off[3] = get_offset(f, cw(i));
1889
1890 // if (orientation(*p[0], *p[1], *p[2], off[0], off[1], off[2]) != RIGHT_TURN)
1891 // return false;
1892 // if (orientation(*p[0], *p[1], *p[3], off[0], off[1], off[3]) != LEFT_TURN)
1893 // return false;
1894 if (orientation(*p[0], *p[1], *p[2], off[0], off[1], off[2]) == LEFT_TURN)
1895 return false;
1896 if (orientation(*p[0], *p[1], *p[3], off[0], off[1], off[3]) == RIGHT_TURN)
1897 return false;
1898 }
1899
1900 return true;
1901 }
1902
1903 template<class Gt, class Tds>
flip(Face_handle f,int i)1904 void Periodic_2_triangulation_2<Gt, Tds>::flip(Face_handle f, int i)
1905 {
1906 if (is_1_cover())
1907 {
1908 flip_single_edge(f, i);
1909 return;
1910 }
1911
1912 Vertex_handle vh1 = f->vertex( cw(i));
1913 Vertex_handle vh2 = f->vertex(ccw(i));
1914 Virtual_vertex_map_it it_vh1 = _virtual_vertices.find(vh1);
1915 Virtual_vertex_map_it it_vh2 = _virtual_vertices.find(vh2);
1916
1917 Offset vh1_offset, vh2_offset;
1918 if (it_vh1 != _virtual_vertices.end())
1919 {
1920 vh1 = it_vh1->second.first;
1921 vh1_offset = it_vh1->second.second;
1922 }
1923 if (it_vh2 != _virtual_vertices.end())
1924 {
1925 vh2 = it_vh2->second.first;
1926 vh2_offset = it_vh2->second.second;
1927 }
1928
1929 CGAL_triangulation_assertion(
1930 virtual_vertices_reverse().find(vh1) != virtual_vertices_reverse().end());
1931 CGAL_triangulation_assertion(
1932 virtual_vertices_reverse().find(vh2) != virtual_vertices_reverse().end());
1933 const std::vector<Vertex_handle> &v1s =
1934 virtual_vertices_reverse().find(vh1)->second;
1935 const std::vector<Vertex_handle> &v2s =
1936 virtual_vertices_reverse().find(vh2)->second;
1937
1938 CGAL_assertion(v1s.size() == 8);
1939 CGAL_assertion(v1s.size() == v2s.size());
1940
1941 Face_handle fh;
1942 int index=0;
1943 Vertex_handle vh1_copy, vh2_copy;
1944
1945 // Virtual copies
1946 for (int x = 0; x < 3; ++x)
1947 {
1948 for (int y = 0; y < 3; ++y)
1949 {
1950 int i1 = 3 * ((x + vh1_offset.x()) % 3) + ((y + vh1_offset.y()) % 3);
1951 int i2 = 3 * ((x + vh2_offset.x()) % 3) + ((y + vh2_offset.y()) % 3);
1952
1953 if (i1 == 0)
1954 vh1_copy = vh1;
1955 else
1956 vh1_copy = v1s[i1 - 1];
1957 if (i2 == 0)
1958 vh2_copy = vh2;
1959 else
1960 vh2_copy = v2s[i2 - 1];
1961
1962 bool found = is_edge(vh1_copy, vh2_copy, fh, index);
1963 CGAL_USE(found);
1964 CGAL_assertion(found);
1965 if (found)
1966 flip_single_edge(fh, index);
1967 }
1968 }
1969
1970 try_to_convert_to_one_cover();
1971 }
1972 template<class Gt, class Tds>
flip_single_edge(Face_handle f,int i)1973 void Periodic_2_triangulation_2<Gt, Tds>::flip_single_edge(Face_handle f, int i)
1974 {
1975 CGAL_triangulation_precondition(f != Face_handle());
1976 CGAL_triangulation_precondition(i == 0 || i == 1 || i == 2);
1977 CGAL_triangulation_precondition(dimension() == 2);
1978
1979 CGAL_triangulation_precondition(flippable(f, i));
1980
1981 if (!is_1_cover())
1982 remove_too_long_edge(f, i);
1983
1984 Face_handle nb = f->neighbor(i);
1985 if (f->has_zero_offsets() && nb->has_zero_offsets())
1986 {
1987 _tds.flip(f, i);
1988
1989 if (!is_1_cover())
1990 insert_too_long_edge(f, ccw(i));
1991
1992 return;
1993 }
1994
1995 int nb_index = nb->index(f);
1996 int offsets[4];
1997 offsets[0] = f->offset(i);
1998 offsets[1] = f->offset(cw(i));
1999 offsets[2] = f->offset(ccw(i));
2000 offsets[3] = nb->offset(nb_index);
2001
2002 // Move the offsets of f and nb in the same space by correcting for nb_offset
2003 Offset nb_offset = get_neighbor_offset(f, i, nb, nb_index);
2004 if (nb_offset.x() != 0)
2005 {
2006 if (nb_offset.x() == 1)
2007 {
2008 CGAL_assertion(((offsets[0] & 2) | (offsets[1] & 2) | (offsets[2] & 2)) == 0);
2009 offsets[0] |= 2;
2010 offsets[1] |= 2;
2011 offsets[2] |= 2;
2012 }
2013 else
2014 {
2015 CGAL_triangulation_assertion(nb_offset.x() == -1);
2016 CGAL_assertion((offsets[3] & 2) == 0);
2017 offsets[3] |= 2;
2018 }
2019 }
2020 if (nb_offset.y() != 0)
2021 {
2022 if (nb_offset.y() == 1)
2023 {
2024 CGAL_assertion(((offsets[0] & 1) | (offsets[1] & 1) | (offsets[2] & 1)) == 0);
2025 offsets[0] |= 1;
2026 offsets[1] |= 1;
2027 offsets[2] |= 1;
2028 }
2029 else
2030 {
2031 CGAL_triangulation_assertion(nb_offset.y() == -1);
2032 CGAL_assertion((offsets[3] & 1) == 0);
2033 offsets[3] |= 1;
2034 }
2035 }
2036 CGAL_assertion((offsets[0] & offsets[1] & offsets[2] & offsets[3]) == 0);
2037 CGAL_triangulation_assertion_code(Vertex_handle vh = f->vertex(i));
2038 CGAL_triangulation_assertion_code(Vertex_handle vh_ccw = f->vertex(ccw(i)));
2039 _tds.flip(f, i);
2040 // Combinatorial checks
2041 CGAL_triangulation_assertion(vh == f->vertex(i));
2042 CGAL_triangulation_assertion(vh_ccw == f->vertex(ccw(i)));
2043 CGAL_triangulation_assertion(f->vertex(i) == nb->vertex(cw(nb_index)));
2044 CGAL_triangulation_assertion(f->vertex(cw(i)) == nb->vertex(nb_index));
2045
2046 // Restore the offsets
2047 int new_off[3];
2048 // For face f
2049 new_off[i] = offsets[0];
2050 new_off[ccw(i)] = offsets[2];
2051 new_off[cw(i)] = offsets[3];
2052 set_offsets(f, new_off[0], new_off[1], new_off[2]);
2053 // For face nb
2054 new_off[nb_index] = offsets[3];
2055 new_off[ccw(nb_index)] = offsets[1];
2056 new_off[cw(nb_index)] = offsets[0];
2057 set_offsets(nb, new_off[0], new_off[1], new_off[2]);
2058
2059 if (!is_1_cover())
2060 insert_too_long_edge(f, ccw(i));
2061 }
2062
2063 template<class Gt, class Tds>
2064 void
remove_from_virtual_copies(Vertex_handle v)2065 Periodic_2_triangulation_2<Gt, Tds>::remove_from_virtual_copies(Vertex_handle v)
2066 {
2067 typename Virtual_vertex_reverse_map::iterator rev_it = _virtual_vertices_reverse.find(v);
2068 CGAL_triangulation_assertion(rev_it != _virtual_vertices_reverse.end());
2069
2070 const std::vector<Vertex_handle> &virtual_copies = rev_it->second;
2071 for (size_t i = 0; i < virtual_copies.size(); ++i)
2072 {
2073 _virtual_vertices.erase(virtual_copies[i]);
2074 }
2075 _virtual_vertices_reverse.erase(rev_it);
2076 }
2077
2078 template<class Gt, class Tds>
2079 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle Periodic_2_triangulation_2 <
insert_first(const Point & p)2080 Gt, Tds >::insert_first(const Point& p)
2081 {
2082 CGAL_assertion(empty());
2083 // The empty triangulation has a single sheeted cover
2084 _cover = make_array(3, 3);
2085
2086 /// Virtual vertices, one per periodic domain
2087 Vertex_handle vir_vertices[3][3];
2088 /// Virtual faces, two per periodic domain
2089 Face_handle faces[3][3][2];
2090
2091 // Initialise vertices:
2092 vir_vertices[0][0] = _tds.create_vertex();
2093 vir_vertices[0][0]->set_point(p);
2094 _virtual_vertices_reverse[vir_vertices[0][0]] = std::vector<Vertex_handle>();
2095 for (int i = 0; i < _cover[0]; i++)
2096 {
2097 for (int j = 0; j < _cover[1]; j++)
2098 {
2099 if ((i != 0) || (j != 0))
2100 {
2101 // Initialise virtual vertices out of the domain for debugging
2102 vir_vertices[i][j] = _tds.create_vertex();
2103 vir_vertices[i][j]->set_point(p); //+Offset(i,j));
2104 _virtual_vertices[vir_vertices[i][j]] = Virtual_vertex(
2105 vir_vertices[0][0], Offset(i, j));
2106 _virtual_vertices_reverse[vir_vertices[0][0]].push_back(
2107 vir_vertices[i][j]);
2108 }
2109 }
2110 }
2111
2112 // Create faces:
2113 for (int i = 0; i < _cover[0]; i++)
2114 {
2115 for (int j = 0; j < _cover[1]; j++)
2116 {
2117 for (int f = 0; f < 2; f++)
2118 {
2119 // f faces per 'rectangle'
2120 faces[i][j][f] = _tds.create_face();
2121 }
2122 }
2123 }
2124
2125 // table containing the vertex information
2126 // index to the right vertex: [number of faces][vertex][offset]
2127 int vertex_ind[2][3][2] = { { { 0, 0 }, { 1, 1 }, { 0, 1 } }, { { 0, 0 }, {
2128 1, 0
2129 }, { 1, 1 }
2130 }
2131 };
2132 // Table containing the neighbor information
2133 // [number of faces][neighbor][offset,face]
2134 int neighb_ind[2][3][3] = { { { 0, 1, 1 }, { -1, 0, 1 }, { 0, 0, 1 } }, { {
2135 1, 0, 0
2136 }, { 0, 0, 0 }, { 0, -1, 0 }
2137 }
2138 };
2139 for (int i = 0; i < _cover[0]; i++)
2140 {
2141 for (int j = 0; j < _cover[1]; j++)
2142 {
2143 int offset =
2144 ((i == _cover[0] - 1 ? 2 : 0) | (j == _cover[1] - 1 ? 1 : 0));
2145 for (int f = 0; f < 2; f++)
2146 {
2147 faces[i][j][f]->set_vertices(vir_vertices[(i + vertex_ind[f][0][0])
2148 % _cover[0]][(j + vertex_ind[f][0][1]) % _cover[1]],
2149 vir_vertices[(i + vertex_ind[f][1][0]) % _cover[0]][(j
2150 + vertex_ind[f][1][1]) % _cover[1]], vir_vertices[(i
2151 + vertex_ind[f][2][0]) % _cover[0]][(j + vertex_ind[f][2][1])
2152 % _cover[1]]);
2153 set_offsets(faces[i][j][f], offset & (vertex_ind[f][0][0] * 2
2154 + vertex_ind[f][0][1] * 1), offset & (vertex_ind[f][1][0] * 2
2155 + vertex_ind[f][1][1] * 1), offset & (vertex_ind[f][2][0] * 2
2156 + vertex_ind[f][2][1] * 1));
2157 faces[i][j][f]->set_neighbors(faces[(i + _cover[0]
2158 + neighb_ind[f][0][0]) % _cover[0]][(j + _cover[1]
2159 + neighb_ind[f][0][1]) % _cover[1]][neighb_ind[f][0][2]], faces[(i
2160 + _cover[0] + neighb_ind[f][1][0]) % _cover[0]][(j + _cover[1]
2161 + neighb_ind[f][1][1]) % _cover[1]][neighb_ind[f][1][2]], faces[(i
2162 + _cover[0] + neighb_ind[f][2][0]) % _cover[0]][(j + _cover[1]
2163 + neighb_ind[f][2][1]) % _cover[1]][neighb_ind[f][2][2]]);
2164 }
2165 }
2166 }
2167 // set pointers from the vertices to incident faces.
2168 for (int i = 0; i < _cover[0]; i++)
2169 {
2170 for (int j = 0; j < _cover[1]; j++)
2171 {
2172 vir_vertices[i][j]->set_face(faces[i][j][0]);
2173 }
2174 }
2175
2176 _tds.set_dimension(2);
2177
2178 // create the base for too_long_edges;
2179 CGAL_triangulation_assertion(_too_long_edges.empty() );
2180 CGAL_triangulation_assertion(_too_long_edge_counter == 0);
2181
2182 // Insert all vertices as the first vertex in the _too_long_edges list
2183 int k = 0;
2184 std::list<Vertex_handle> empty_list;
2185 for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit)
2186 {
2187 _too_long_edges[vit] = empty_list;
2188 k++;
2189 }
2190
2191 // Insert all edges as all edges are too long
2192 _too_long_edge_counter = 0;
2193 for (Edge_iterator eit = edges_begin(); eit != edges_end(); eit++)
2194 {
2195 Vertex_handle vh1 = eit->first->vertex(ccw(eit->second));
2196 Vertex_handle vh2 = eit->first->vertex(cw(eit->second));
2197 if (&*vh1 < &*vh2)
2198 {
2199 _too_long_edges[vh1].push_back(vh2);
2200 }
2201 else
2202 {
2203 _too_long_edges[vh2].push_back(vh1);
2204 }
2205 _too_long_edge_counter++;
2206 }
2207
2208 return vir_vertices[0][0];
2209 }
2210
2211 template<class Gt, class Tds>
2212 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert_in_edge(const Point & p,Face_handle f,int i)2213 Periodic_2_triangulation_2<Gt, Tds>::insert_in_edge(const Point& p,
2214 Face_handle f, int i)
2215 {
2216 return insert(p, EDGE, f, i);
2217 }
2218 template<class Gt, class Tds>
2219 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert_in_edge(const Point & p,const Offset & o,Face_handle f,int i,Vertex_handle vh)2220 Periodic_2_triangulation_2<Gt, Tds>::insert_in_edge(const Point& p, const Offset &o,
2221 Face_handle f, int i,
2222 Vertex_handle vh)
2223 {
2224 // Insert in edge calls an insert_in_face and a flip.
2225 // Therefore there is no need to update the too_long_edges bookkeeping directly.
2226
2227 CGAL_triangulation_assertion(number_of_vertices() != 0);
2228 CGAL_triangulation_assertion((!is_1_cover()) || (o == Offset()));
2229
2230 // Backup of the neighbor and its relative offset
2231 Face_handle nb = f->neighbor(i);
2232 int j = nb->index(f);
2233 CGAL_triangulation_assertion_code(Offset current_offset = get_location_offset(f, p, o));
2234 CGAL_triangulation_assertion
2235 (orientation(f->vertex(cw(i))->point(), p, f->vertex(ccw(i))->point(),
2236 get_offset(f, cw(i)), combine_offsets(o, current_offset), get_offset(f, ccw(i))) == COLLINEAR &&
2237 collinear_between(f->vertex(cw(i))->point(), p, f->vertex(ccw(i))->point(),
2238 get_offset(f, cw(i)), combine_offsets(o, current_offset), get_offset(f, ccw(i))) );
2239
2240 /// Insert in the face and flip an edge
2241 Vertex_handle v = insert_in_face(p, o, f, vh);
2242 flip_single_edge(nb, j);
2243
2244 return v;
2245 }
2246
2247 template<class Gt, class Tds>
2248 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert_in_face(const Point & p,Face_handle f)2249 Periodic_2_triangulation_2<Gt, Tds>::insert_in_face(const Point& p, Face_handle f)
2250 {
2251 return insert(p, FACE, f, 0);
2252 }
2253 template<class Gt, class Tds>
2254 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert_in_face(const Point & p,const Offset & o,Face_handle f,Vertex_handle vh)2255 Periodic_2_triangulation_2<Gt, Tds>::insert_in_face(const Point& p, const Offset &o,
2256 Face_handle f,
2257 Vertex_handle vh)
2258 {
2259 CGAL_triangulation_assertion(f != Face_handle());
2260 CGAL_triangulation_assertion(number_of_vertices() != 0);
2261 CGAL_triangulation_assertion((!is_1_cover()) || (o == Offset()));
2262
2263 const bool simplicity_criterion = f->has_zero_offsets() && o.is_zero();
2264
2265
2266 Offset current_off;
2267
2268 // Save the neighbors and the offsets
2269 Face_handle nb[3];
2270 int nb_index[3];
2271 int offsets[3];
2272 CGAL_triangulation_assertion_code( Vertex_handle vertices[3]; )
2273
2274 if (!simplicity_criterion)
2275 {
2276 // Choose the periodic copy of tester.point() that is inside c.
2277 current_off = get_location_offset(f, p, o);
2278
2279 CGAL_triangulation_assertion(oriented_side(f, p, combine_offsets(o, current_off)) != ON_NEGATIVE_SIDE);
2280
2281 for (int i = 0; i < 3; ++i)
2282 {
2283 nb[i] = f->neighbor(i);
2284 nb_index[i] = nb[i]->index(f);
2285 offsets[i] = f->offset(i);
2286 CGAL_triangulation_assertion_code( vertices[i] = f->vertex(i); );
2287 }
2288 }
2289
2290 // Insert the new vertex
2291 Vertex_handle v = _tds.insert_in_face(f);
2292 v->set_point(p);
2293
2294 if (!simplicity_criterion)
2295 {
2296 // Update the offsets
2297 int v_offset = off_to_int(current_off);
2298 int new_offsets[3];
2299 for (int i = 0; i < 3; ++i)
2300 {
2301 Face_handle new_face = nb[i]->neighbor(nb_index[i]);
2302 int v_index = new_face->index(v);
2303
2304 CGAL_triangulation_assertion(new_face->vertex(ccw(v_index)) == vertices[ccw(i)]);
2305 CGAL_triangulation_assertion(new_face->vertex(cw(v_index)) == vertices[cw(i)]);
2306
2307 new_offsets[v_index] = v_offset;
2308 new_offsets[ccw(v_index)] = offsets[ccw(i)];
2309 new_offsets[cw(v_index)] = offsets[cw(i)];
2310 set_offsets(new_face, new_offsets[0], new_offsets[1], new_offsets[2]);
2311 }
2312 }
2313
2314 if (!is_1_cover())
2315 {
2316 // update the book-keeping in case of a periodic copy
2317 if (vh != Vertex_handle())
2318 {
2319 _virtual_vertices[v] = Virtual_vertex(vh, o);
2320 _virtual_vertices_reverse[vh].push_back(v);
2321 }
2322
2323 insert_too_long_edges_in_star(v);
2324 }
2325
2326 return v;
2327 }
2328
2329 template<class Gt, class Tds>
2330 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert(const Point & p,Face_handle start)2331 Periodic_2_triangulation_2<Gt, Tds>::insert(const Point &p, Face_handle start)
2332 {
2333 CGAL_triangulation_assertion((_domain.xmin() <= p.x()) &&
2334 (p.x() < _domain.xmax()));
2335 CGAL_triangulation_assertion((_domain.ymin() <= p.y()) &&
2336 (p.y() < _domain.ymax()));
2337
2338 if (number_of_stored_vertices() == 0)
2339 {
2340 return insert_first(p);
2341 }
2342
2343 if (start == Face_handle())
2344 {
2345 start = faces_begin();
2346 }
2347
2348 Locate_type lt;
2349 int li;
2350 Face_handle loc = locate(p, lt, li, start);
2351
2352 if (start != Face_handle())
2353 {
2354 CGAL_assertion(start->vertex(0) != Vertex_handle());
2355 }
2356 return insert(p, lt, loc, li);
2357 }
2358
2359 template<class Gt, class Tds>
2360 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle
insert(const Point & p,Locate_type lt,Face_handle loc,int li)2361 Periodic_2_triangulation_2<Gt, Tds>::insert(const Point& p,
2362 Locate_type lt, Face_handle loc, int li)
2363 {
2364 if (number_of_stored_vertices() == 0)
2365 {
2366 return insert_first(p);
2367 }
2368
2369 // vstart is a vertex incident to the Face_handle start that will be used as
2370 // for creating a start point for the virtual vertices.
2371 // We use the virtual copies of a vertex incident to loc.
2372 Vertex_handle vstart;
2373 if (!is_1_cover())
2374 {
2375 Virtual_vertex_map_it vvmit = _virtual_vertices.find(loc->vertex(0));
2376 if (vvmit == _virtual_vertices.end())
2377 {
2378 vstart = loc->vertex(0);
2379 }
2380 else
2381 {
2382 vstart = vvmit->second.first;
2383 }
2384
2385 // vstart should be non-virtual, but should have virtual copies
2386 CGAL_triangulation_assertion(_virtual_vertices.find(vstart)
2387 == _virtual_vertices.end());
2388 CGAL_triangulation_assertion(_virtual_vertices_reverse.find(vstart)
2389 != _virtual_vertices_reverse.end());
2390 }
2391
2392 Vertex_handle vh = insert(p, Offset(), lt, loc, li, Vertex_handle());
2393
2394 // Don't add periodic copies if we are on the 1-cover
2395 if (is_1_cover())
2396 return vh;
2397
2398 // Don't continue if the point lies on a vertex as this will break the
2399 // start_vertices array below.
2400 if (lt == VERTEX)
2401 return vh;
2402
2403 const std::vector<Vertex_handle> &start_vertices =
2404 _virtual_vertices_reverse.find(vstart)->second;
2405 CGAL_assertion(start_vertices.size() == size_t(number_of_sheets()[0] * number_of_sheets()[1] - 1));
2406
2407 for (int i = 0; i < number_of_sheets()[0]; i++)
2408 {
2409 for (int j = 0; j < number_of_sheets()[1]; j++)
2410 {
2411 if ((i != 0) || (j != 0))
2412 {
2413 loc = start_vertices[i * 3 + j - 1]->face();
2414 Offset offset(i, j);
2415
2416 loc = locate(p, offset, lt, li, loc);
2417
2418 insert(p, offset, lt, loc, li, vh);
2419 }
2420 }
2421 }
2422
2423 return vh;
2424 }
2425
2426 template<class Gt, class Tds>
2427 typename Periodic_2_triangulation_2<Gt, Tds>::Vertex_handle Periodic_2_triangulation_2 <
insert(const Point & p,const Offset & o,Locate_type lt,Face_handle loc,int li,Vertex_handle vh)2428 Gt, Tds >::insert(const Point& p, const Offset &o, Locate_type lt,
2429 Face_handle loc, int li, Vertex_handle vh)
2430 // insert a point p, whose localization is known (lt, f, i)
2431 {
2432 Vertex_handle result;
2433 switch (lt)
2434 {
2435 case FACE:
2436 {
2437 result = insert_in_face(p, o, loc, vh);
2438 break;
2439 }
2440 case EDGE:
2441 {
2442 result = insert_in_edge(p, o, loc, li, vh);
2443 break;
2444 }
2445 case VERTEX:
2446 {
2447 // The vertex is a special case, we can return immediately
2448 CGAL_assertion(vh == Vertex_handle());
2449 return loc->vertex(li);
2450 }
2451 case EMPTY:
2452 {
2453 result = insert_first(p);
2454 break;
2455 }
2456 default:
2457 {
2458 CGAL_triangulation_assertion(false); // locate step failed
2459 return Vertex_handle();
2460 }
2461 }
2462
2463 if (!is_1_cover() && (vh == Vertex_handle()))
2464 {
2465 _virtual_vertices_reverse[result] = std::vector<Vertex_handle>();
2466 }
2467
2468 return result;
2469 }
2470
2471 template<class Gt, class Tds>
remove_degree_3(Vertex_handle v)2472 inline void Periodic_2_triangulation_2<Gt, Tds>::remove_degree_3(Vertex_handle v)
2473 {
2474 CGAL_assertion(number_of_vertices() > 1);
2475 CGAL_assertion(degree(v) == 3);
2476
2477 if (is_1_cover())
2478 {
2479 remove_degree_3_single_copy(v);
2480 return;
2481 }
2482
2483 {
2484 Virtual_vertex_map_it it = _virtual_vertices.find(v);
2485 if (it != _virtual_vertices.end())
2486 {
2487 v = it->second.first;
2488 }
2489 }
2490
2491 remove_too_long_edges_in_star(v);
2492
2493 typename Virtual_vertex_reverse_map::iterator reverse_it =
2494 _virtual_vertices_reverse.find(v);
2495 CGAL_assertion(reverse_it != _virtual_vertices_reverse.end());
2496
2497 const std::vector<Vertex_handle> &virtual_copies = reverse_it->second;
2498 for (typename std::vector<Vertex_handle>::const_iterator it = virtual_copies.begin();
2499 it != virtual_copies.end(); ++it)
2500 {
2501 _virtual_vertices.erase(*it);
2502 remove_degree_3_single_copy(*it);
2503 }
2504
2505 _virtual_vertices_reverse.erase(reverse_it);
2506 remove_degree_3_single_copy(v);
2507
2508 }
2509
2510 template<class Gt, class Tds>
remove_degree_3_single_copy(Vertex_handle vh)2511 inline void Periodic_2_triangulation_2<Gt, Tds>::remove_degree_3_single_copy(Vertex_handle vh)
2512 {
2513 Face_handle f = vh->face();
2514 int i = ccw(f->index(vh));
2515 Face_handle f2 = f->neighbor(i);
2516 int j = f2->index(f);
2517 // Get the offsets in ccw order
2518 Offset off[3];
2519 off[i] = get_offset(f, i);
2520 off[ccw(i)] = get_offset(f, ccw(i));
2521 off[cw(i)] = combine_offsets(get_offset(f2, j), get_neighbor_offset(f2, j, f, i));
2522 if (off[0].x() < 0 || off[1].x() < 0 || off[2].x() < 0)
2523 {
2524 Offset o(number_of_sheets()[0], 0);
2525 off[0] += o;
2526 off[1] += o;
2527 off[2] += o;
2528 }
2529 if (off[0].y() < 0 || off[1].y() < 0 || off[2].y() < 0)
2530 {
2531 Offset o(0, number_of_sheets()[1]);
2532 off[0] += o;
2533 off[1] += o;
2534 off[2] += o;
2535 }
2536
2537 // Remove the vertex, keep face f
2538 _tds.remove_degree_3(vh, f);
2539
2540 // Reset the offsets
2541 set_offsets(f,
2542 (off[0].x() >= number_of_sheets()[0] ? 2 : 0) + (off[0].y() >= number_of_sheets()[1] ? 1 : 0),
2543 (off[1].x() >= number_of_sheets()[0] ? 2 : 0) + (off[1].y() >= number_of_sheets()[1] ? 1 : 0),
2544 (off[2].x() >= number_of_sheets()[0] ? 2 : 0) + (off[2].y() >= number_of_sheets()[1] ? 1 : 0));
2545 }
2546
2547 template<class Gt, class Tds>
remove_first(Vertex_handle)2548 inline void Periodic_2_triangulation_2<Gt, Tds>::remove_first(Vertex_handle)
2549 {
2550 CGAL_assertion(number_of_vertices() == 1);
2551 clear();
2552 return;
2553 }
2554
2555
2556 template < class Gt, class Tds >
2557 bool
2558 Periodic_2_triangulation_2<Gt, Tds>::
remove_degree_init(Vertex_handle v,const Offset & v_o,std::vector<Face_handle> & f,std::vector<Vertex_handle> & w,std::vector<Offset> & offset_w,std::vector<int> & i,int & d,int & maxd,bool & simplicity_criterion)2559 remove_degree_init(Vertex_handle v, const Offset &v_o,
2560 std::vector<Face_handle> &f,
2561 std::vector<Vertex_handle> &w,
2562 std::vector<Offset> &offset_w,
2563 std::vector<int> &i,
2564 int &d, int &maxd,
2565 bool &simplicity_criterion)
2566 {
2567 Bbox_2 bbox = v->point().bbox();
2568 simplicity_criterion = is_1_cover();
2569
2570 f[0] = v->face();
2571 d = 0;
2572
2573 do
2574 {
2575 i[d] = f[d]->index(v);
2576 w[d] = f[d]->vertex( ccw(i[d]) );
2577 offset_w[d] = get_offset(f[d], ccw(i[d])) - get_offset(f[d], i[d]) + v_o;
2578 w[d]->set_face( f[d]->neighbor(i[d])); // do no longer bother about set_face
2579 simplicity_criterion &= (offset_w[d] == offset_w[0]);
2580
2581 bbox = bbox + this->construct_point(w[d]->point(), offset_w[d]).bbox();
2582
2583 ++d;
2584 if ( d == maxd)
2585 {
2586 maxd *= 2;
2587 f.resize(maxd);
2588 w.resize(maxd);
2589 offset_w.resize(maxd);
2590 i.resize(maxd);
2591 }
2592
2593 f[d] = f[d - 1]->neighbor( ccw(i[d - 1]) );
2594
2595 }
2596 while(f[d] != f[0]);
2597
2598 return is_1_cover() &&
2599 this->edge_is_too_long(Point(bbox.xmin(), bbox.ymin()), Point(bbox.xmax(), bbox.ymax()));
2600 }
2601
2602 template<class Gt, class Tds>
make_hole(Vertex_handle v,std::list<Edge> & hole)2603 void Periodic_2_triangulation_2<Gt, Tds>::make_hole(Vertex_handle v, std::list<Edge> & hole)
2604 {
2605 remove_too_long_edges_in_star(v);
2606
2607 std::list<Face_handle> to_delete;
2608
2609 Face_handle f, fn;
2610 int i, in;
2611 Vertex_handle vv;
2612
2613 Face_circulator fc = incident_faces(v);
2614 Face_circulator done(fc);
2615 do
2616 {
2617 f = fc;
2618 fc++;
2619 i = f->index(v);
2620 fn = f->neighbor(i);
2621 in = fn->index(f);
2622 vv = f->vertex(cw(i));
2623 if (vv->face() == f)
2624 vv->set_face(fn);
2625 vv = f->vertex(ccw(i));
2626 if (vv->face() == f)
2627 vv->set_face(fn);
2628 fn->set_neighbor(in, Face_handle());
2629 hole.push_back(Edge(fn, in));
2630 to_delete.push_back(f);
2631 }
2632 while (fc != done);
2633
2634 while (!to_delete.empty())
2635 {
2636 delete_face(to_delete.front());
2637 to_delete.pop_front();
2638 }
2639 return;
2640 }
2641
2642 template<class Gt, class Tds>
2643 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Face_handle f1,int i1,Face_handle f2,int i2,Face_handle f3,int i3)2644 Gt, Tds >::create_face(Face_handle f1, int i1, Face_handle f2, int i2,
2645 Face_handle f3, int i3)
2646 {
2647 return _tds.create_face(f1, i1, f2, i2, f3, i3);
2648 }
2649
2650 template<class Gt, class Tds>
2651 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Face_handle f1,int i1,Face_handle f2,int i2)2652 Gt, Tds >::create_face(Face_handle f1, int i1, Face_handle f2, int i2)
2653 {
2654 return _tds.create_face(f1, i1, f2, i2);
2655 }
2656
2657 template<class Gt, class Tds>
2658 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Face_handle f,int i,Vertex_handle v)2659 Gt, Tds >::create_face(Face_handle f, int i, Vertex_handle v)
2660 {
2661 return _tds.create_face(f, i, v);
2662 }
2663
2664 template<class Gt, class Tds>
2665 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Vertex_handle v1,Vertex_handle v2,Vertex_handle v3)2666 Gt, Tds >::create_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3)
2667 {
2668 return _tds.create_face(v1, v2, v3);
2669 }
2670
2671 template<class Gt, class Tds>
2672 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Vertex_handle v1,Vertex_handle v2,Vertex_handle v3,Face_handle f1,Face_handle f2,Face_handle f3)2673 Gt, Tds >::create_face(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3,
2674 Face_handle f1, Face_handle f2, Face_handle f3)
2675 {
2676 return _tds.create_face(v1, v2, v3, f1, f2, f3);
2677 }
2678
2679 template<class Gt, class Tds>
2680 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face(Face_handle fh)2681 Gt, Tds >::create_face(Face_handle fh)
2682 {
2683 return _tds.create_face(fh);
2684 }
2685
2686 template<class Gt, class Tds>
2687 inline typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
create_face()2688 Gt, Tds >::create_face()
2689 {
2690 return _tds.create_face();
2691 }
2692
2693 template<class Gt, class Tds>
2694 inline
delete_face(Face_handle f)2695 void Periodic_2_triangulation_2<Gt, Tds>::delete_face(Face_handle f)
2696 {
2697 _tds.delete_face(f);
2698 }
2699
2700 template<class Gt, class Tds>
2701 inline
delete_vertex(Vertex_handle v)2702 void Periodic_2_triangulation_2<Gt, Tds>::delete_vertex(Vertex_handle v)
2703 {
2704 _tds.delete_vertex(v);
2705 }
2706
2707 template<class Gt, class Tds>
compare_walks(const Point & p,Face_handle c1,Face_handle c2,Locate_type & lt1,Locate_type & lt2,int li1,int li2)2708 bool Periodic_2_triangulation_2<Gt, Tds>::compare_walks(const Point& p,
2709 Face_handle c1, Face_handle c2, Locate_type& lt1, Locate_type& lt2,
2710 int li1, int li2) const
2711 {
2712 bool b = true;
2713 b &= (lt1 == lt2);
2714 if ((lt1 == lt2) && (lt1 == VERTEX))
2715 {
2716 b &= (c1->vertex(li1) == c2->vertex(li2));
2717 }
2718 else if ((lt1 == lt2) && (lt1 == EDGE))
2719 {
2720 b &= ((c1 == c2)
2721 || ((c1->neighbor(li1) == c2) && (c2->neighbor(li2) == c1)));
2722 }
2723 else if ((lt1 == lt2) && (lt1 == EMPTY))
2724 {
2725 // Skip
2726 }
2727 else
2728 {
2729 b &= (lt1 == lt2);
2730 b &= (lt1 == FACE);
2731 b &= (c1 == c2);
2732 }
2733
2734 if (!b)
2735 {
2736 std::cerr << "from compare_walks " << std::endl;
2737 std::cerr << "point " << p << std::endl;
2738 std::cerr << "locate 1 " << &*c1 << "\t" << lt1 << "\t" << li1 << std::endl;
2739 std::cerr << "locate 2 " << &*c2 << "\t" << lt2 << "\t" << li2 << std::endl;
2740 std::cerr << std::endl;
2741 show_face(c1);
2742 std::cerr << std::endl;
2743 show_face(c2);
2744 std::cerr << std::endl;
2745 }
2746
2747 CGAL_triangulation_assertion(b);
2748 return b;
2749 }
2750
2751 template<class Gt, class Tds>
2752 typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle
2753 Periodic_2_triangulation_2<Gt, Tds>::
march_locate_2D(Face_handle f,const Point & query,const Offset & o_p,Locate_type & lt,int & li)2754 march_locate_2D(Face_handle f, const Point& query,
2755 const Offset& o_p, Locate_type& lt, int& li) const
2756 {
2757 CGAL_assertion(!empty());
2758
2759 Offset off_query = o_p;
2760
2761 // Random generator
2762 boost::rand48 rng;
2763 boost::uniform_smallint<> two(0, 1);
2764 boost::variate_generator<boost::rand48&, boost::uniform_smallint<> > coin(rng, two);
2765
2766 // Give the point the best start-offset possible
2767 if (is_1_cover() && !f->has_zero_offsets())
2768 {
2769 int cumm_off = f->offset(0) | f->offset(1) | f->offset(2);
2770 if (((cumm_off & 2) == 2) &&
2771 (FT(2) * query.x() < (_domain.xmax() + _domain.xmin())))
2772 off_query += Offset(1, 0);
2773 if (((cumm_off & 1) == 1) &&
2774 (FT(2) * query.y() < (_domain.ymax() + _domain.ymin())))
2775 off_query += Offset(0, 1);
2776 }
2777
2778 Face_handle prev = Face_handle();
2779 int prev_index = 0;
2780 Offset off[3];
2781 Orientation o[3];
2782 while (1)
2783 {
2784 // Instead of testing its edges in a random order we do the following
2785 // until we find a neighbor to go further:
2786 // As we come from prev we do not have to check the edge leading to prev
2787 // Now we flip a coin in order to decide if we start checking the
2788 // edge before or the edge after the edge leading to prev
2789 int left_first = coin() % 2;
2790
2791 bool simplicity_criterion =
2792 f->has_zero_offsets() && off_query.is_null() && is_1_cover();
2793
2794 const Point *p[3] =
2795 {
2796 &f->vertex(0)->point(),
2797 &f->vertex(1)->point(),
2798 &f->vertex(2)->point()
2799 };
2800
2801 // Get the offsets
2802 if (!simplicity_criterion)
2803 {
2804 if (!is_1_cover())
2805 {
2806 // Just fetch the vertices of c as points with offsets
2807 for (int i = 0; i < 3; i++)
2808 {
2809 off[i] = get_offset(f, i);
2810 }
2811 }
2812 else
2813 {
2814 // We are on the one cover and on the boundary between domains
2815 // Hence, we need to check predicates with offsets
2816 for (int i = 0; i < 3; i++)
2817 {
2818 off[i] = int_to_off(f->offset(i));
2819 }
2820 }
2821 }
2822
2823 if (prev == Face_handle())
2824 {
2825 prev = f;
2826 // First step, also check the prev_index
2827 if (simplicity_criterion)
2828 {
2829 o[ccw(prev_index)] =
2830 orientation(*p[ccw(prev_index)], *p[cw(prev_index)], query);
2831 }
2832 else
2833 {
2834 o[ccw(prev_index)] =
2835 orientation(*p[ccw(prev_index)], *p[cw(prev_index)], query,
2836 off[ccw(prev_index)], off[cw(prev_index)], off_query);
2837 }
2838 if (o[ccw(prev_index)] == NEGATIVE)
2839 {
2840 // This assignment is already done: prev = f
2841 f = f->neighbor(prev_index);
2842 int new_index = f->index(prev);
2843 if (!(simplicity_criterion && f->has_zero_offsets()))
2844 off_query = combine_offsets(off_query,
2845 get_neighbor_offset(prev, prev_index,
2846 f, new_index));
2847 prev_index = new_index;
2848 continue;
2849 }
2850 }
2851 else
2852 {
2853 o[ccw(prev_index)] = POSITIVE;
2854 }
2855
2856 if (left_first)
2857 {
2858 if (simplicity_criterion)
2859 {
2860 o[prev_index] =
2861 orientation(*p[prev_index], *p[ccw(prev_index)], query);
2862 }
2863 else
2864 {
2865 o[prev_index] =
2866 orientation(*p[prev_index], *p[ccw(prev_index)], query,
2867 off[prev_index], off[ccw(prev_index)], off_query);
2868 }
2869 if (o[prev_index] == NEGATIVE)
2870 {
2871 prev = f;
2872 f = f->neighbor(cw(prev_index));
2873 int new_index = f->index(prev);
2874 if (!(simplicity_criterion && f->has_zero_offsets()))
2875 off_query = combine_offsets(off_query,
2876 get_neighbor_offset(prev, cw(prev_index), f, new_index));
2877 prev_index = new_index;
2878 continue;
2879 }
2880 }
2881 {
2882 // Do right side
2883 if (simplicity_criterion)
2884 {
2885 o[cw(prev_index)] =
2886 orientation(*p[cw(prev_index)], *p[prev_index], query);
2887 }
2888 else
2889 {
2890 o[cw(prev_index)] =
2891 orientation(*p[cw(prev_index)], *p[prev_index], query,
2892 off[cw(prev_index)], off[prev_index], off_query);
2893 }
2894 if (o[cw(prev_index)] == NEGATIVE)
2895 {
2896 prev = f;
2897 f = f->neighbor(ccw(prev_index));
2898 int new_index = f->index(prev);
2899 if (!(simplicity_criterion && f->has_zero_offsets()))
2900 off_query = combine_offsets(off_query,
2901 get_neighbor_offset(prev, ccw(prev_index), f, new_index));
2902 prev_index = new_index;
2903 continue;
2904 }
2905 }
2906 if (!left_first)
2907 {
2908 if (simplicity_criterion)
2909 {
2910 o[prev_index] = orientation(*p[prev_index], *p[ccw(prev_index)], query);
2911 }
2912 else
2913 {
2914 o[prev_index] = orientation(*p[prev_index], *p[ccw(prev_index)], query,
2915 off[prev_index], off[ccw(prev_index)], off_query);
2916 }
2917 if (o[prev_index] == NEGATIVE)
2918 {
2919 prev = f;
2920 f = f->neighbor(cw(prev_index));
2921 int new_index = f->index(prev);
2922 if (!(simplicity_criterion && f->has_zero_offsets()))
2923 off_query = combine_offsets(off_query,
2924 get_neighbor_offset(prev, cw(prev_index), f, new_index));
2925 prev_index = new_index;
2926 continue;
2927 }
2928 }
2929
2930 // now p is in c or on its boundary
2931 int sum = (o[0] == COLLINEAR) + (o[1] == COLLINEAR) + (o[2] == COLLINEAR);
2932 switch (sum)
2933 {
2934 case 0:
2935 {
2936 lt = FACE;
2937 li = 4;
2938 break;
2939 }
2940 case 1:
2941 {
2942 lt = EDGE;
2943 li = (o[0] == COLLINEAR) ? 2 : (o[1] == COLLINEAR) ? 0 : 1;
2944 break;
2945 }
2946 case 2:
2947 {
2948 lt = VERTEX;
2949 li = (o[0] != COLLINEAR) ? 2 : (o[1] != COLLINEAR) ? 0 : 1;
2950 break;
2951 }
2952 }
2953 return f;
2954 }
2955 }
2956
2957 template<class Gt, class Tds>
2958 typename Periodic_2_triangulation_2<Gt, Tds>::Face_handle Periodic_2_triangulation_2 <
locate(const Point & p,const Offset & o,Locate_type & lt,int & li,Face_handle start)2959 Gt, Tds >::locate(const Point& p, const Offset &o, Locate_type& lt, int& li,
2960 Face_handle start) const
2961 {
2962 CGAL_triangulation_assertion((_domain.xmin() <= p.x()) &&
2963 (p.x() < _domain.xmax()));
2964 CGAL_triangulation_assertion((_domain.ymin() <= p.y()) &&
2965 (p.y() < _domain.ymax()));
2966
2967 if (dimension() <= 0)
2968 {
2969 lt = EMPTY;
2970 li = 4;
2971 return Face_handle();
2972 }
2973
2974 // Triangulation is not empty
2975 if (start == Face_handle())
2976 {
2977 start = faces_begin();
2978 }
2979
2980 return march_locate_2D(start, p, o, lt, li);
2981 }
2982
2983 /** Delete each redundant face and the not anymore needed data
2984 * structures.
2985 *
2986 * This function consists of four iterations over all faces and one
2987 * iteration over all vertices:
2988 * -# Face iteration: mark all faces that are to delete
2989 * -# Face iteration: redirect neighbors of remaining faces
2990 * -# Face iteration: redirect vertices of remaining faces
2991 * -# Face iteration: delete all faces marked in the 1. iteration
2992 * -# Vertex iteration: delete all vertices outside the original domain.
2993 */
2994 template<class Gt, class Tds>
convert_to_1_sheeted_covering()2995 void Periodic_2_triangulation_2<Gt, Tds>::convert_to_1_sheeted_covering()
2996 {
2997 // ###################################################################
2998 // ### First face iteration ##########################################
2999 // ###################################################################
3000 {
3001 if (is_1_cover())
3002 return;
3003 CGAL_triangulation_expensive_assertion(is_triangulation_in_1_sheet());
3004
3005 bool to_delete, has_simplifiable_offset;
3006 Virtual_vertex_map_it vvmit;
3007 // First iteration over all faces: Mark the faces that are to delete.
3008 // Faces are to delete if they cannot be translated anymore in the
3009 // direction of one of the axes without yielding negative offsets.
3010 for (Face_iterator it = all_faces_begin(); it != all_faces_end(); ++it)
3011 {
3012 to_delete = false;
3013 // for all directions in 2D Space
3014 for (int j = 0; j < 2; j++)
3015 {
3016 has_simplifiable_offset = true;
3017 // for all vertices of face it
3018 for (int i = 0; i < 3; i++)
3019 {
3020 vvmit = _virtual_vertices.find(it->vertex(i));
3021 if (vvmit == _virtual_vertices.end())
3022 {
3023 // if it->vertex(i) lies inside the original domain:
3024
3025 // the face cannot be moved any more because if we did, then
3026 // it->vertex(i) will get at least one negative offset.
3027 has_simplifiable_offset = false;
3028 }
3029 else
3030 {
3031 // if it->vertex(i) lies outside the original domain:
3032
3033 // The face can certainly be deleted if the offset contains a 2
3034 to_delete |= (vvmit->second.second[j] == 2);
3035 // The face can be moved into one direction only if the offset of
3036 // all for vertices is >=1 for this direction. Since we already
3037 // tested for 2 it is sufficient to test here for 1.
3038 has_simplifiable_offset &= (vvmit->second.second[j] == 1);
3039 }
3040 }
3041 // if the offset can be simplified, i.e. the face can be moved, then
3042 // it can be deleted.
3043 if (has_simplifiable_offset)
3044 to_delete = true;
3045 }
3046 // Mark all faces that are to delete. They cannot be deleted yet,
3047 // because neighboring information still needs to be extracted.
3048 it->set_additional_flag(to_delete ? 1 : 0);
3049 }
3050 }
3051
3052 // ###################################################################
3053 // ### Second face iteration #########################################
3054 // ###################################################################
3055 {
3056 Vertex_handle vert[3], nbv[3];
3057 Offset off[3];
3058 Face_handle nb, new_neighbor;
3059 std::vector<Triple<Face_handle, int, Face_handle> > new_neighbor_relations;
3060
3061 // Second iteration over all faces: redirect neighbors where necessary
3062 for (Face_iterator it = all_faces_begin(); it != all_faces_end(); ++it)
3063 {
3064 // Skip all faces that are to delete.
3065 if (it->get_additional_flag() == 1)
3066 continue;
3067
3068 // Redirect neighbors: Only neighbors that are marked by the
3069 // additional_flag have to be substituted by one of their periodic
3070 // copies. The unmarked neighbors stay the same.
3071 for (int i = 0; i < 3; i++)
3072 {
3073 if (it->neighbor(i)->get_additional_flag() != 1)
3074 continue;
3075
3076 nb = it->neighbor(i);
3077
3078 for (int j = 0; j < 3; j++)
3079 {
3080 off[j] = Offset();
3081 get_vertex(nb, j, vert[j], off[j]);
3082 }
3083 int x, y;
3084 x = (std::min)((std::min)(off[0][0], off[1][0]), off[2][0]);
3085 y = (std::min)((std::min)(off[0][1], off[1][1]), off[2][1]);
3086
3087 // The vector from nb to the "original" periodic copy of nb, that is
3088 // the copy that will not be deleted.
3089 Offset difference_offset(x, y);
3090 CGAL_triangulation_assertion( !difference_offset.is_null() );
3091
3092 // We now have to find the "original" periodic copy of nb from
3093 // its vertices. Therefore, we first have to find the vertices.
3094 for (int j = 0; j < 3; j++)
3095 {
3096 CGAL_triangulation_assertion( (off[j] - difference_offset)[0] >= 0);
3097 CGAL_triangulation_assertion( (off[j] - difference_offset)[1] >= 0);
3098 CGAL_triangulation_assertion( (off[j] - difference_offset)[0] < 3);
3099 CGAL_triangulation_assertion( (off[j] - difference_offset)[1] < 3);
3100
3101 // find the Vertex_handles of the vertices of the "original"
3102 // periodic copy of nb. If the vertex is inside the original
3103 // domain, there is nothing to do
3104 if ((off[j] - difference_offset).is_null())
3105 {
3106 nbv[j] = vert[j];
3107 // If the vertex is outside the original domain, we have to search
3108 // in _virtual_vertices in the "wrong" direction. That means, we
3109 // cannot use _virtual_vertices.find but have to use
3110 // _virtual_vertices_reverse.
3111 }
3112 else
3113 {
3114 Offset nbo = off[j] - difference_offset;
3115 nbv[j] = _virtual_vertices_reverse.find(vert[j]) ->second[nbo[0]
3116 * 3 + nbo[1] - 1];
3117 }
3118 }
3119 // Find the new neighbor by its 4 vertices
3120 new_neighbor = get_face(nbv);
3121
3122 // Store the new neighbor relation. This cannot be applied yet because
3123 // it would disturb the functioning of get_face( ... )
3124 new_neighbor_relations.push_back(make_triple(it, i, new_neighbor));
3125 }
3126 }
3127 // Apply the new neighbor relations now.
3128 for (unsigned int i = 0; i < new_neighbor_relations.size(); i++)
3129 {
3130 new_neighbor_relations[i].first->set_neighbor(
3131 new_neighbor_relations[i].second, new_neighbor_relations[i].third);
3132 }
3133 }
3134
3135 // ###################################################################
3136 // ### Third face iteration ##########################################
3137 // ###################################################################
3138 {
3139 Vertex_handle vert[3];
3140 Offset off[3];
3141 // Third iteration over all faces: redirect vertices where necessary
3142 for (Face_iterator it = all_faces_begin(); it != all_faces_end(); ++it)
3143 {
3144 // Skip all faces that are marked to delete
3145 if (it->get_additional_flag() == 1)
3146 continue;
3147 // Find the corresponding vertices of it in the original domain
3148 // and set them as new vertices of it.
3149 for (int i = 0; i < 3; i++)
3150 {
3151 off[i] = Offset();
3152 get_vertex(it, i, vert[i], off[i]);
3153 it->set_vertex(i, vert[i]);
3154 CGAL_triangulation_assertion(vert[i]->point()[0] < _domain.xmax());
3155 CGAL_triangulation_assertion(vert[i]->point()[1] < _domain.ymax());
3156 CGAL_triangulation_assertion(vert[i]->point()[0] >= _domain.xmin());
3157 CGAL_triangulation_assertion(vert[i]->point()[1] >= _domain.ymin());
3158
3159 // redirect also the face pointer of the vertex.
3160 it->vertex(i)->set_face(it);
3161 }
3162 // Set the offsets.
3163 set_offsets(it, off[0], off[1], off[2]);
3164 CGAL_triangulation_assertion( int_to_off(it->offset(0)) == off[0] );
3165 CGAL_triangulation_assertion( int_to_off(it->offset(1)) == off[1] );
3166 CGAL_triangulation_assertion( int_to_off(it->offset(2)) == off[2] );
3167 }
3168 }
3169
3170 // ###################################################################
3171 // ### Fourth face iteration #########################################
3172 // ###################################################################
3173 {
3174 // Delete the marked faces.
3175 std::vector<Face_handle> faces_to_delete;
3176 for (Face_iterator fit = all_faces_begin(); fit != all_faces_end(); ++fit)
3177 {
3178 if (fit->get_additional_flag() == 1)
3179 faces_to_delete.push_back(fit);
3180 }
3181 for (typename std::vector<Face_handle>::iterator it =
3182 faces_to_delete.begin(); it != faces_to_delete.end(); ++it)
3183 {
3184 _tds.delete_face(*it);
3185 }
3186 }
3187
3188 // ###################################################################
3189 // ### Vertex iteration ##############################################
3190 // ###################################################################
3191 {
3192 // Delete all the vertices in _virtual_vertices, that is all vertices
3193 // outside the original domain.
3194 std::vector<Vertex_handle> vertices_to_delete;
3195 for (Vertex_iterator vit = all_vertices_begin(); vit != all_vertices_end(); ++vit)
3196 {
3197 if (_virtual_vertices.count(vit) != 0)
3198 {
3199 CGAL_triangulation_assertion( _virtual_vertices.count( vit ) == 1 );
3200 vertices_to_delete.push_back(vit);
3201 }
3202 }
3203 for (typename std::vector<Vertex_handle>::iterator it =
3204 vertices_to_delete.begin(); it != vertices_to_delete.end(); ++it)
3205 {
3206 _tds.delete_vertex(*it);
3207 }
3208 }
3209
3210 _cover = make_array(1, 1);
3211 _virtual_vertices.clear();
3212 _virtual_vertices_reverse.clear();
3213 _too_long_edge_counter = 0;
3214 _too_long_edges.clear();
3215
3216 CGAL_triangulation_assertion(is_1_cover());
3217 }
3218
3219 template<class Gt, class Tds>
convert_to_9_sheeted_covering()3220 void Periodic_2_triangulation_2<Gt, Tds>::convert_to_9_sheeted_covering()
3221 {
3222 if (_cover == make_array(3, 3))
3223 return;
3224 CGAL_triangulation_precondition(is_1_cover());
3225
3226 // Create 9 copies of each vertex and write virtual_vertices and
3227 // virtual_vertices_reverse
3228 std::list<Vertex_handle> original_vertices;
3229 // try to use std::copy instead of the following loop.
3230 for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit)
3231 original_vertices.push_back(vit);
3232 for (typename std::list<Vertex_handle>::iterator vit =
3233 original_vertices.begin(); vit != original_vertices.end(); ++vit)
3234 {
3235 Vertex_handle v_cp;
3236 std::vector<Vertex_handle> copies;
3237 for (int i = 0; i < 3; i++)
3238 for (int j = 0; j < 3; j++)
3239 {
3240 if (i == 0 && j == 0)
3241 continue;
3242 v_cp = _tds.create_vertex(*vit);
3243 copies.push_back(v_cp);
3244 _virtual_vertices.insert(std::make_pair(v_cp, std::make_pair(*vit,
3245 Offset(i, j))));
3246 }
3247 _virtual_vertices_reverse.insert(std::make_pair(*vit, copies));
3248 }
3249
3250 // Create 9 copies of each face from the respective copies of the
3251 // vertices and write virtual_faces and virtual_faces_reverse.
3252 typedef std::map<Face_handle, std::pair<Face_handle, Offset> >
3253 Virtual_face_map;
3254 typedef std::map<Face_handle, std::vector<Face_handle> >
3255 Virtual_face_reverse_map;
3256 typedef typename Virtual_face_reverse_map::const_iterator VCRMIT;
3257
3258 Virtual_face_map virtual_faces;
3259 Virtual_face_reverse_map virtual_faces_reverse;
3260
3261 std::list<Face_handle> original_faces;
3262 for (Face_iterator fit = faces_begin(); fit != faces_end(); ++fit)
3263 original_faces.push_back(fit);
3264
3265 // Store vertex offsets in a separate data structure
3266 std::list<Offset> off_v;
3267 for (typename std::list<Vertex_handle>::iterator vit =
3268 original_vertices.begin(); vit != original_vertices.end(); ++vit)
3269 {
3270 Face_handle ccc = (*vit)->face();
3271 int v_index = ccc->index(*vit);
3272 off_v.push_back(int_to_off(ccc->offset(v_index)));
3273 }
3274
3275 // Store neighboring offsets in a separate data structure
3276 std::list<std::array<Offset, 3> > off_nb;
3277 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3278 != original_faces.end(); ++fit)
3279 {
3280 std::array<Offset, 3> off_nb_f;
3281 for (int i = 0; i < 3; i++)
3282 {
3283 Face_handle fff = *fit;
3284 Face_handle nnn = fff->neighbor(i);
3285 off_nb_f[i] = get_neighbor_offset(fff, i, nnn, nnn->index(fff));
3286 }
3287 off_nb.push_back(off_nb_f);
3288 }
3289
3290 // Create copies of faces
3291 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3292 != original_faces.end(); ++fit)
3293 {
3294 Face_handle c_cp;
3295 Vertex_handle v0, v1, v2;
3296 std::vector<Face_handle> copies;
3297 Virtual_vertex_reverse_map_it vvrmit[3];
3298 Offset vvoff[3];
3299 for (int i = 0; i < 3; i++)
3300 {
3301 vvrmit[i] = _virtual_vertices_reverse.find((*fit)->vertex(i));
3302 CGAL_triangulation_assertion(
3303 vvrmit[i] != _virtual_vertices_reverse.end());
3304 vvoff[i] = int_to_off((*fit)->offset(i));
3305 }
3306 Vertex_handle vvh[3];
3307 for (int n = 0; n < 8; n++) // iterate over faces
3308 {
3309 for (int i = 0; i < 3; i++) // iterate over vertices of the face
3310 {
3311 // Decomposition of n into an offset (nx,ny):
3312 // nx = ((n+1)/3)%3, ny = (n+1)%3
3313 int o_i = ((n + 1) / 3 + vvoff[i].x() + 3) % 3;
3314 int o_j = ((n + 1) + vvoff[i].y() + 3) % 3;
3315 int n_c = 3 * o_i + o_j - 1;
3316 CGAL_triangulation_assertion(n_c >= -1);
3317 if (n_c == -1)
3318 vvh[i] = (*fit)->vertex(i);
3319 else
3320 vvh[i] = vvrmit[i]->second[n_c];
3321 }
3322 c_cp = _tds.create_face(vvh[0], vvh[1], vvh[2]);
3323 copies.push_back(c_cp);
3324 }
3325 virtual_faces_reverse.insert(std::make_pair(*fit, copies));
3326 }
3327
3328 // Set new vertices of boundary faces of the original domain.
3329 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3330 != original_faces.end(); ++fit)
3331 {
3332 for (int i = 0; i < 3; i++)
3333 {
3334 Virtual_vertex_reverse_map_it vvrmit = _virtual_vertices_reverse.find(
3335 (*fit)->vertex(i));
3336 CGAL_triangulation_assertion(vvrmit != _virtual_vertices_reverse.end());
3337 Offset vvoff = int_to_off((*fit)->offset(i));
3338 if (!vvoff.is_null())
3339 {
3340 int n_f = 3 * vvoff.x() + vvoff.y() - 1;
3341 CGAL_triangulation_assertion(n_f >= 0);
3342 CGAL_triangulation_assertion(static_cast<unsigned int>(n_f) < vvrmit->second.size());
3343 (*fit)->set_vertex(i, vvrmit->second[n_f]);
3344 }
3345 }
3346 }
3347
3348 // Set neighboring relations of face copies
3349 typename std::list<std::array<Offset, 3> >::iterator oit = off_nb.begin();
3350 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3351 != original_faces.end(); ++fit, ++oit)
3352 {
3353 CGAL_triangulation_assertion( oit != off_nb.end() );
3354 VCRMIT c_cp = virtual_faces_reverse.find(*fit);
3355 CGAL_triangulation_assertion(c_cp != virtual_faces_reverse.end());
3356 for (int i = 0; i < 3; i++)
3357 {
3358 Face_handle fit_nb = (*fit)->neighbor(i);
3359 VCRMIT c_cp_nb = virtual_faces_reverse.find(fit_nb);
3360 CGAL_triangulation_assertion(c_cp_nb != virtual_faces_reverse.end());
3361 Offset nboff = (*oit)[i];
3362 for (int n = 0; n < 8; n++)
3363 {
3364 int n_nb;
3365 if (nboff.is_null())
3366 n_nb = n;
3367 else
3368 {
3369 int o_i = ((n + 1) / 3 - nboff.x() + 3) % 3;
3370 int o_j = (n + 1 - nboff.y() + 3) % 3;
3371 n_nb = 3 * o_i + o_j - 1;
3372 }
3373 if (n_nb == -1)
3374 {
3375 CGAL_triangulation_assertion(fit_nb->has_vertex(c_cp->second[n]->vertex(ccw(i))) );
3376 CGAL_triangulation_assertion(fit_nb->has_vertex(c_cp->second[n]->vertex( cw(i))) );
3377 c_cp->second[n]->set_neighbor(i, fit_nb);
3378 }
3379 else
3380 {
3381 CGAL_triangulation_assertion(n_nb >= 0);
3382 CGAL_triangulation_assertion(static_cast<unsigned int>(n_nb) <= c_cp_nb->second.size());
3383 CGAL_triangulation_assertion(c_cp_nb->second[n_nb]->has_vertex(c_cp->second[n]->vertex(ccw(i))) );
3384 CGAL_triangulation_assertion(c_cp_nb->second[n_nb]->has_vertex(c_cp->second[n]->vertex( cw(i))) );
3385 c_cp->second[n]->set_neighbor(i, c_cp_nb->second[n_nb]);
3386 }
3387 }
3388 }
3389 }
3390
3391 // Set neighboring relations of original faces
3392 oit = off_nb.begin();
3393 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3394 != original_faces.end(); ++fit, ++oit)
3395 {
3396 CGAL_triangulation_assertion( oit != off_nb.end() );
3397 for (int i = 0; i < 3; i++)
3398 {
3399 Offset nboff = (*oit)[i];
3400 if (!nboff.is_null())
3401 {
3402 Face_handle fit_nb = (*fit)->neighbor(i);
3403 VCRMIT c_cp_nb = virtual_faces_reverse.find(fit_nb);
3404 CGAL_triangulation_assertion(c_cp_nb != virtual_faces_reverse.end());
3405 int o_i = (3 - nboff.x()) % 3;
3406 int o_j = (3 - nboff.y()) % 3;
3407 int n_nb = 3 * o_i + o_j - 1;
3408 CGAL_triangulation_assertion(n_nb >= 0);
3409 CGAL_triangulation_assertion(static_cast<unsigned int>(n_nb) <= c_cp_nb->second.size());
3410 CGAL_triangulation_assertion(c_cp_nb->second[n_nb]->has_vertex((*fit)->vertex(ccw(i))) );
3411 CGAL_triangulation_assertion(c_cp_nb->second[n_nb]->has_vertex((*fit)->vertex( cw(i))) );
3412 (*fit)->set_neighbor(i, c_cp_nb->second[n_nb]);
3413 }
3414 }
3415 }
3416
3417 // Set incident faces
3418 for (Face_iterator fit = faces_begin(); fit != faces_end(); ++fit)
3419 {
3420 for (int i = 0; i < 3; i++)
3421 {
3422 fit->vertex(i)->set_face(fit);
3423 }
3424 }
3425
3426 // Set offsets where necessary
3427 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3428 != original_faces.end(); ++fit)
3429 {
3430 VCRMIT c_cp = virtual_faces_reverse.find(*fit);
3431 CGAL_triangulation_assertion( c_cp != virtual_faces_reverse.end());
3432 Offset off[3];
3433 for (int i = 0; i < 3; i++)
3434 off[i] = int_to_off((*fit)->offset(i));
3435 if (off[0].is_null() && off[1].is_null() && off[2].is_null())
3436 continue;
3437 for (int n = 0; n < 8; n++)
3438 {
3439 Offset off_cp[4];
3440 int o_i = ((n + 1) / 3) % 3;
3441 int o_j = (n + 1) % 3;
3442 if (o_i != 2 && o_j != 2)
3443 continue;
3444 for (int i = 0; i < 3; i++)
3445 {
3446 off_cp[i] = Offset((o_i == 2) ? off[i].x() : 0, (o_j == 2) ? off[i].y()
3447 : 0);
3448 CGAL_triangulation_assertion(off_cp[i].x() == 0 || off_cp[i].x() == 1);
3449 CGAL_triangulation_assertion(off_cp[i].y() == 0 || off_cp[i].y() == 1);
3450 }
3451 set_offsets(c_cp->second[n], off_cp[0], off_cp[1], off_cp[2]);
3452 }
3453 }
3454
3455 // Iterate over all original faces and reset offsets.
3456 for (typename std::list<Face_handle>::iterator fit = original_faces.begin(); fit
3457 != original_faces.end(); ++fit)
3458 {
3459 //This statement does not seem to have any effect
3460 set_offsets(*fit, 0, 0, 0);
3461 CGAL_triangulation_assertion((*fit)->offset(0) == 0);
3462 CGAL_triangulation_assertion((*fit)->offset(1) == 0);
3463 CGAL_triangulation_assertion((*fit)->offset(2) == 0);
3464 }
3465
3466 _cover = make_array(3, 3);
3467
3468 // Set up too long edges data structure
3469 int i = 0;
3470 for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit)
3471 {
3472 _too_long_edges[vit] = std::list<Vertex_handle>();
3473 ++i;
3474 }
3475 _too_long_edge_counter = find_too_long_edges(_too_long_edges);
3476
3477 CGAL_triangulation_expensive_assertion(is_valid());
3478
3479 CGAL_triangulation_assertion(!is_1_cover());
3480 }
3481
3482 // iterate over all edges and store the ones that are longer than
3483 // edge_length_threshold in edges. Return the number of too long edges.
3484 template<class GT, class Tds>
find_too_long_edges(std::map<Vertex_handle,std::list<Vertex_handle>> & edges)3485 inline int Periodic_2_triangulation_2<GT, Tds>::find_too_long_edges(std::map <
3486 Vertex_handle, std::list<Vertex_handle> > & edges) const
3487 {
3488 Point p1, p2;
3489 int counter = 0;
3490 Vertex_handle v_no, vh;
3491 for (Edge_iterator eit = edges_begin(); eit != edges_end(); eit++)
3492 {
3493 p1 = construct_point(eit->first->vertex(ccw(eit->second))->point(),
3494 get_offset(eit->first, ccw(eit->second)));
3495 p2 = construct_point(eit->first->vertex(cw(eit->second))->point(),
3496 get_offset(eit->first, cw(eit->second)));
3497 if (edge_is_too_long(p1, p2))
3498 {
3499 if (&*(eit->first->vertex(ccw(eit->second))) < &*(eit->first->vertex(cw(
3500 eit->second))))
3501 {
3502 v_no = eit->first->vertex(ccw(eit->second));
3503 vh = eit->first->vertex(cw(eit->second));
3504 }
3505 else
3506 {
3507 v_no = eit->first->vertex(cw(eit->second));
3508 vh = eit->first->vertex(ccw(eit->second));
3509 }
3510 edges[v_no].push_back(vh);
3511 counter++;
3512 }
3513 }
3514 return counter;
3515 }
3516
3517 /**
3518 * - fh->offset(i) is an bit tuple encapsulated in an integer. Each bit
3519 * represents the offset in one direction --> 2-cover!
3520 * - int_to_off(int) decodes this again.
3521 * - Finally the offset vector is multiplied by cover.
3522 * So if we are working in 3-cover we translate it to the neighboring
3523 * 3-cover and not only to the neighboring domain.
3524 */
3525 template<class GT, class Tds>
get_vertex(Face_handle fh,int i,Vertex_handle & vh,Offset & off)3526 inline void Periodic_2_triangulation_2<GT, Tds>::get_vertex(Face_handle fh,
3527 int i, Vertex_handle &vh, Offset &off) const
3528 {
3529 off = combine_offsets(Offset(), int_to_off(fh->offset(i)));
3530 vh = fh->vertex(i);
3531
3532 if (is_1_cover())
3533 return;
3534 Vertex_handle vh_i = vh;
3535 get_vertex(vh_i, vh, off);
3536 return;
3537 }
3538
3539 template<class GT, class Tds>
get_vertex(Vertex_handle vh_i,Vertex_handle & vh,Offset & off)3540 inline void Periodic_2_triangulation_2<GT, Tds>::get_vertex(Vertex_handle vh_i,
3541 Vertex_handle &vh, Offset &off) const
3542 {
3543 Virtual_vertex_map_it it = _virtual_vertices.find(vh_i);
3544
3545 if (it == _virtual_vertices.end())
3546 {
3547 // if vh_i is not contained in virtual_vertices, then it is in the
3548 // original domain.
3549 vh = vh_i;
3550 CGAL_triangulation_assertion(vh != Vertex_handle());
3551 }
3552 else
3553 {
3554 // otherwise it has to be looked up as well as its offset.
3555 vh = it->second.first;
3556 off += it->second.second;
3557 CGAL_triangulation_assertion(vh->point().x() < _domain.xmax());
3558 CGAL_triangulation_assertion(vh->point().y() < _domain.ymax());
3559 CGAL_triangulation_assertion(vh->point().x() >= _domain.xmin());
3560 CGAL_triangulation_assertion(vh->point().y() >= _domain.ymin());
3561 }
3562 }
3563
3564 /** Find the Face that consists of the three given vertices
3565 *
3566 * Iterates over all faces and compare the three vertices of each face
3567 * with the three vertices in vh.
3568 */
3569 template<class GT, class Tds>
3570 inline typename Periodic_2_triangulation_2<GT, Tds>::Face_handle Periodic_2_triangulation_2 <
get_face(const Vertex_handle * vh)3571 GT, Tds >::get_face(const Vertex_handle* vh) const
3572 {
3573 bool contains_v[2];
3574 Face_circulator fc = incident_faces(vh[2]);
3575 Face_circulator done(fc);
3576 do
3577 {
3578 CGAL_triangulation_assertion(
3579 fc->vertex(0) == vh[2] ||
3580 fc->vertex(1) == vh[2] ||
3581 fc->vertex(2) == vh[2]);
3582 for (int j = 0; j < 2; j++)
3583 {
3584 contains_v[j] = (fc->vertex(0) == vh[j]) || (fc->vertex(1) == vh[j])
3585 || (fc->vertex(2) == vh[j]);
3586 }
3587 if (contains_v[0] && contains_v[1])
3588 {
3589 return fc;
3590 }
3591 }
3592 while (++fc != done);
3593
3594 CGAL_triangulation_assertion(false);
3595 return Face_handle();
3596 }
3597
3598 template<class Gt, class Tds>
side_of_face(const Point & q,const Offset & off,Face_handle f,Locate_type & lt,int & li)3599 Bounded_side Periodic_2_triangulation_2<Gt, Tds>::side_of_face(const Point &q,
3600 const Offset &off, Face_handle f, Locate_type <, int &li) const
3601 {
3602 CGAL_triangulation_precondition(number_of_vertices() != 0);
3603
3604 Orientation o0, o1, o2;
3605 o0 = o1 = o2 = ZERO;
3606 int cumm_off = f->offset(0) | f->offset(1) | f->offset(2);
3607
3608 if ((cumm_off == 0) && is_1_cover())
3609 {
3610 CGAL_triangulation_assertion(off == Offset());
3611
3612 const Point &p0 = f->vertex(0)->point();
3613 const Point &p1 = f->vertex(1)->point();
3614 const Point &p2 = f->vertex(2)->point();
3615
3616 if (((o0 = orientation(q, p1, p2)) == NEGATIVE) || ((o1 = orientation(p0,
3617 q, p2)) == NEGATIVE) || ((o2 = orientation(p0, p1, q)) == NEGATIVE))
3618 {
3619 return ON_UNBOUNDED_SIDE;
3620 }
3621 }
3622 else // Special case for the periodic space.
3623 {
3624 Offset off_q;
3625 Offset offs[3];
3626 const Point *p[3];
3627 for (int i = 0; i < 3; i++)
3628 {
3629 p[i] = &(f->vertex(i)->point());
3630 offs[i] = get_offset(f, i);
3631 }
3632 CGAL_triangulation_assertion(orientation(*p[0], *p[1], *p[2],
3633 offs[0], offs[1], offs[2]) == POSITIVE);
3634 bool found = false;
3635 for (int i = 0; (i < 4) && (!found); i++)
3636 {
3637 if ((cumm_off | ((~i) & 3)) == 3)
3638 {
3639 o0 = o1 = o2 = NEGATIVE;
3640 off_q = combine_offsets(off, int_to_off(i));
3641
3642 if (((o0 = orientation(q, *p[1], *p[2], off_q, offs[1], offs[2]))
3643 != NEGATIVE) && ((o1 = orientation(*p[0], q, *p[2], offs[0], off_q,
3644 offs[2])) != NEGATIVE) && ((o2 = orientation(*p[0], *p[1], q,
3645 offs[0], offs[1], off_q)) != NEGATIVE))
3646 {
3647 found = true;
3648 }
3649 }
3650 }
3651 if (!found)
3652 return ON_UNBOUNDED_SIDE;
3653 }
3654
3655 // now all the oi's are >=0
3656 // sum gives the number of facets p lies on
3657 int sum = ((o0 == ZERO) ? 1 : 0) + ((o1 == ZERO) ? 1 : 0) + ((o2 == ZERO) ? 1
3658 : 0);
3659
3660 switch (sum)
3661 {
3662 case 0:
3663 {
3664 lt = FACE;
3665 return ON_BOUNDED_SIDE;
3666 }
3667 case 1:
3668 {
3669 lt = EDGE;
3670 // i = index such that q lies on edge (f,li)
3671 li = (o0 == ZERO) ? 0 : (o1 == ZERO) ? 1 : 2;
3672 return ON_BOUNDARY;
3673 }
3674 case 2:
3675 {
3676 lt = VERTEX;
3677 // i = index such that q lies on vertex li
3678 li = (o0 != ZERO) ? 0 : (o1 != ZERO) ? 1 : 2;
3679 return ON_BOUNDARY;
3680 }
3681 default:
3682 {
3683 // impossible : cannot be on 3 edges for a real triangle
3684 CGAL_triangulation_assertion(false);
3685 return ON_BOUNDARY;
3686 }
3687 }
3688 }
3689
3690 template<class Gt, class Tds>
oriented_side(Face_handle f,const Point & p,const Offset & o)3691 Oriented_side Periodic_2_triangulation_2<Gt, Tds>::oriented_side(Face_handle f,
3692 const Point& p, const Offset &o) const
3693 {
3694 Point &p0 = f->vertex(0)->point();
3695 Point &p1 = f->vertex(1)->point();
3696 Point &p2 = f->vertex(2)->point();
3697
3698 int cumm_off = f->offset(0) | f->offset(1) | f->offset(2);
3699
3700 if ((cumm_off == 0) && is_1_cover())
3701 {
3702 CGAL_precondition(o == Offset());
3703
3704 // return position of point p with respect to the oriented triangle p0p1p2
3705 // the orientation of the vertices is assumed to be counter clockwise
3706 CGAL_assertion(orientation(p0, p1, p2) == LEFT_TURN);
3707
3708 Bounded_side bs = bounded_side(p0, p1, p2, p);
3709 switch (bs)
3710 {
3711 case ON_BOUNDARY:
3712 return ON_ORIENTED_BOUNDARY;
3713 case ON_BOUNDED_SIDE:
3714 return ON_POSITIVE_SIDE;
3715 case ON_UNBOUNDED_SIDE:
3716 return ON_NEGATIVE_SIDE;
3717 }
3718 }
3719 else // Special case for the periodic space.
3720 {
3721 Offset off_q;
3722 Offset off0 = get_offset(f, 0);
3723 Offset off1 = get_offset(f, 1);
3724 Offset off2 = get_offset(f, 2);
3725
3726 // return position of point p with respect to the oriented triangle p0p1p2
3727 // the orientation of the vertices is assumed to be counter clockwise
3728 CGAL_assertion(orientation(p0, p1, p2, off0, off1, off2) == LEFT_TURN);
3729
3730 Bounded_side bs;
3731 for (int i = 0; (i <= 7); i++)
3732 {
3733 if ((cumm_off | ((~i) & 3)) == 3)
3734 {
3735 off_q = combine_offsets(o, int_to_off(i));
3736 bs = bounded_side(p0, p1, p2, p, off0, off1, off2, off_q);
3737 if (bs != ON_UNBOUNDED_SIDE)
3738 {
3739 return (bs == ON_BOUNDARY ? ON_ORIENTED_BOUNDARY : ON_POSITIVE_SIDE);
3740 }
3741 }
3742 }
3743
3744 return ON_NEGATIVE_SIDE;
3745 }
3746
3747 CGAL_assertion(false);
3748 return ON_NEGATIVE_SIDE;
3749 }
3750
3751 template<class Gt, class Tds>
bounded_side(const Point & p0,const Point & p1,const Point & p2,const Point & p)3752 Bounded_side Periodic_2_triangulation_2<Gt, Tds>::bounded_side(const Point &p0, const Point &p1, const Point &p2, const Point &p) const
3753 {
3754
3755 // return position of point p with respect to triangle p0p1p2
3756 CGAL_triangulation_precondition( orientation(p0, p1, p2) != COLLINEAR);
3757
3758 Orientation o1 = orientation(p0, p1, p);
3759 Orientation o2 = orientation(p1, p2, p);
3760 Orientation o3 = orientation(p2, p0, p);
3761
3762 if (o1 == COLLINEAR)
3763 {
3764 if (o2 == COLLINEAR || o3 == COLLINEAR)
3765 return ON_BOUNDARY;
3766 if (collinear_between(p0, p, p1))
3767 return ON_BOUNDARY;
3768 return ON_UNBOUNDED_SIDE;
3769 }
3770
3771 if (o2 == COLLINEAR)
3772 {
3773 if (o3 == COLLINEAR)
3774 return ON_BOUNDARY;
3775 if (collinear_between(p1, p, p2))
3776 return ON_BOUNDARY;
3777 return ON_UNBOUNDED_SIDE;
3778 }
3779
3780 if (o3 == COLLINEAR)
3781 {
3782 if (collinear_between(p2, p, p0))
3783 return ON_BOUNDARY;
3784 return ON_UNBOUNDED_SIDE;
3785 }
3786
3787 // from here none ot, o1, o2 and o3 are known to be non null
3788 if (o1 == o2 && o2 == o3)
3789 return ON_BOUNDED_SIDE;
3790 return ON_UNBOUNDED_SIDE;
3791 }
3792
3793 template<class Gt, class Tds>
bounded_side(const Point & p0,const Point & p1,const Point & p2,const Point & p,const Offset & o0,const Offset & o1,const Offset & o2,const Offset & o)3794 Bounded_side Periodic_2_triangulation_2<Gt, Tds>::bounded_side(const Point &p0,
3795 const Point &p1, const Point &p2, const Point &p, const Offset &o0,
3796 const Offset &o1, const Offset &o2, const Offset &o) const
3797 {
3798 // return position of point p with respect to triangle p0p1p2
3799 CGAL_triangulation_precondition( orientation(p0, p1, p2, o0, o1, o2) != COLLINEAR);
3800 Orientation orient1 = orientation(p0, p1, p, o0, o1, o);
3801 Orientation orient2 = orientation(p1, p2, p, o1, o2, o);
3802 Orientation orient3 = orientation(p2, p0, p, o2, o0, o);
3803
3804 if (orient1 == COLLINEAR)
3805 {
3806 if (orient2 == COLLINEAR || orient3 == COLLINEAR)
3807 return ON_BOUNDARY;
3808 if (collinear_between(p0, p, p1, o0, o, o1))
3809 return ON_BOUNDARY;
3810 return ON_UNBOUNDED_SIDE;
3811 }
3812
3813 if (orient2 == COLLINEAR)
3814 {
3815 if (orient3 == COLLINEAR)
3816 return ON_BOUNDARY;
3817 if (collinear_between(p1, p, p2, o1, o, o2))
3818 return ON_BOUNDARY;
3819 return ON_UNBOUNDED_SIDE;
3820 }
3821
3822 if (orient3 == COLLINEAR)
3823 {
3824 if (collinear_between(p2, p, p0, o2, o, o0))
3825 return ON_BOUNDARY;
3826 return ON_UNBOUNDED_SIDE;
3827 }
3828
3829 // from here none ot, o1, o2 and o3 are known to be non null
3830 if (orient1 == orient2 && orient2 == orient3)
3831 return ON_BOUNDED_SIDE;
3832 return ON_UNBOUNDED_SIDE;
3833 }
3834
3835
3836 template<class Gt, class Tds>
collinear_between(const Point & p,const Point & q,const Point & r)3837 bool Periodic_2_triangulation_2<Gt, Tds>::collinear_between(const Point& p,
3838 const Point& q, const Point& r) const
3839 {
3840 // return true if point q is strictly between p and r
3841 // p,q and r are supposed to be collinear points
3842 Comparison_result c_pr = compare_x(p, r);
3843 Comparison_result c_pq;
3844 Comparison_result c_qr;
3845 if(c_pr == EQUAL)
3846 {
3847 c_pq = compare_y(p, q);
3848 c_qr = compare_y(q, r);
3849 }
3850 else
3851 {
3852 c_pq = compare_x(p, q);
3853 c_qr = compare_x(q, r);
3854 }
3855 return ( (c_pq == SMALLER) && (c_qr == SMALLER) ) ||
3856 ( (c_pq == LARGER) && (c_qr == LARGER) );
3857
3858 }
3859
3860 template<class Gt, class Tds>
collinear_between(const Point & p,const Point & q,const Point & r,const Offset & o_p,const Offset & o_q,const Offset & o_r)3861 bool Periodic_2_triangulation_2<Gt, Tds>::collinear_between(const Point& p,
3862 const Point& q, const Point& r, const Offset& o_p, const Offset& o_q,
3863 const Offset& o_r) const
3864 {
3865 // return true if point q is strictly between p and r
3866 // p,q and r are supposed to be collinear points
3867 Comparison_result c_pr = compare_x(p, r, o_p, o_r);
3868 Comparison_result c_pq;
3869 Comparison_result c_qr;
3870 if (c_pr == EQUAL)
3871 {
3872 c_pq = compare_y(p, q, o_p, o_q);
3873 c_qr = compare_y(q, r, o_q, o_r);
3874 }
3875 else
3876 {
3877 c_pq = compare_x(p, q, o_p, o_q);
3878 c_qr = compare_x(q, r, o_q, o_r);
3879 }
3880 return (((c_pq == SMALLER) && (c_qr == SMALLER)) ||
3881 ((c_pq == LARGER) && (c_qr == LARGER)));
3882 }
3883
3884 template<class Gt, class Tds>
compare_x(const Point & p,const Point & q)3885 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_x(
3886 const Point& p, const Point& q) const
3887 {
3888 return geom_traits().compare_x_2_object()(p, q);
3889 }
3890
3891 template<class Gt, class Tds>
compare_x(const Point & p,const Point & q,const Offset & o_p,const Offset & o_q)3892 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_x(
3893 const Point& p, const Point& q, const Offset &o_p, const Offset &o_q) const
3894 {
3895 return geom_traits().compare_x_2_object()(p, q, o_p, o_q);
3896 }
3897
3898 template<class Gt, class Tds>
compare_xy(const Point & p,const Point & q)3899 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_xy(
3900 const Point& p, const Point& q) const
3901 {
3902 Comparison_result res = geom_traits().compare_x_2_object()(p, q);
3903 if (res == EQUAL)
3904 {
3905 return geom_traits().compare_y_2_object()(p, q);
3906 }
3907 return res;
3908 }
3909
3910 template<class Gt, class Tds>
compare_xy(const Point & p,const Point & q,const Offset & o_p,const Offset & o_q)3911 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_xy(
3912 const Point& p, const Point& q, const Offset &o_p, const Offset &o_q) const
3913 {
3914 Comparison_result res = geom_traits().compare_x_2_object()(p, q, o_p, o_q);
3915 if (res == EQUAL)
3916 {
3917 return geom_traits().compare_y_2_object()(p, q, o_p, o_q);
3918 }
3919 return res;
3920 }
3921
3922 template<class Gt, class Tds>
compare_y(const Point & p,const Point & q)3923 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_y(
3924 const Point& p, const Point& q) const
3925 {
3926 return geom_traits().compare_y_2_object()(p, q);
3927 }
3928
3929 template<class Gt, class Tds>
compare_y(const Point & p,const Point & q,const Offset & o_p,const Offset & o_q)3930 inline Comparison_result Periodic_2_triangulation_2<Gt, Tds>::compare_y(
3931 const Point& p, const Point& q, const Offset &o_p, const Offset &o_q) const
3932 {
3933 return geom_traits().compare_y_2_object()(p, q, o_p, o_q);
3934 }
3935
3936 template<class Gt, class Tds>
3937 inline
xy_equal(const Point & p,const Point & q)3938 bool Periodic_2_triangulation_2<Gt, Tds>::xy_equal(const Point& p,
3939 const Point& q) const
3940 {
3941 return compare_xy(p, q) == EQUAL;
3942 }
3943
3944 template<class Gt, class Tds>
orientation(const Point & p0,const Point & p1,const Point & p2)3945 inline Orientation Periodic_2_triangulation_2<Gt, Tds>::orientation(
3946 const Point& p0, const Point& p1, const Point& p2) const
3947 {
3948 return geom_traits().orientation_2_object()(p0, p1, p2);
3949 }
3950 template<class Gt, class Tds>
orientation(const Point & p0,const Point & p1,const Point & p2,const Offset & o0,const Offset & o1,const Offset & o2)3951 inline Orientation Periodic_2_triangulation_2<Gt, Tds>::orientation(
3952 const Point& p0, const Point& p1, const Point& p2, const Offset& o0,
3953 const Offset& o1, const Offset& o2) const
3954 {
3955 return geom_traits().orientation_2_object()(p0, p1, p2, o0, o1, o2);
3956 }
3957
3958 template<class Gt, class Tds>
insert_too_long_edges_in_star(Vertex_handle vh)3959 void Periodic_2_triangulation_2<Gt, Tds>::insert_too_long_edges_in_star(Vertex_handle vh)
3960 {
3961 // Insert the too long edges in the star of vh
3962 Face_handle f = vh->face();
3963 Face_handle f_start = f;
3964
3965 do
3966 {
3967 int i = ccw(f->index(vh));
3968
3969 insert_too_long_edge(f, i);
3970
3971 // Proceed to the next face
3972 f = f->neighbor(i);
3973 }
3974 while (f != f_start);
3975 }
3976
3977 template<class Gt, class Tds>
insert_too_long_edge(Face_handle f,int i)3978 void Periodic_2_triangulation_2<Gt, Tds>::insert_too_long_edge(Face_handle f, int i)
3979 {
3980 Vertex_handle vh1 = f->vertex(ccw(i));
3981 Vertex_handle vh2 = f->vertex(cw(i));
3982 CGAL_assertion(vh1 != Vertex_handle());
3983 CGAL_assertion(vh2 != Vertex_handle());
3984 Point p1 = construct_point(vh1->point(), get_offset(f, ccw(i)));
3985 Point p2 = construct_point(vh2->point(), get_offset(f, cw(i)));
3986
3987 if (&*vh1 < &*vh2)
3988 {
3989 if (edge_is_too_long(p1, p2) &&
3990 (find(_too_long_edges[vh1].begin(), _too_long_edges[vh1].end(), vh2) == _too_long_edges[vh1].end()))
3991 {
3992 _too_long_edges[vh1].push_back(vh2);
3993 _too_long_edge_counter++;
3994 }
3995 }
3996 else
3997 {
3998 CGAL_triangulation_precondition(&*vh2 < &*vh1);
3999 if (edge_is_too_long(p2, p1) &&
4000 (find(_too_long_edges[vh2].begin(), _too_long_edges[vh2].end(), vh1) == _too_long_edges[vh2].end()))
4001 {
4002 _too_long_edges[vh2].push_back(vh1);
4003 _too_long_edge_counter++;
4004 }
4005 }
4006 }
4007
4008 template<class Gt, class Tds>
remove_too_long_edges_in_star(Vertex_handle vh)4009 void Periodic_2_triangulation_2<Gt, Tds>::remove_too_long_edges_in_star(
4010 Vertex_handle vh)
4011 {
4012 if (is_1_cover())
4013 return;
4014
4015 // Insert the too long edges in the star of vh
4016 Face_handle f = vh->face();
4017 Face_handle f_start = f;
4018
4019 do
4020 {
4021 int i = f->index(vh);
4022 int i2 = ccw(i);
4023 Vertex_handle vh2 = f->vertex(i2);
4024
4025 // Point corresponding to v
4026 Point p1 = construct_point(vh->point(), get_offset(f, f->index(vh)));
4027 // Point corresponding to the other vertex
4028 Point p2 = construct_point(vh2->point(), get_offset(f, i2));
4029
4030 if (&*vh < &*vh2)
4031 {
4032 if (edge_is_too_long(p1, p2) &&
4033 (find(_too_long_edges[vh].begin(), _too_long_edges[vh].end(), vh2) !=
4034 _too_long_edges[vh].end()))
4035 {
4036 _too_long_edges[vh].remove(vh2);
4037 _too_long_edge_counter--;
4038 }
4039 }
4040 else
4041 {
4042 CGAL_triangulation_precondition(&*vh2 < &*vh);
4043 if (edge_is_too_long(p1, p2) &&
4044 (find(_too_long_edges[vh2].begin(), _too_long_edges[vh2].end(), vh) !=
4045 _too_long_edges[vh2].end()))
4046 {
4047 _too_long_edges[vh2].remove(vh);
4048 _too_long_edge_counter--;
4049 }
4050 }
4051
4052 // Proceed to the next face
4053 f = f->neighbor(i2);
4054 }
4055 while (f != f_start);
4056 }
4057
4058 template<class Gt, class Tds>
remove_too_long_edge(Face_handle f,int i)4059 void Periodic_2_triangulation_2<Gt, Tds>::remove_too_long_edge(Face_handle f,
4060 int i)
4061 {
4062 Vertex_handle vh1 = f->vertex(cw(i));
4063 Vertex_handle vh2 = f->vertex(ccw(i));
4064 Point p1 = construct_point(vh1->point(), get_offset(f, cw(i)));
4065 Point p2 = construct_point(vh2->point(), get_offset(f, ccw(i)));
4066 if (edge_is_too_long(p1, p2))
4067 {
4068 if (&*vh1 < &*vh2)
4069 {
4070 typename std::list<Vertex_handle>::iterator it = find(
4071 _too_long_edges[vh1].begin(), _too_long_edges[vh1].end(), vh2);
4072 if (it != _too_long_edges[vh1].end())
4073 {
4074 _too_long_edges[vh1].erase(it);
4075 _too_long_edge_counter--;
4076 }
4077 }
4078 else
4079 {
4080 typename std::list<Vertex_handle>::iterator it = find(
4081 _too_long_edges[vh2].begin(), _too_long_edges[vh2].end(), vh1);
4082 if (it != _too_long_edges[vh2].end())
4083 {
4084 _too_long_edges[vh2].erase(it);
4085 _too_long_edge_counter--;
4086 }
4087 }
4088 }
4089 }
4090
4091 template<class Gt, class Tds>
edge_is_too_long(const Point & p1,const Point & p2)4092 bool Periodic_2_triangulation_2<Gt, Tds>::edge_is_too_long(const Point &p1,
4093 const Point &p2) const
4094 {
4095 return squared_distance(p1, p2) > _edge_length_threshold;
4096 }
4097
4098 template<class GT, class Tds>
is_triangulation_in_1_sheet()4099 inline bool Periodic_2_triangulation_2<GT, Tds>::is_triangulation_in_1_sheet() const
4100 {
4101 if (is_1_cover())
4102 return true;
4103 for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit)
4104 {
4105 if (_virtual_vertices.find(vit) == _virtual_vertices.end())
4106 continue;
4107 std::set<Vertex_handle> nb_v_odom;
4108 Vertex_handle vh;
4109 Offset off;
4110 Vertex_circulator vcir = adjacent_vertices(vit);
4111 Vertex_circulator vstart = vcir;
4112 size_t degree = 0;
4113 do
4114 {
4115 get_vertex(vcir, vh, off);
4116 nb_v_odom.insert(vh);
4117 degree++;
4118 }
4119 while (++vcir != vstart);
4120 if (degree != nb_v_odom.size())
4121 return false;
4122 }
4123 return true;
4124 }
4125
4126 template<class Gt, class Tds>
4127 std::ostream&
save(std::ostream & os)4128 Periodic_2_triangulation_2<Gt, Tds>::save(std::ostream& os) const
4129 {
4130 // writes :
4131 // the number of vertices
4132 // the domain as four coordinates: xmin ymin ymax zmax
4133 // the current covering that guarantees the triangulation to be a
4134 // simplicial complex
4135 // the non combinatorial information on vertices (points in case of 1-sheeted
4136 // covering, point-offset pairs otherwise)
4137 // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY
4138 // the number of faces
4139 // the faces by the indices of their vertices in the preceding list
4140 // of vertices, plus the non combinatorial information on each face
4141 // the neighbors of each face by their index in the preceding list of faces
4142
4143 // outputs dimension, domain and number of vertices
4144 Covering_sheets cover = number_of_sheets();
4145 size_type n = number_of_vertices();
4146
4147
4148 if (IO::is_ascii(os))
4149 os << domain() << std::endl
4150 << cover[0] << " " << cover[1] << std::endl
4151 << n*cover[0]*cover[1] << std::endl;
4152 else
4153 {
4154 os << domain();
4155 write(os, cover[0]);
4156 write(os, cover[1]);
4157 write(os, n * cover[0]*cover[1]);
4158 }
4159 std::cout << "Line:" << __LINE__ << " cover[0]:" << cover[0] << " cover[1]:" << cover[1] << " n*c0*c1:" << (n * cover[0]*cover[1]) << std::endl;
4160
4161 std::cout << "save, #Vertices: " << n << std::endl;
4162
4163 if (n == 0)
4164 return os;
4165
4166 // write the vertices
4167 Unique_hash_map<Vertex_handle, std::size_t > V;
4168 std::size_t i = 0;
4169 if (is_1_cover())
4170 {
4171 for (Vertex_iterator it = vertices_begin(); it != vertices_end(); ++it)
4172 {
4173 V[it] = i++;
4174 os << it->point();
4175 if (IO::is_ascii(os))
4176 os << std::endl;
4177 }
4178 }
4179 else
4180 {
4181 Virtual_vertex_map_it vit, vvit;
4182 std::vector<Vertex_handle> vv;
4183 for (Vertex_iterator it = vertices_begin(); it != vertices_end(); ++it)
4184 {
4185 vit = _virtual_vertices.find(it);
4186 if (vit != _virtual_vertices.end()) continue;
4187 V[it] = i++;
4188 if (IO::is_ascii(os))
4189 os << it->point() << std::endl
4190 << Offset(0, 0) << std::endl;
4191 else
4192 os << it->point() << Offset(0, 0);
4193 CGAL_triangulation_assertion(_virtual_vertices_reverse.find(it)
4194 != _virtual_vertices_reverse.end());
4195 vv = _virtual_vertices_reverse.find(it)->second;
4196 CGAL_triangulation_assertion(vv.size() == 8);
4197 for (std::size_t j = 0; j < vv.size(); j++)
4198 {
4199 vvit = _virtual_vertices.find(vv[j]);
4200 CGAL_triangulation_assertion(vvit != _virtual_vertices.end());
4201 V[vv[j]] = i++;
4202 if (IO::is_ascii(os))
4203 os << vv[j]->point() << std::endl
4204 << vvit->second.second << std::endl;
4205 else os << vv[j]->point() << vvit->second.second;
4206 }
4207 }
4208 }
4209 CGAL_triangulation_postcondition(i == _cover[0]*_cover[1]*n);
4210
4211 Unique_hash_map<Face_handle, std::size_t> F;
4212 int inum = 0;
4213 // asks the tds for the combinatorial information
4214 // vertices of the faces
4215 size_type m = _tds.number_of_faces();
4216 if (IO::is_ascii(os)) os << std::endl << m << std::endl;
4217 else write(os, m);
4218 std::cout << "save, #Faces: " << m << std::endl;
4219
4220 for( Face_iterator ib = faces_begin();
4221 ib != faces_end(); ++ib)
4222 {
4223 F[ib] = inum++;
4224 for(int j = 0; j < 3 ; ++j)
4225 {
4226 if(IO::is_ascii(os)) os << V[ib->vertex(j)] << " ";
4227 else write(os, V[ib->vertex(j)]);
4228 }
4229 os << *ib ;
4230 if(IO::is_ascii(os)) os << "\n";
4231 }
4232 if(IO::is_ascii(os)) os << "\n";
4233
4234 std::cout << "save, face check: " << inum << " == " << m << std::endl;
4235 CGAL_assertion(m == (size_type)inum);
4236
4237 // neighbor pointers of the faces
4238 for( Face_iterator it = faces_begin();
4239 it != faces_end(); ++it)
4240 {
4241 for(int j = 0; j < 3; ++j)
4242 {
4243 CGAL_assertion(F.is_defined(it->neighbor(j)));
4244 if(IO::is_ascii(os)) os << F[it->neighbor(j)] << " ";
4245 else write(os, F[it->neighbor(j)]);
4246 }
4247 if(IO::is_ascii(os)) os << "\n";
4248 }
4249
4250 // write offsets
4251 //for (unsigned int i=0 ; i<number_of_faces() ; i++) {
4252 for (Face_iterator it = faces_begin(); it != faces_end(); ++it)
4253 {
4254 //Face_handle ch = std::find(faces_begin(), faces_end(), i);
4255 Face_handle ch(it);
4256 for (int j = 0; j < 3; j++)
4257 {
4258 if(IO::is_ascii(os))
4259 {
4260 os << ch->offset(j);
4261 if ( j == 3 )
4262 os << std::endl;
4263 else
4264 os << ' ';
4265 }
4266 else write(os, ch->offset(j));
4267 }
4268 }
4269
4270 // write the non combinatorial information on the faces
4271 // using the << operator of Face
4272 // works because the iterator of the tds traverses the faces in the
4273 // same order as the iterator of the triangulation
4274 if(number_of_vertices() != 0)
4275 {
4276 for(Face_iterator it = faces_begin(); it != faces_end(); ++it)
4277 {
4278 os << *it; // other information
4279 if(IO::is_ascii(os))
4280 os << std::endl;
4281 }
4282 }
4283
4284 return os;
4285 }
4286
4287 template<class Gt, class Tds>
4288 std::istream&
load(std::istream & is)4289 Periodic_2_triangulation_2<Gt, Tds>::load(std::istream& is)
4290 {
4291 // reads
4292 // the current covering that guarantees the triangulation to be a
4293 // simplicial complex
4294 // the number of vertices
4295 // the non combinatorial information on vertices (points in case of 1-sheeted
4296 // covering, point-offset pairs otherwise)
4297 // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY
4298 // the number of faces
4299 // the faces by the indices of their vertices in the preceding list
4300 // of vertices, plus the non combinatorial information on each face
4301 // the neighbors of each face by their index in the preceding list of face
4302 CGAL_triangulation_precondition(is.good());
4303
4304 clear();
4305
4306 Iso_rectangle domain(0, 0, 1, 1);
4307 int cx = 0, cy = 0;
4308 size_type n = 0;
4309
4310 if (IO::is_ascii(is))
4311 {
4312 is >> domain;
4313 is >> cx >> cy >> n;
4314 }
4315 else
4316 {
4317 is >> domain;
4318 read(is, cx);
4319 read(is, cy);
4320 read(is, n);
4321 }
4322 std::cout << "Line:" << __LINE__ << " cx:" << cx << " cy:" << cy << " n:" << n << std::endl;
4323
4324 CGAL_triangulation_assertion((n / (cx * cy))*cx*cy == n);
4325
4326 tds().set_dimension((n == 0 ? -2 : 2));
4327 _domain = domain;
4328 _gt.set_domain(domain);
4329 _cover = make_array(cx, cy);
4330
4331 if ( n == 0 ) return is;
4332
4333 std::map< std::size_t, Vertex_handle > V;
4334
4335 if (cx == 1 && cy == 1)
4336 {
4337 Point p;
4338 for (std::size_t i = 0; i < n; i++)
4339 {
4340 V[i] = tds().create_vertex();
4341 is >> p;
4342 V[i]->set_point(p);
4343 }
4344 }
4345 else
4346 {
4347 Vertex_handle v, w;
4348 std::vector<Vertex_handle> vv;
4349 Offset off;
4350 Point p;
4351 for (std::size_t i = 0; i < n; i++)
4352 {
4353 v = tds().create_vertex();
4354 V[i] = v;
4355 is >> p >> off;
4356 V[i]->set_point(p);
4357 vv.clear();
4358 for (int j = 1; j < cx * cy; j++)
4359 {
4360 i++;
4361 w = tds().create_vertex();
4362 V[i] = w;
4363 is >> p >> off;
4364 V[i]->set_point(p);
4365 vv.push_back(w);
4366 _virtual_vertices[w] = std::make_pair(v, off);
4367 }
4368 _virtual_vertices_reverse[v] = vv;
4369 }
4370 }
4371
4372 // Creation of the faces
4373 std::size_t index;
4374 size_type m;
4375 if (IO::is_ascii(is)) is >> m;
4376 else read(is, m);
4377 std::vector<Face_handle> F(m);
4378 std::cout << "load, #Faces: " << m << std::endl;
4379 {
4380 for(size_t i = 0; i < m; ++i)
4381 {
4382 F[i] = _tds.create_face() ;
4383 for(int j = 0; j < 3 ; ++j)
4384 {
4385 if (IO::is_ascii(is)) is >> index;
4386 else read(is, index);
4387 CGAL_assertion(index < V.size());
4388 F[i]->set_vertex(j, V[index]);
4389 // The face pointer of vertices is set too often,
4390 // but otherwise we had to use one more map
4391 V[index]->set_face(F[i]);
4392 }
4393 // read in non combinatorial info of the face
4394 is >> *(F[i]) ;
4395 }
4396 }
4397
4398 // Setting the neighbor pointers
4399 {
4400 for(size_t i = 0; i < m; ++i)
4401 {
4402 for(int j = 0; j < 3; ++j)
4403 {
4404 if (IO::is_ascii(is)) is >> index;
4405 else read(is, index);
4406 if (index >= F.size()) {
4407 std::cout << __FILE__ << ", " << __FUNCTION__ << ", l:" << __LINE__ << " f="
4408 << i << "<" << m << ", index=" << j << " nb=" << index << " #F=" << F.size()
4409 << std::endl;
4410 }
4411 CGAL_assertion(i < F.size());
4412 CGAL_assertion(index < F.size());
4413 F[i]->set_neighbor(j, F[index]);
4414 }
4415 }
4416 }
4417
4418 // read offsets
4419 int off[3] = {0, 0, 0};
4420 for (std::size_t j = 0 ; j < m; j++)
4421 {
4422 if (IO::is_ascii(is))
4423 is >> off[0] >> off[1] >> off[2];
4424 else
4425 {
4426 read(is, off[0]);
4427 read(is, off[1]);
4428 read(is, off[2]);
4429 }
4430 set_offsets(F[j], off[0], off[1], off[2]);
4431 }
4432
4433 // read potential other information
4434 for (std::size_t j = 0 ; j < m; j++)
4435 is >> *(F[j]);
4436
4437 int i = 0;
4438 for (Vertex_iterator vi = vertices_begin();
4439 vi != vertices_end(); ++vi)
4440 {
4441 _too_long_edges[vi] = std::list<Vertex_handle>();
4442 ++i;
4443 }
4444
4445 _edge_length_threshold = FT(0.166) * (_domain.xmax() - _domain.xmin())
4446 * (_domain.xmax() - _domain.xmin());
4447 _too_long_edge_counter = find_too_long_edges(_too_long_edges);
4448
4449 CGAL_triangulation_expensive_assertion( is_valid() );
4450 return is;
4451 }
4452
4453
4454 namespace internal
4455 {
4456
4457 /// Internal function used by operator==.
4458 //TODO: introduce offsets
4459 template <class GT, class Tds1, class Tds2>
4460 bool
test_next(const Periodic_2_triangulation_2<GT,Tds1> & t1,const Periodic_2_triangulation_2<GT,Tds2> & t2,typename Periodic_2_triangulation_2<GT,Tds1>::Face_handle c1,typename Periodic_2_triangulation_2<GT,Tds2>::Face_handle c2,std::map<typename Periodic_2_triangulation_2<GT,Tds1>::Face_handle,typename Periodic_2_triangulation_2<GT,Tds2>::Face_handle> & Cmap,std::map<typename Periodic_2_triangulation_2<GT,Tds1>::Vertex_handle,typename Periodic_2_triangulation_2<GT,Tds2>::Vertex_handle> & Vmap)4461 test_next(const Periodic_2_triangulation_2<GT, Tds1> &t1,
4462 const Periodic_2_triangulation_2<GT, Tds2> &t2,
4463 typename Periodic_2_triangulation_2<GT, Tds1>::Face_handle c1,
4464 typename Periodic_2_triangulation_2<GT, Tds2>::Face_handle c2,
4465 std::map < typename Periodic_2_triangulation_2<GT, Tds1>::Face_handle,
4466 typename Periodic_2_triangulation_2<GT, Tds2>::Face_handle > &Cmap,
4467 std::map < typename Periodic_2_triangulation_2<GT, Tds1>::Vertex_handle,
4468 typename Periodic_2_triangulation_2<GT, Tds2>::Vertex_handle > &Vmap)
4469 {
4470 // This function tests and registers the 4 neighbors of c1/c2,
4471 // and recursively calls itself over them.
4472 // Returns false if an inequality has been found.
4473
4474 // Precondition: c1, c2 have been registered as well as their 4 vertices.
4475 CGAL_triangulation_precondition(t1.number_of_vertices() != 0);
4476 CGAL_triangulation_precondition(Cmap[c1] == c2);
4477 CGAL_triangulation_precondition(Vmap.find(c1->vertex(0)) != Vmap.end());
4478 CGAL_triangulation_precondition(Vmap.find(c1->vertex(1)) != Vmap.end());
4479 CGAL_triangulation_precondition(Vmap.find(c1->vertex(2)) != Vmap.end());
4480
4481 typedef Periodic_2_triangulation_2<GT, Tds1> Tr1;
4482 typedef Periodic_2_triangulation_2<GT, Tds2> Tr2;
4483 typedef typename Tr1::Vertex_handle Vertex_handle1;
4484 typedef typename Tr1::Face_handle Face_handle1;
4485 typedef typename Tr2::Vertex_handle Vertex_handle2;
4486 typedef typename Tr2::Face_handle Face_handle2;
4487 typedef typename std::map<Face_handle1, Face_handle2>::const_iterator Cit;
4488 typedef typename std::map < Vertex_handle1,
4489 Vertex_handle2 >::const_iterator Vit;
4490
4491 for (int i = 0; i <= 2; ++i)
4492 {
4493 Face_handle1 n1 = c1->neighbor(i);
4494 Cit cit = Cmap.find(n1);
4495 Vertex_handle1 v1 = c1->vertex(i);
4496 Vertex_handle2 v2 = Vmap[v1];
4497 Face_handle2 n2 = c2->neighbor(c2->index(v2));
4498 if (cit != Cmap.end())
4499 {
4500 // n1 was already registered.
4501 if (cit->second != n2)
4502 return false;
4503 continue;
4504 }
4505 // n1 has not yet been registered.
4506 // We check that the new vertices match geometrically.
4507 // And we register them.
4508 Vertex_handle1 vn1 = n1->vertex(n1->index(c1));
4509 Vertex_handle2 vn2 = n2->vertex(n2->index(c2));
4510 Vit vit = Vmap.find(vn1);
4511 if (vit != Vmap.end())
4512 {
4513 // vn1 already registered
4514 if (vit->second != vn2)
4515 return false;
4516 }
4517 else
4518 {
4519 if (t1.geom_traits().compare_xy_2_object()(vn1->point(),
4520 vn2->point()) != 0)
4521 return false;
4522
4523 // We register vn1/vn2.
4524 Vmap.insert(std::make_pair(vn1, vn2));
4525 }
4526
4527 // We register n1/n2.
4528 Cmap.insert(std::make_pair(n1, n2));
4529 // We recurse on n1/n2.
4530 if (!test_next(t1, t2, n1, n2, Cmap, Vmap))
4531 return false;
4532 }
4533
4534 return true;
4535 }
4536
4537 } // namespace internal
4538
4539
4540 template<class Gt, class Tds>
4541 std::istream&
4542 operator>>(std::istream& is, Periodic_2_triangulation_2<Gt, Tds> &tr)
4543 {
4544 return tr.load(is);
4545 }
4546 template<class Gt, class Tds>
4547 std::ostream&
4548 operator<<(std::ostream& os, Periodic_2_triangulation_2<Gt, Tds> &tr)
4549 {
4550 return tr.save(os);
4551 }
4552
4553 template < class GT, class Tds1, class Tds2 >
4554 bool
4555 operator==(const Periodic_2_triangulation_2<GT, Tds1> &t1,
4556 const Periodic_2_triangulation_2<GT, Tds2> &t2)
4557 {
4558 typedef typename Periodic_2_triangulation_2<GT, Tds1>::Vertex_handle
4559 Vertex_handle1;
4560 typedef typename Periodic_2_triangulation_2<GT, Tds1>::Face_handle
4561 Face_handle1;
4562 typedef typename Periodic_2_triangulation_2<GT, Tds2>::Vertex_handle
4563 Vertex_handle2;
4564 typedef typename Periodic_2_triangulation_2<GT, Tds2>::Vertex_handle
4565 Vertex_iterator2;
4566 typedef typename Periodic_2_triangulation_2<GT, Tds2>::Face_handle
4567 Face_handle2;
4568 typedef typename Periodic_2_triangulation_2<GT, Tds2>::Face_circulator
4569 Face_circulator2;
4570
4571 typedef typename Periodic_2_triangulation_2<GT, Tds1>::Point Point;
4572 typedef typename Periodic_2_triangulation_2<GT, Tds1>::Offset Offset;
4573
4574 // Some quick checks.
4575 if ( t1.domain() != t2.domain()
4576 || t1.number_of_sheets() != t2.number_of_sheets())
4577 return false;
4578
4579 if ( t1.number_of_vertices() != t2.number_of_vertices()
4580 || t1.number_of_faces() != t2.number_of_faces())
4581 return false;
4582
4583 // Special case for empty triangulations
4584 if (t1.number_of_vertices() == 0)
4585 return true;
4586
4587 // We will store the mapping between the 2 triangulations vertices and
4588 // faces in 2 maps.
4589 std::map<Vertex_handle1, Vertex_handle2> Vmap;
4590 std::map<Face_handle1, Face_handle2> Cmap;
4591
4592 // find a common point
4593 Vertex_handle1 v1 = static_cast<Vertex_handle1>(t1.vertices_begin());
4594 Vertex_handle2 iv2;
4595 for (Vertex_iterator2 vit2 = t2.vertices_begin() ;
4596 vit2 != t2.vertices_end(); ++vit2)
4597 {
4598 if (t1.compare_xy(vit2->point(), v1->point(),
4599 t2.get_offset(vit2), t1.get_offset(v1)) != EQUAL)
4600 continue;
4601 iv2 = static_cast<Vertex_handle2>(vit2);
4602 break;
4603 }
4604 if (iv2 == Vertex_handle2())
4605 return false;
4606 Vmap.insert(std::make_pair(v1, iv2));
4607
4608 // We pick one face of t1, and try to match it against the
4609 // faces of t2.
4610 Face_handle1 c = v1->face();
4611 Vertex_handle1 v2 = c->vertex(t1.cw(c->index(v1)));
4612 Vertex_handle1 v3 = c->vertex(t1.ccw(c->index(v1)));
4613 Point p2 = v2->point();
4614 Point p3 = v3->point();
4615 Offset o2 = t1.get_offset(v2);
4616 Offset o3 = t1.get_offset(v3);
4617
4618 Face_circulator2 fc = t2.incident_faces(iv2), done(fc);
4619 do
4620 {
4621 int inf = fc->index(iv2);
4622
4623 if (t1.compare_xy(p2, fc->vertex((inf + 1) % 3)->point(),
4624 o2, t2.get_offset(fc->vertex((inf + 1) % 3))) == EQUAL)
4625 Vmap.insert(std::make_pair(v2, fc->vertex((inf + 1) % 3)));
4626 else if (t1.compare_xy(p2, fc->vertex((inf + 2) % 3)->point(),
4627 o2, t2.get_offset(fc->vertex((inf + 2) % 3))) == EQUAL)
4628 Vmap.insert(std::make_pair(v2, fc->vertex((inf + 2) % 3)));
4629 else
4630 continue; // None matched v2.
4631
4632 if (t1.compare_xy(p3, fc->vertex((inf + 1) % 3)->point(),
4633 o3, t2.get_offset(fc->vertex((inf + 1) % 3))) == EQUAL)
4634 Vmap.insert(std::make_pair(v3, fc->vertex((inf + 1) % 3)));
4635 else if (t1.compare_xy(p3, fc->vertex((inf + 2) % 3)->point(),
4636 o3, t2.get_offset(fc->vertex((inf + 2) % 3))) == EQUAL)
4637 Vmap.insert(std::make_pair(v3, fc->vertex((inf + 2) % 3)));
4638 else
4639 continue; // None matched v3.
4640
4641 // Found it !
4642 Cmap.insert(std::make_pair(c, fc));
4643 break;
4644 }
4645 while (++fc != done);
4646
4647 if (Cmap.size() == 0)
4648 return false;
4649
4650 // We now have one face, we need to propagate recursively.
4651 return internal::test_next(t1, t2,
4652 Cmap.begin()->first, Cmap.begin()->second, Cmap, Vmap);
4653 }
4654
4655 template < class GT, class Tds1, class Tds2 >
4656 inline
4657 bool
4658 operator!=(const Periodic_2_triangulation_2<GT, Tds1> &t1,
4659 const Periodic_2_triangulation_2<GT, Tds2> &t2)
4660 {
4661 return ! (t1 == t2);
4662 }
4663
4664 #define CGAL_INCLUDE_FROM_PERIODIC_2_TRIANGULATION_2_H
4665 #include <CGAL/Periodic_2_triangulation_dummy_12.h>
4666 #undef CGAL_INCLUDE_FROM_PERIODIC_2_TRIANGULATION_2_H
4667 } //namespace CGAL
4668
4669
4670 #endif //CGAL_PERIODIC_2_TRIANGULATION_2_H
4671