1 // Copyright (c) 1999-2003,2006-2009,2014-2015,2017 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_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h $ 7 // $Id: Periodic_3_triangulation_3.h d6b2c8d 2021-05-18T18:13:38+02:00 Laurent Rineau 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s) : Monique Teillaud <Monique.Teillaud@sophia.inria.fr> 11 // Sylvain Pion <Sylvain.Pion@sophia.inria.fr> 12 // Nico Kruithof <Nico.Kruithof@sophia.inria.fr> 13 // Manuel Caroli <Manuel.Caroli@sophia.inria.fr> 14 // Aymeric Pelle <Aymeric.Pelle@sophia.inria.fr> 15 // Mael Rouxel-Labbé 16 17 #ifndef CGAL_PERIODIC_3_TRIANGULATION_3_H 18 #define CGAL_PERIODIC_3_TRIANGULATION_3_H 19 20 #include <CGAL/license/Periodic_3_triangulation_3.h> 21 22 #include <CGAL/basic.h> 23 24 #include <CGAL/internal/Periodic_3_triangulation_iterators_3.h> 25 #include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h> 26 #include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h> 27 #include <CGAL/Periodic_3_triangulation_traits_3.h> 28 #include <CGAL/Triangulation_data_structure_3.h> 29 #include <CGAL/Triangulation_cell_base_3.h> 30 #include <CGAL/Triangulation_vertex_base_3.h> 31 #include <CGAL/triangulation_assertions.h> 32 #include <CGAL/internal/canonicalize_helper.h> 33 34 #include <CGAL/array.h> 35 #include <CGAL/internal/Exact_type_selector.h> 36 #include <CGAL/NT_converter.h> 37 #include <CGAL/Unique_hash_map.h> 38 #include <CGAL/use.h> 39 40 #ifndef CGAL_NO_STRUCTURAL_FILTERING 41 #include <CGAL/internal/Static_filters/tools.h> 42 #include <CGAL/Triangulation_structural_filtering_traits.h> 43 #include <CGAL/determinant.h> 44 #endif // no CGAL_NO_STRUCTURAL_FILTERING 45 46 #include <boost/random/linear_congruential.hpp> 47 #include <boost/random/uniform_smallint.hpp> 48 #include <boost/random/variate_generator.hpp> 49 #include <boost/tuple/tuple.hpp> 50 #include <boost/unordered_map.hpp> 51 52 #include <iostream> 53 #include <algorithm> 54 #include <cmath> 55 #include <functional> 56 #include <list> 57 #include <utility> 58 59 namespace CGAL { 60 61 template < class GT, class TDS > class Periodic_3_triangulation_3; 62 63 template < class GT, class TDS > std::istream& operator>> 64 (std::istream& is, Periodic_3_triangulation_3<GT,TDS>& tr); 65 template < class GT, class TDS > std::ostream& operator<< 66 (std::ostream& os, const Periodic_3_triangulation_3<GT,TDS>& tr); 67 68 #ifndef CGAL_NO_STRUCTURAL_FILTERING 69 namespace internal { 70 // structural filtering is performed only for EPIC 71 struct Periodic_structural_filtering_3_tag {}; 72 struct No_periodic_structural_filtering_3_tag {}; 73 74 template <bool filter> 75 struct Periodic_structural_filtering_selector_3 76 { 77 #ifdef FORCE_STRUCTURAL_FILTERING 78 typedef Periodic_structural_filtering_3_tag Tag; 79 #else 80 typedef No_periodic_structural_filtering_3_tag Tag; 81 #endif 82 }; 83 84 template <> 85 struct Periodic_structural_filtering_selector_3<true> 86 { 87 typedef Periodic_structural_filtering_3_tag Tag; 88 }; 89 } 90 #endif // no CGAL_NO_STRUCTURAL_FILTERING 91 92 /**\class Periodic_3_triangulation_3 93 * 94 * \brief implements functionality for computing in periodic space. 95 * 96 * There are several things that are special to computing in $\mathbb{T}^3$ 97 * such as 98 * - periodicity --> offsets 99 * - no infinite vertex 100 * - no degenerate dimensions 101 * All functions that are affected can be found in this class. When it is 102 * necessary to provide different implementations for Delaunay and regular 103 * triangulation, we work with visitors. 104 */ 105 106 template < class GT, 107 class TDS = Triangulation_data_structure_3 < 108 Triangulation_vertex_base_3<GT, 109 Periodic_3_triangulation_ds_vertex_base_3<> >, 110 Triangulation_cell_base_3<GT, 111 Periodic_3_triangulation_ds_cell_base_3<> > > > 112 class Periodic_3_triangulation_3 113 : public Triangulation_utils_3 114 { 115 friend std::istream& operator>> <> 116 (std::istream& is, Periodic_3_triangulation_3<GT, TDS>& tr); 117 118 typedef Periodic_3_triangulation_3<GT,TDS> Self; 119 120 public: 121 typedef GT Geometric_traits; 122 typedef TDS Triangulation_data_structure; 123 124 typedef typename GT::Periodic_3_offset_3 Offset; 125 typedef typename GT::Iso_cuboid_3 Iso_cuboid; 126 typedef std::array<int, 3> Covering_sheets; 127 128 // point types 129 typedef typename TDS::Vertex::Point Point; 130 typedef typename GT::Point_3 Point_3; 131 typedef std::pair<Point, Offset> Periodic_point; 132 typedef std::pair<Point_3, Offset> Periodic_point_3; 133 134 typedef typename GT::Segment_3 Segment; 135 typedef typename GT::Triangle_3 Triangle; 136 typedef typename GT::Tetrahedron_3 Tetrahedron; 137 138 typedef std::array<Periodic_point, 2> Periodic_segment; 139 typedef std::array<Periodic_point, 3> Periodic_triangle; 140 typedef std::array<Periodic_point, 4> Periodic_tetrahedron; 141 142 typedef std::array<Periodic_point_3, 2> Periodic_segment_3; 143 typedef std::array<Periodic_point_3, 3> Periodic_triangle_3; 144 typedef std::array<Periodic_point_3, 4> Periodic_tetrahedron_3; 145 146 typedef typename TDS::Vertex Vertex; 147 typedef typename TDS::Cell Cell; 148 typedef typename TDS::Facet Facet; 149 typedef typename TDS::Edge Edge; 150 151 typedef typename TDS::Vertex_handle Vertex_handle; 152 typedef typename TDS::Cell_handle Cell_handle; 153 154 typedef typename TDS::size_type size_type; 155 typedef typename TDS::difference_type difference_type; 156 157 typedef typename TDS::Cell_iterator Cell_iterator; 158 typedef typename TDS::Facet_iterator Facet_iterator; 159 typedef typename TDS::Edge_iterator Edge_iterator; 160 typedef typename TDS::Vertex_iterator Vertex_iterator; 161 162 typedef typename TDS::Cell_circulator Cell_circulator; 163 typedef typename TDS::Facet_circulator Facet_circulator; 164 165 typedef Cell_iterator All_cells_iterator; 166 typedef Facet_iterator All_facets_iterator; 167 typedef Edge_iterator All_edges_iterator; 168 typedef Vertex_iterator All_vertices_iterator; 169 typedef Periodic_3_triangulation_unique_vertex_iterator_3<Self> 170 Unique_vertex_iterator; 171 172 private: 173 typedef typename GT::FT FT; 174 typedef std::pair< Vertex_handle, Offset > Virtual_vertex; 175 typedef boost::unordered_map<Vertex_handle, Virtual_vertex> 176 Virtual_vertex_map; 177 typedef typename Virtual_vertex_map::const_iterator 178 Virtual_vertex_map_it; 179 typedef boost::unordered_map<Vertex_handle, std::vector<Vertex_handle > > 180 Virtual_vertex_reverse_map; 181 typedef typename Virtual_vertex_reverse_map::const_iterator 182 Virtual_vertex_reverse_map_it; 183 typedef Triple< Vertex_handle, Vertex_handle, Vertex_handle > 184 Vertex_triple; 185 186 public: 187 typedef Periodic_3_triangulation_tetrahedron_iterator_3<Self> 188 Periodic_tetrahedron_iterator; 189 typedef Periodic_3_triangulation_triangle_iterator_3<Self> 190 Periodic_triangle_iterator; 191 typedef Periodic_3_triangulation_segment_iterator_3<Self> 192 Periodic_segment_iterator; 193 typedef Periodic_3_triangulation_point_iterator_3<Self> 194 Periodic_point_iterator; 195 196 typedef Point value_type; 197 typedef const value_type& const_reference; 198 199 //Tag to distinguish regular triangulations from others 200 typedef Tag_false Weighted_tag; 201 202 // Tag to distinguish periodic triangulations from others 203 typedef Tag_true Periodic_tag; 204 205 public: 206 enum Iterator_type { 207 STORED = 0, 208 UNIQUE, //1 209 STORED_COVER_DOMAIN, //2 210 UNIQUE_COVER_DOMAIN };//3 211 212 enum Locate_type { 213 VERTEX = 0, 214 EDGE, //1 215 FACET, //2 216 CELL, //3 217 EMPTY , //4 218 OUTSIDE_CONVEX_HULL, // unused, for compatibility with Alpha_shape_3 219 OUTSIDE_AFFINE_HULL }; // unused, for compatibility with Alpha_shape_3 220 221 // unused and undocumented types and functions required to be 222 // compatible with Alpha_shape_3 / Periodic_3_mesh_3 223 public: 224 typedef Cell_iterator Finite_cells_iterator; 225 typedef Facet_iterator Finite_facets_iterator; 226 typedef Edge_iterator Finite_edges_iterator; 227 typedef Vertex_iterator Finite_vertices_iterator; 228 229 int dimension() const { return (number_of_vertices() == 0) ? -2 : 3; } 230 231 template<class T> 232 bool is_infinite(const T&, int = 0, int = 0) const { return false; } 233 234 size_type number_of_finite_cells() const { return number_of_cells(); } 235 size_type number_of_finite_facets() const { return number_of_facets(); } 236 size_type number_of_finite_edges() const { return number_of_edges(); } 237 size_type number_of_finite_vertices() const { return number_of_vertices(); } 238 239 private: 240 Geometric_traits _gt; 241 Triangulation_data_structure _tds; 242 243 protected: 244 /// map of offsets for periodic copies of vertices 245 Virtual_vertex_map virtual_vertices; 246 Virtual_vertex_reverse_map virtual_vertices_reverse; 247 248 protected: 249 /// v_offsets temporarily stores all the vertices on the border of a 250 /// conflict region. 251 mutable std::vector<Vertex_handle> v_offsets; 252 253 /// Determines if we currently compute in 3-cover or 1-cover. 254 Covering_sheets _cover; 255 256 public: 257 /** @name Creation */ 258 Periodic_3_triangulation_3(const Iso_cuboid& domain = Iso_cuboid(0,0,0, 1,1,1), 259 const Geometric_traits& gt = Geometric_traits()) 260 : _gt(gt), _tds() 261 { 262 typedef typename internal::Exact_field_selector<FT>::Type EFT; 263 typedef NT_converter<FT,EFT> NTC; 264 CGAL_USE_TYPE(NTC); 265 CGAL_triangulation_precondition_code( NTC ntc; ) 266 CGAL_triangulation_precondition(ntc(domain.xmax())-ntc(domain.xmin()) 267 == ntc(domain.ymax())-ntc(domain.ymin())); 268 CGAL_triangulation_precondition(ntc(domain.ymax())-ntc(domain.ymin()) 269 == ntc(domain.zmax())-ntc(domain.zmin())); 270 CGAL_triangulation_precondition(ntc(domain.zmax())-ntc(domain.zmin()) 271 == ntc(domain.xmax())-ntc(domain.xmin())); 272 273 _gt.set_domain(domain); 274 _cover = CGAL::make_array(3,3,3); 275 init_tds(); 276 } 277 278 protected: 279 // Copy constructor helpers 280 class Finder; 281 282 public: 283 // Copy constructor duplicates vertices and cells 284 Periodic_3_triangulation_3(const Periodic_3_triangulation_3& tr) 285 : _gt(tr.geom_traits()), 286 _cover(tr.number_of_sheets()) 287 { 288 if(is_1_cover()) 289 tds() = tr.tds(); 290 else 291 copy_multiple_covering(tr); 292 293 CGAL_triangulation_expensive_postcondition(*this == tr); 294 CGAL_triangulation_expensive_postcondition(is_valid()); 295 } 296 297 virtual ~Periodic_3_triangulation_3() {} 298 299 void copy_multiple_covering(const Periodic_3_triangulation_3& tr) 300 { 301 // Write the respective offsets in the vertices to make them 302 // automatically copy with the tds. 303 for(Vertex_iterator vit = tr.vertices_begin(); vit != tr.vertices_end(); ++vit) 304 vit->set_offset(tr.get_offset(vit)); 305 306 // copy the tds 307 tds() = tr.tds(); 308 309 // make a list of all vertices that belong to the original 310 // domain and initialize the basic structure of 311 // virtual_vertices_reverse 312 std::list<Vertex_handle> vlist; 313 for(Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) 314 { 315 if(vit->offset() == Offset()) 316 { 317 vlist.push_back(vit); 318 virtual_vertices_reverse.insert(std::make_pair(vit,std::vector<Vertex_handle>(26))); 319 CGAL_triangulation_assertion(virtual_vertices_reverse.find(vit)->second.size() == 26); 320 } 321 } 322 323 // Iterate over all vertices that are not in the original domain 324 // and construct the respective entries to virtual_vertices and 325 // virtual_vertices_reverse 326 for(Vertex_iterator vit2 = vertices_begin(); vit2 != vertices_end(); ++vit2) 327 { 328 if(vit2->offset() != Offset()) 329 { 330 typename std::list<Vertex_handle>::iterator vlist_it 331 = std::find_if(vlist.begin(), vlist.end(), Finder(this,vit2->point())); 332 Offset off = vit2->offset(); 333 virtual_vertices.insert(std::make_pair(vit2, std::make_pair(*vlist_it,off))); 334 virtual_vertices_reverse.find(*vlist_it)->second[9*off[0]+3*off[1]+off[2]-1] = vit2; 335 CGAL_triangulation_assertion(get_offset(vit2) == off); 336 } 337 } 338 339 // Cleanup vertex offsets 340 for(Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) 341 vit->clear_offset(); 342 for(Vertex_iterator vit = tr.vertices_begin(); vit != tr.vertices_end(); ++vit) 343 vit->clear_offset(); 344 } 345 346 /** @name Assignment */ 347 Periodic_3_triangulation_3& operator=(Periodic_3_triangulation_3 tr) 348 { 349 swap(tr); 350 return *this; 351 } 352 353 void swap(Periodic_3_triangulation_3& tr) 354 { 355 std::swap(tr._gt, _gt); 356 _tds.swap(tr._tds); 357 std::swap(virtual_vertices,tr.virtual_vertices); 358 std::swap(virtual_vertices_reverse,tr.virtual_vertices_reverse); 359 std::swap(_cover, tr._cover); 360 } 361 362 virtual void clear_covering_data() { } 363 364 /// Clears the triangulation and initializes it again. 365 void clear() 366 { 367 _tds.clear(); 368 clear_covering_data(); 369 virtual_vertices.clear(); 370 virtual_vertices_reverse.clear(); 371 v_offsets.clear(); 372 373 _cover = CGAL::make_array(3,3,3); 374 init_tds(); 375 } 376 377 private: 378 /// Initializes the triangulation data structure 379 void init_tds() 380 { 381 _tds.set_dimension(-2); 382 v_offsets.reserve(48); 383 } 384 385 public: 386 /** @name Access functions */ 387 const Geometric_traits& geom_traits() const { return _gt; } 388 const TDS& tds() const { return _tds; } 389 TDS& tds() { return _tds; } 390 391 bool is_parallel() const { return false; } 392 393 virtual void gather_cell_hidden_points(const Cell_handle /*cit*/, std::vector<Point>& /* hidden_points*/) { } 394 virtual void reinsert_hidden_points_after_converting_to_1_sheeted(const std::vector<Point>& /* hidden_points*/) { } 395 396 const Iso_cuboid& domain() const { return _gt.get_domain(); } 397 398 virtual void update_cover_data_after_setting_domain() {} 399 400 void set_domain(const Iso_cuboid& domain) 401 { 402 clear(); 403 _gt.set_domain(domain); 404 update_cover_data_after_setting_domain(); 405 } 406 407 const Covering_sheets& number_of_sheets() const { return _cover; } 408 const std::pair<Vertex_handle, Offset> original_vertex(const Vertex_handle v) const 409 { 410 return (virtual_vertices.find(v) == virtual_vertices.end()) ? 411 std::make_pair(v,Offset()) : virtual_vertices.find(v)->second; 412 } 413 const std::vector<Vertex_handle>& periodic_copies(const Vertex_handle v) const 414 { 415 CGAL_triangulation_precondition(number_of_sheets() != CGAL::make_array(1,1,1)); 416 CGAL_triangulation_precondition(virtual_vertices.find(v) == virtual_vertices.end()); 417 CGAL_triangulation_assertion( 418 virtual_vertices_reverse.find(v) != virtual_vertices_reverse.end()); 419 return virtual_vertices_reverse.find(v)->second; 420 } 421 422 bool is_triangulation_in_1_sheet() const; 423 424 void convert_to_1_sheeted_covering(); 425 virtual void update_cover_data_after_converting_to_27_sheeted_covering() { } 426 void convert_to_27_sheeted_covering(); 427 428 size_type number_of_cells() const { 429 if(is_1_cover()) return _tds.number_of_cells(); 430 else return _tds.number_of_cells()/27; 431 } 432 size_type number_of_facets() const { 433 if(is_1_cover()) return _tds.number_of_facets(); 434 else return _tds.number_of_facets()/27; 435 } 436 size_type number_of_edges() const { 437 if(is_1_cover()) return _tds.number_of_edges(); 438 else return _tds.number_of_edges()/27; 439 } 440 size_type number_of_vertices() const { 441 if(is_1_cover()) return _tds.number_of_vertices(); 442 else return _tds.number_of_vertices()/27; 443 } 444 445 size_type number_of_stored_cells() const { 446 return _tds.number_of_cells(); 447 } 448 size_type number_of_stored_facets() const { 449 return _tds.number_of_facets(); 450 } 451 size_type number_of_stored_edges() const { 452 return _tds.number_of_edges(); 453 } 454 size_type number_of_stored_vertices() const { 455 return _tds.number_of_vertices(); 456 } 457 458 public: 459 bool is_1_cover() const 460 { 461 bool flag; 462 flag = ((_cover[0] == 1) && (_cover[1] == 1) && (_cover[2] == 1)); 463 return flag; 464 } 465 466 void set_cover(const Covering_sheets& cover) 467 { 468 _cover = cover; 469 } 470 471 public: 472 bool is_virtual(Vertex_handle v) 473 { 474 if(is_1_cover()) 475 return false; 476 return (virtual_vertices.find(v) != virtual_vertices.end()); 477 } 478 479 public: 480 // Offset converters 481 int off_to_int(const Offset& off) const 482 { 483 CGAL_triangulation_assertion( off.x() == 0 || off.x() == 1 ); 484 CGAL_triangulation_assertion( off.y() == 0 || off.y() == 1 ); 485 CGAL_triangulation_assertion( off.z() == 0 || off.z() == 1 ); 486 int i = ((off.x()&1)<<2) + ((off.y()&1)<<1) + ((off.z()&1)); 487 return i; 488 } 489 490 Offset int_to_off(int i) const 491 { 492 return Offset((i>>2)&1,(i>>1)&1,i&1); 493 } 494 495 void set_offsets(Cell_handle c, int o0,int o1,int o2,int o3) 496 { 497 int off0[3] = {(o0>>2)&1,(o0>>1)&1,(o0&1)}; 498 int off1[3] = {(o1>>2)&1,(o1>>1)&1,(o1&1)}; 499 int off2[3] = {(o2>>2)&1,(o2>>1)&1,(o2&1)}; 500 int off3[3] = {(o3>>2)&1,(o3>>1)&1,(o3&1)}; 501 for(int i=0; i<3; i++) { 502 int min_off = (std::min)((std::min)(off0[i],off1[i]), 503 (std::min)(off2[i],off3[i])); 504 if(min_off != 0) { 505 off0[i] -= min_off; off1[i] -= min_off; 506 off2[i] -= min_off; off3[i] -= min_off; 507 } 508 } 509 o0 = ((off0[0]&1)<<2)+((off0[1]&1)<<1)+(off0[2]&1); 510 o1 = ((off1[0]&1)<<2)+((off1[1]&1)<<1)+(off1[2]&1); 511 o2 = ((off2[0]&1)<<2)+((off2[1]&1)<<1)+(off2[2]&1); 512 o3 = ((off3[0]&1)<<2)+((off3[1]&1)<<1)+(off3[2]&1); 513 c->set_offsets(o0,o1,o2,o3); 514 } 515 516 template <class Offset> 517 void set_offsets(Cell_handle c, Offset o0,Offset o1,Offset o2,Offset o3) 518 { 519 int off0[3] = {o0.x(),o0.y(),o0.z()}; 520 int off1[3] = {o1.x(),o1.y(),o1.z()}; 521 int off2[3] = {o2.x(),o2.y(),o2.z()}; 522 int off3[3] = {o3.x(),o3.y(),o3.z()}; 523 for(int i=0; i<3; i++) { 524 int min_off = (std::min)((std::min)(off0[i],off1[i]), 525 (std::min)(off2[i],off3[i])); 526 if(min_off != 0) { 527 off0[i] -= min_off; off1[i] -= min_off; 528 off2[i] -= min_off; off3[i] -= min_off; 529 } 530 } 531 532 CGAL_triangulation_assertion((std::min)((std::min)(off0[0],off1[0]), 533 (std::min)(off2[0],off3[0])) == 0); 534 CGAL_triangulation_assertion((std::min)((std::min)(off0[1],off1[1]), 535 (std::min)(off2[1],off3[1])) == 0); 536 CGAL_triangulation_assertion((std::min)((std::min)(off0[2],off1[2]), 537 (std::min)(off2[2],off3[2])) == 0); 538 CGAL_triangulation_assertion((0 <= off0[0]) && (off0[0] < 2)); 539 CGAL_triangulation_assertion((0 <= off1[0]) && (off1[0] < 2)); 540 CGAL_triangulation_assertion((0 <= off2[0]) && (off2[0] < 2)); 541 CGAL_triangulation_assertion((0 <= off3[0]) && (off3[0] < 2)); 542 CGAL_triangulation_assertion((0 <= off0[1]) && (off0[1] < 2)); 543 CGAL_triangulation_assertion((0 <= off1[1]) && (off1[1] < 2)); 544 CGAL_triangulation_assertion((0 <= off2[1]) && (off2[1] < 2)); 545 CGAL_triangulation_assertion((0 <= off3[1]) && (off3[1] < 2)); 546 CGAL_triangulation_assertion((0 <= off0[2]) && (off0[2] < 2)); 547 CGAL_triangulation_assertion((0 <= off1[2]) && (off1[2] < 2)); 548 CGAL_triangulation_assertion((0 <= off2[2]) && (off2[2] < 2)); 549 CGAL_triangulation_assertion((0 <= off3[2]) && (off3[2] < 2)); 550 551 int o0i = ((off0[0]&1)<<2)+((off0[1]&1)<<1)+(off0[2]&1); 552 int o1i = ((off1[0]&1)<<2)+((off1[1]&1)<<1)+(off1[2]&1); 553 int o2i = ((off2[0]&1)<<2)+((off2[1]&1)<<1)+(off2[2]&1); 554 int o3i = ((off3[0]&1)<<2)+((off3[1]&1)<<1)+(off3[2]&1); 555 c->set_offsets(o0i,o1i,o2i,o3i); 556 } 557 558 public: 559 /** @name Wrapping the traits */ 560 // Note that calling functors with "construct_point(p), offset" and not 561 // construct_point(p, offset) is done on purpose. Indeed, construct_point(p) 562 // is not a real construction: it's either the identity (when `Point` is `Point_3`) 563 // or getting the bare point within a weighted point. However, construct_point(p, offset) 564 // is a real construction. When using filters, construct_point() must be done 565 // within the filtered predicates with the appropriate construct_point_3_object. 566 567 template<typename P> // can be Point or Point_3 568 Comparison_result compare_xyz(const P& p1, const P& p2) const { 569 return geom_traits().compare_xyz_3_object()(construct_point(p1), 570 construct_point(p2)); 571 } 572 template<typename P> // can be Point or Point_3 573 Comparison_result compare_xyz(const P& p1, const P& p2, 574 const Offset& o1, const Offset& o2) const { 575 return geom_traits().compare_xyz_3_object()( 576 construct_point(p1), construct_point(p2), o1, o2); 577 } 578 template<typename P> // can be Point or Point_3 579 Comparison_result compare_xyz(const std::pair<P, Offset>& pp1, 580 const std::pair<P, Offset>& pp2) const { 581 return geom_traits().compare_xyz_3_object()( 582 construct_point(pp1.first), construct_point(pp2.first), pp1.second, pp2.second); 583 } 584 585 template<typename P> // can be Point or Point_3 586 Orientation orientation(const P& p1, const P& p2, const P& p3, const P& p4) const { 587 return geom_traits().orientation_3_object()(construct_point(p1), construct_point(p2), 588 construct_point(p3), construct_point(p4)); 589 } 590 template<typename P> // can be Point or Point_3 591 Orientation orientation(const P& p1, const P& p2, const P& p3, const P& p4, 592 const Offset& o1, const Offset& o2, 593 const Offset& o3, const Offset& o4) const { 594 return geom_traits().orientation_3_object()( 595 construct_point(p1), construct_point(p2), construct_point(p3), construct_point(p4), 596 o1, o2, o3, o4); 597 } 598 template<typename P> // can be Point or Point_3 599 Orientation orientation(const std::pair<P, Offset>& pp1, const std::pair<P, Offset>& pp2, 600 const std::pair<P, Offset>& pp3, const std::pair<P, Offset>& pp4) const { 601 return geom_traits().orientation_3_object()( 602 construct_point(pp1.first), construct_point(pp2.first), 603 construct_point(pp3.first), construct_point(pp4.first), 604 pp1.second, pp2.second, pp3.second, pp4.second); 605 } 606 607 template<typename P> // can be Point or Point_3 608 bool equal(const P& p1, const P& p2) const { 609 return compare_xyz(construct_point(p1), construct_point(p2)) == EQUAL; 610 } 611 template<typename P> // can be Point or Point_3 612 bool equal(const P& p1, const P& p2, const Offset& o1, const Offset& o2) const { 613 return compare_xyz(construct_point(p1), construct_point(p2), o1, o2) == EQUAL; 614 } 615 template<typename P> // can be Point or Point_3 616 bool equal(const std::pair<P, Offset>& pp1, const std::pair<P, Offset>& pp2) const { 617 return compare_xyz(construct_point(pp1.first), construct_point(pp1.second), 618 pp2.first, pp2.second) == EQUAL; 619 } 620 621 template<typename P> // can be Point or Point_3 622 bool coplanar(const P& p1, const P& p2, const P& p3, const P& p4) const { 623 return orientation(construct_point(p1), construct_point(p2), 624 construct_point(p3), construct_point(p4)) == COPLANAR; 625 } 626 template<typename P> // can be Point or Point_3 627 bool coplanar(const P& p1, const P& p2, const P& p3, const P& p4, 628 const Offset& o1, const Offset& o2, const Offset& o3, const Offset& o4) const { 629 return orientation(construct_point(p1), construct_point(p2), 630 construct_point(p3), construct_point(p4), 631 o1, o2, o3, o4) == COPLANAR; 632 } 633 template<typename P> // can be Point or Point_3 634 bool coplanar(const std::pair<P, Offset>& pp1, const std::pair<P, Offset>& pp2, 635 const std::pair<P, Offset>& pp3, const std::pair<P, Offset>& pp4) const { 636 return orientation(construct_point(pp1.first), construct_point(pp2.first), 637 construct_point(pp3.first), construct_point(pp4.first), 638 pp1.second, pp2.second, pp3.second, pp4.second) == COPLANAR; 639 } 640 641 public: 642 /** @name Geometric access functions */ 643 644 // ************************************************************************* 645 // -*-*-*-*-*-*-*-*-*-*-*-*-*-* POINT -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 646 // ************************************************************************* 647 648 // Note that construct_point()-like functions have return type Point_3, 649 // but point() -like functions have return type Point 650 651 template<typename P> // can be Point or Point_3 652 decltype(auto) 653 construct_point(const P& p) const { 654 return geom_traits().construct_point_3_object()(p); 655 } 656 657 template<typename P> // can be Point or Point_3 658 Point_3 construct_point(const P& p, const Offset& o) const { 659 return geom_traits().construct_point_3_object()(p, o); 660 } 661 template<typename P> // can be Point or Point_3 662 Point_3 construct_point(const std::pair<P, Offset>& pp) const { 663 return construct_point(pp.first, pp.second); 664 } 665 666 Periodic_point_3 construct_periodic_point(const Point_3& p, 667 bool& had_to_use_exact) const 668 { 669 // The function is a different file to be able to be used where there is 670 // no triangulation (namely, the domains of Periodic_3_mesh_3). 671 return ::CGAL::P3T3::internal::construct_periodic_point(p, had_to_use_exact, geom_traits()); 672 } 673 674 Periodic_point_3 construct_periodic_point(const Point_3& p) const 675 { 676 bool useless = false; 677 return construct_periodic_point(p, useless); 678 } 679 680 // --------------------------------------------------------------------------- 681 // The following functions return objects of type Point and Periodic_point, 682 // _not_ Point_3 and Periodic_point_3. 683 // They are templated by `construct_point` to distingush between Delaunay and 684 // regular triangulations 685 // --------------------------------------------------------------------------- 686 687 template <class ConstructPoint> 688 Point point(const Periodic_point& pp, ConstructPoint cp) const 689 { 690 return cp(pp.first, pp.second); 691 } 692 693 // The following functions return the "real" position in space (unconstrained 694 // to the fundamental domain) of the vertices v and c->vertex(idx), 695 // respectively 696 template <class ConstructPoint> 697 Point point(Vertex_handle v, ConstructPoint cp) const { 698 return point(periodic_point(v), cp); 699 } 700 701 template <class ConstructPoint> 702 Point point(Cell_handle c, int idx, ConstructPoint cp) const 703 { 704 // if(is_1_cover()) 705 return point(periodic_point(c, idx), cp); 706 707 Offset vec_off[4]; 708 for(int i=0; i<4; i++) 709 vec_off[i] = periodic_point(c,i).second; 710 711 int ox = vec_off[0].x(); 712 int oy = vec_off[0].y(); 713 int oz = vec_off[0].z(); 714 for(int i=1; i<4; i++) { 715 ox = (std::min)(ox, vec_off[i].x()); 716 oy = (std::min)(oy, vec_off[i].y()); 717 oz = (std::min)(oz, vec_off[i].z()); 718 } 719 Offset diff_off(-ox, -oy, -oz); 720 if(diff_off.is_null()) 721 return point(periodic_point(c, idx), cp); 722 723 for(unsigned int i=0; i<4; i++) 724 vec_off[i] += diff_off; 725 Vertex_handle canonic_vh[4]; 726 for(int i=0; i<4; i++) { 727 Virtual_vertex_map_it vvmit = virtual_vertices.find(c->vertex(i)); 728 Vertex_handle orig_vh; 729 if(vvmit == virtual_vertices.end()) 730 orig_vh = c->vertex(i); 731 else 732 orig_vh = vvmit->second.first; 733 734 if(vec_off[i].is_null()) { 735 canonic_vh[i] = orig_vh; 736 } 737 else { 738 CGAL_assertion(virtual_vertices_reverse.find(orig_vh) 739 != virtual_vertices_reverse.end()); 740 canonic_vh[i] = virtual_vertices_reverse.find(orig_vh)-> 741 second[9*vec_off[i][0]+3*vec_off[i][1]+vec_off[i][2]-1]; 742 } 743 } 744 745 std::vector<Cell_handle> cells; 746 incident_cells(canonic_vh[0], std::back_inserter(cells)); 747 for(unsigned int i=0; i<cells.size(); i++) { 748 CGAL_assertion(cells[i]->has_vertex(canonic_vh[0])); 749 if(cells[i]->has_vertex(canonic_vh[1]) 750 && cells[i]->has_vertex(canonic_vh[2]) 751 && cells[i]->has_vertex(canonic_vh[3]) ) 752 return point(periodic_point(cells[i], cells[i]->index(canonic_vh[idx])), cp); 753 } 754 CGAL_assertion(false); 755 return Point(); 756 } 757 758 Periodic_point periodic_point(const Vertex_handle v) const 759 { 760 if(is_1_cover()) 761 return std::make_pair(v->point(), Offset(0,0,0)); 762 763 Virtual_vertex_map_it it = virtual_vertices.find(v); 764 if(it == virtual_vertices.end()) { 765 // if v is not contained in virtual_vertices, then it is in the 766 // original domain. 767 return std::make_pair(v->point(), Offset(0,0,0)); 768 } else { 769 // otherwise it has to be looked up as well as its offset. 770 return std::make_pair(it->second.first->point(), it->second.second); 771 } 772 } 773 774 Periodic_point periodic_point(const Cell_handle c, int i) const 775 { 776 if(is_1_cover()) 777 return std::make_pair(c->vertex(i)->point(), int_to_off(c->offset(i))); 778 779 Virtual_vertex_map_it it = virtual_vertices.find(c->vertex(i)); 780 if(it == virtual_vertices.end()) { 781 // if c->vertex(i) is not contained in virtual_vertices, then it 782 // is in the original domain. 783 return std::make_pair(c->vertex(i)->point(), 784 combine_offsets(Offset(),int_to_off(c->offset(i))) ); 785 } else { 786 // otherwise it has to be looked up as well as its offset. 787 return std::make_pair(it->second.first->point(), 788 combine_offsets(it->second.second, int_to_off(c->offset(i))) ); 789 } 790 } 791 792 // ************************************************************************* 793 // -*-*-*-*-*-*-*-*-*-*-*-*-*-* SEGMENT -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 794 // ************************************************************************* 795 796 // Note that an object of type `Segment` is constructed from two objects of type `Point_3` 797 // while an object of type `Periodic_segment` is constructed from two objects of type `Point` 798 799 template<typename P> // can be Point or Point_3 800 Segment construct_segment(const P& p1, const P& p2) const { 801 return geom_traits().construct_segment_3_object()(construct_point(p1), 802 construct_point(p2)); 803 } 804 template<typename P> // can be Point or Point_3 805 Segment construct_segment(const P& p1, const P& p2, 806 const Offset& o1, const Offset& o2) const { 807 return geom_traits().construct_segment_3_object()(construct_point(p1, o1), 808 construct_point(p2, o2)); 809 } 810 template<typename PS> // can be Periodic_segment or Periodic_segment_3 811 Segment construct_segment(const PS& ps) const { 812 return geom_traits().construct_segment_3_object()( 813 construct_point(ps[0].first, ps[0].second), 814 construct_point(ps[1].first, ps[1].second)); 815 } 816 template<typename P> // can be Point or Point_3 817 Periodic_segment construct_periodic_segment(const P& p1, const P& p2, 818 const Offset& o1 = Offset(), 819 const Offset& o2 = Offset()) const { 820 return CGAL::make_array(std::make_pair(construct_point(p1), o1), 821 std::make_pair(construct_point(p2), o2)); 822 } 823 824 Periodic_segment periodic_segment(const Cell_handle c, int i, int j) const { 825 CGAL_triangulation_precondition( i != j ); 826 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 827 CGAL_triangulation_precondition( i >= 0 && i <= 3 && j >= 0 && j <= 3 ); 828 return CGAL::make_array(std::make_pair(c->vertex(i)->point(), get_offset(c,i)), 829 std::make_pair(c->vertex(j)->point(), get_offset(c,j))); 830 } 831 Periodic_segment periodic_segment(const Edge& e) const { 832 return periodic_segment(e.first,e.second,e.third); 833 } 834 Periodic_segment periodic_segment(const Cell_handle c, Offset offset, int i, int j) const { 835 Periodic_segment result = periodic_segment(c,i,j); 836 offset.x() *= _cover[0]; 837 offset.y() *= _cover[1]; 838 offset.z() *= _cover[2]; 839 result[0].second += offset; 840 result[1].second += offset; 841 return result; 842 } 843 844 // ************************************************************************* 845 // -*-*-*-*-*-*-*-*-*-*-*-*-*-* TRIANGLE -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 846 // ************************************************************************* 847 848 // Note that an object of type `Triangle` is constructed from three objects of type `Point_3` 849 // while an object of type `Periodic_triangle` is constructed from three objects of type `Point` 850 851 template<typename P> // can be Point or Point_3 852 Triangle construct_triangle(const P& p1, const P& p2, const P& p3) const { 853 return geom_traits().construct_triangle_3_object()(construct_point(p1), 854 construct_point(p2), 855 construct_point(p3)); 856 } 857 template<typename P> // can be Point or Point_3 858 Triangle construct_triangle(const P& p1, const P& p2, const P& p3, 859 const Offset& o1, const Offset& o2, const Offset& o3) const { 860 return geom_traits().construct_triangle_3_object()(construct_point(p1, o1), 861 construct_point(p2, o2), 862 construct_point(p3, o3)); 863 } 864 template<typename PT> // can be Periodic_triangle or Periodic_triangle_3 865 Triangle construct_triangle(const PT& tri) const { 866 return geom_traits().construct_triangle_3_object()( 867 construct_point(tri[0].first, tri[0].second), 868 construct_point(tri[1].first, tri[1].second), 869 construct_point(tri[2].first, tri[2].second)); 870 } 871 template<typename P> // can be Point or Point_3 872 Periodic_triangle construct_periodic_triangle(const P& p1, const P& p2, 873 const P& p3, 874 const Offset& o1 = Offset(), 875 const Offset& o2 = Offset(), 876 const Offset& o3 = Offset()) const { 877 return CGAL::make_array(std::make_pair(construct_point(p1), o1), 878 std::make_pair(construct_point(p2), o2), 879 std::make_pair(construct_point(p3), o3)); 880 } 881 882 Periodic_triangle periodic_triangle(const Cell_handle c, int i) const { 883 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 884 CGAL_triangulation_precondition( i >= 0 && i <= 3 ); 885 if( (i&1)==0 ) 886 return CGAL::make_array(std::make_pair(c->vertex((i+2)&3)->point(), get_offset(c,(i+2)&3)), 887 std::make_pair(c->vertex((i+1)&3)->point(), get_offset(c,(i+1)&3)), 888 std::make_pair(c->vertex((i+3)&3)->point(), get_offset(c,(i+3)&3)) ); 889 890 return CGAL::make_array(std::make_pair(c->vertex((i+1)&3)->point(), get_offset(c,(i+1)&3)), 891 std::make_pair(c->vertex((i+2)&3)->point(), get_offset(c,(i+2)&3)), 892 std::make_pair(c->vertex((i+3)&3)->point(), get_offset(c,(i+3)&3)) ); 893 } 894 Periodic_triangle periodic_triangle(const Facet& f) const { 895 return periodic_triangle(f.first, f.second); 896 } 897 Periodic_triangle periodic_triangle(const Cell_handle c, Offset offset, int i) const { 898 Periodic_triangle result = periodic_triangle(c,i); 899 offset.x() *= _cover[0]; 900 offset.y() *= _cover[1]; 901 offset.z() *= _cover[2]; 902 result[0].second += offset; 903 result[1].second += offset; 904 result[2].second += offset; 905 return result; 906 } 907 908 // ************************************************************************* 909 // -*-*-*-*-*-*-*-*-*-*-*-*-*-* TETRAHEDRON -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 910 // ************************************************************************* 911 912 // Note that an object of type `Tetrahedron` is constructed from four objects of type `Point_3` 913 // while an object of type `Periodic_tetrahedron` is constructed from four objects of type `Point` 914 915 template<typename P> // can be Point or Point_3 916 Tetrahedron construct_tetrahedron(const P& p1, const P& p2, 917 const P& p3, const P& p4) const { 918 return geom_traits().construct_tetrahedron_3_object()(construct_point(p1), 919 construct_point(p2), 920 construct_point(p3), 921 construct_point(p4)); 922 } 923 template<typename P> // can be Point or Point_3 924 Tetrahedron construct_tetrahedron(const P& p1, const P& p2, 925 const P& p3, const P& p4, 926 const Offset& o1, const Offset& o2, 927 const Offset& o3, const Offset& o4) const { 928 return geom_traits().construct_tetrahedron_3_object()( 929 construct_point(p1, o1), construct_point(p2, o2), 930 construct_point(p3, o3), construct_point(p4, o4)); 931 } 932 template<typename PT> // can be Periodic_tetrahedron or Periodic_tetrahedron_3 933 Tetrahedron construct_tetrahedron(const PT& tet) const { 934 return geom_traits().construct_tetrahedron_3_object()( 935 construct_point(tet[0].first, tet[0].second), 936 construct_point(tet[1].first, tet[1].second), 937 construct_point(tet[2].first, tet[2].second), 938 construct_point(tet[3].first, tet[3].second)); 939 } 940 template<typename P> // can be Point or Point_3 941 Periodic_tetrahedron construct_periodic_tetrahedron(const P& p1, const P& p2, 942 const P& p3, const P& p4, 943 const Offset& o1 = Offset(), 944 const Offset& o2 = Offset(), 945 const Offset& o3 = Offset(), 946 const Offset& o4 = Offset()) const { 947 return CGAL::make_array(std::make_pair(construct_point(p1), o1), 948 std::make_pair(construct_point(p2), o2), 949 std::make_pair(construct_point(p3), o3), 950 std::make_pair(construct_point(p4), o4)); 951 } 952 953 Periodic_tetrahedron periodic_tetrahedron(const Cell_handle c) const 954 { 955 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 956 return CGAL::make_array(std::make_pair(c->vertex(0)->point(), get_offset(c,0)), 957 std::make_pair(c->vertex(1)->point(), get_offset(c,1)), 958 std::make_pair(c->vertex(2)->point(), get_offset(c,2)), 959 std::make_pair(c->vertex(3)->point(), get_offset(c,3)) ); 960 } 961 Periodic_tetrahedron periodic_tetrahedron(const Cell_handle c, Offset offset) const 962 { 963 Periodic_tetrahedron result = periodic_tetrahedron(c); 964 offset.x() *= _cover[0]; 965 offset.y() *= _cover[1]; 966 offset.z() *= _cover[2]; 967 result[0].second += offset; 968 result[1].second += offset; 969 result[2].second += offset; 970 result[3].second += offset; 971 return result; 972 } 973 974 // end of geometric functions 975 976 public: 977 /** @name Queries */ 978 bool is_vertex(const Point& p, Vertex_handle& v) const; 979 980 bool is_vertex(Vertex_handle v) const { 981 return _tds.is_vertex(v); 982 } 983 984 bool is_edge(Vertex_handle u, Vertex_handle v, Cell_handle& c, int& i, int& j) const { 985 return _tds.is_edge(u, v, c, i, j); 986 } 987 988 bool is_edge(Vertex_handle u, const Offset& off_u, 989 Vertex_handle v, const Offset& off_v, 990 Cell_handle& c, int& i, int& j) const 991 { 992 if(!_tds.is_edge(u,v,c,i,j)) 993 return false; 994 995 if((get_offset(c,i) == off_u) && (get_offset(c,j) == off_v)) 996 return true; 997 998 // it might be that different cells containing (u,v) yield 999 // different offsets, which forces us to test for all possibilities. 1000 else { 1001 Cell_circulator ccirc = incident_cells(c,i,j,c); 1002 while(++ccirc != c) { 1003 i = ccirc->index(u); 1004 j = ccirc->index(v); 1005 if((get_offset(ccirc,i) == off_u) && (get_offset(ccirc,j) == off_v)) { 1006 c = ccirc; 1007 return true; 1008 } 1009 } 1010 return false; 1011 } 1012 } 1013 1014 bool is_facet(Vertex_handle u, Vertex_handle v, Vertex_handle w, 1015 Cell_handle& c, int& i, int& j, int& k) const { 1016 return _tds.is_facet(u, v, w, c, i, j, k); 1017 } 1018 1019 bool is_facet(Vertex_handle u, const Offset& off_u, 1020 Vertex_handle v, const Offset& off_v, 1021 Vertex_handle w, const Offset& off_w, 1022 Cell_handle& c, int& i, int& j, int& k) const 1023 { 1024 if(!_tds.is_facet(u,v,w,c,i,j,k)) return false; 1025 if((get_offset(c,i) == off_u) 1026 && (get_offset(c,j) == off_v) 1027 && (get_offset(c,k) == off_w) ) 1028 return true; 1029 // it might be that c and c->neighbor(l) yield different offsets 1030 // which forces us to test for both possibilities. 1031 else { 1032 int l = 6-i-j-k; 1033 c = c->neighbor(l); 1034 i = c->index(u); 1035 j = c->index(v); 1036 k = c->index(w); 1037 return ((get_offset(c,i) == off_u) 1038 && (get_offset(c,j) == off_v) 1039 && (get_offset(c,k) == off_w) ); 1040 } 1041 } 1042 1043 bool is_cell(Cell_handle c) const { 1044 return _tds.is_cell(c); 1045 } 1046 1047 bool is_cell(Vertex_handle u, Vertex_handle v, 1048 Vertex_handle w, Vertex_handle t, 1049 Cell_handle& c, int& i, int& j, int& k, int& l) const { 1050 return _tds.is_cell(u, v, w, t, c, i, j, k, l); 1051 } 1052 1053 bool is_cell(Vertex_handle u, Vertex_handle v, Vertex_handle w, 1054 Vertex_handle t, Cell_handle& c) const 1055 { 1056 int i,j,k,l; 1057 return _tds.is_cell(u, v, w, t, c, i, j, k, l); 1058 } 1059 1060 bool is_cell(Vertex_handle u, const Offset& off_u, 1061 Vertex_handle v, const Offset& off_v, 1062 Vertex_handle w, const Offset& off_w, 1063 Vertex_handle t, const Offset& off_t, 1064 Cell_handle& c, int& i, int& j, int& k, int& l) const 1065 { 1066 if(!_tds.is_cell(u,v,w,t,c,i,j,k,l)) return false; 1067 return ((get_offset(c,i) == off_u) 1068 && (get_offset(c,j) == off_v) 1069 && (get_offset(c,k) == off_w) 1070 && (get_offset(c,l) == off_t) ); 1071 return false; 1072 } 1073 1074 bool is_cell(Vertex_handle u, const Offset& off_u, 1075 Vertex_handle v, const Offset& off_v, 1076 Vertex_handle w, const Offset& off_w, 1077 Vertex_handle t, const Offset& off_t, 1078 Cell_handle& c) const 1079 { 1080 int i, j, k, l; 1081 return is_cell(u,off_u,v,off_v,w,off_w,t,off_t,c,i,j,k,l); 1082 } 1083 1084 bool has_vertex(const Facet& f, Vertex_handle v, int& j) const { 1085 return _tds.has_vertex(f.first, f.second, v, j); 1086 } 1087 bool has_vertex(Cell_handle c, int i, Vertex_handle v, int& j) const { 1088 return _tds.has_vertex(c, i, v, j); 1089 } 1090 bool has_vertex(const Facet& f, Vertex_handle v) const { 1091 return _tds.has_vertex(f.first, f.second, v); 1092 } 1093 bool has_vertex(Cell_handle c, int i, Vertex_handle v) const { 1094 return _tds.has_vertex(c, i, v); 1095 } 1096 1097 bool are_equal(Cell_handle c, int i, Cell_handle n, int j) const { 1098 return _tds.are_equal(c, i, n, j); 1099 } 1100 bool are_equal(const Facet& f, const Facet& g) const { 1101 return _tds.are_equal(f.first, f.second, g.first, g.second); 1102 } 1103 bool are_equal(const Facet& f, Cell_handle n, int j) const { 1104 return _tds.are_equal(f.first, f.second, n, j); 1105 } 1106 1107 #ifdef CGAL_NO_STRUCTURAL_FILTERING 1108 Cell_handle 1109 periodic_locate(const Point& p, const Offset& o_p, Offset& lo, 1110 Locate_type& lt, int& li, int& lj, 1111 Cell_handle start = Cell_handle()) const; 1112 #else // no CGAL_NO_STRUCTURAL_FILTERING 1113 # ifndef CGAL_P3T3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS 1114 # define CGAL_P3T3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS 2500 1115 # endif // no CGAL_P3T3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS 1116 1117 public: 1118 Cell_handle 1119 inexact_periodic_locate(const Point& p, const Offset& o_p, 1120 Cell_handle start = Cell_handle(), 1121 int max_num_cells = CGAL_P3T3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const; 1122 protected: 1123 Cell_handle 1124 exact_periodic_locate(const Point& p, const Offset& o_p, Offset& lo, 1125 Locate_type& lt, 1126 int& li, int& lj, 1127 Cell_handle start) const; 1128 1129 Cell_handle 1130 generic_periodic_locate(const Point& p, const Offset& o_p, Offset& lo, 1131 Locate_type& lt, 1132 int& li, int& lj, 1133 Cell_handle start, 1134 internal::Periodic_structural_filtering_3_tag) const { 1135 return exact_periodic_locate(p, o_p, lo, lt, li, lj, inexact_periodic_locate(p, o_p, start)); 1136 } 1137 1138 Cell_handle 1139 generic_periodic_locate(const Point& p, const Offset& o_p, Offset& lo, 1140 Locate_type& lt, 1141 int& li, int& lj, 1142 Cell_handle start, 1143 internal::No_periodic_structural_filtering_3_tag) const { 1144 return exact_periodic_locate(p, o_p, lo, lt, li, lj, start); 1145 } 1146 1147 Orientation 1148 inexact_orientation(const Point& p, const Point& q, 1149 const Point& r, const Point& s, 1150 const Offset& o_p = Offset(), const Offset& o_q = Offset(), 1151 const Offset& o_r = Offset(), const Offset& o_s = Offset()) const 1152 { 1153 const Point_3 bp = construct_point(p, o_p); 1154 const Point_3 bq = construct_point(q, o_q); 1155 const Point_3 br = construct_point(r, o_r); 1156 const Point_3 bs = construct_point(s, o_s); 1157 1158 // So that this code works well with Lazy_kernel 1159 internal::Static_filters_predicates::Get_approx<Point_3> get_approx; 1160 const double px = to_double(get_approx(bp).x()); 1161 const double py = to_double(get_approx(bp).y()); 1162 const double pz = to_double(get_approx(bp).z()); 1163 const double qx = to_double(get_approx(bq).x()); 1164 const double qy = to_double(get_approx(bq).y()); 1165 const double qz = to_double(get_approx(bq).z()); 1166 const double rx = to_double(get_approx(br).x()); 1167 const double ry = to_double(get_approx(br).y()); 1168 const double rz = to_double(get_approx(br).z()); 1169 const double sx = to_double(get_approx(bs).x()); 1170 const double sy = to_double(get_approx(bs).y()); 1171 const double sz = to_double(get_approx(bs).z()); 1172 1173 const double pqx = qx - px; 1174 const double pqy = qy - py; 1175 const double pqz = qz - pz; 1176 const double prx = rx - px; 1177 const double pry = ry - py; 1178 const double prz = rz - pz; 1179 const double psx = sx - px; 1180 const double psy = sy - py; 1181 const double psz = sz - pz; 1182 1183 const double det = determinant(pqx, pqy, pqz, 1184 prx, pry, prz, 1185 psx, psy, psz); 1186 if(det > 0) return POSITIVE; 1187 if(det < 0) return NEGATIVE; 1188 return ZERO; 1189 } 1190 1191 public: 1192 Cell_handle 1193 periodic_locate(const Point& p, const Offset& o_p, Offset& lo, 1194 Locate_type& lt, int& li, int& lj, 1195 Cell_handle start = Cell_handle()) const 1196 { 1197 typedef Triangulation_structural_filtering_traits<Geometric_traits> TSFT; 1198 typedef typename internal::Periodic_structural_filtering_selector_3< 1199 TSFT::Use_structural_filtering_tag::value >::Tag Should_filter_tag; 1200 1201 return generic_periodic_locate(p, o_p, lo, lt, li, lj, start, Should_filter_tag()); 1202 } 1203 1204 Cell_handle 1205 inexact_locate(const Point& p, 1206 Cell_handle start = Cell_handle(), 1207 int max_num_cells = CGAL_P3T3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const 1208 { 1209 return inexact_periodic_locate(p, Offset(), start, max_num_cells); 1210 } 1211 #endif // no CGAL_NO_STRUCTURAL_FILTERING 1212 1213 protected: 1214 /** @name Location helpers */ 1215 // Cell_handle periodic_locate(const Point& p, const Offset& o_p, 1216 // Locate_type& lt, int& li, int& lj, Cell_handle start) const; 1217 1218 Bounded_side side_of_cell(const Point& p, const Offset& off, Cell_handle c, 1219 Locate_type& lt, int& i, int& j) const; 1220 1221 public: 1222 /** @name Point Location */ 1223 /** Wrapper function for locate if only the request point is given. 1224 */ 1225 Cell_handle locate(const Point& p, Cell_handle start = Cell_handle()) const 1226 { 1227 Locate_type lt; 1228 int li, lj; 1229 return locate( p, lt, li, lj, start); 1230 } 1231 1232 /** Wrapper function calling locate with an empty offset if there was no 1233 * offset given. 1234 */ 1235 Cell_handle locate(const Point& p, 1236 Locate_type& lt, int& li, int& lj, 1237 Cell_handle start = Cell_handle()) const 1238 { 1239 Offset lo; 1240 return locate( p, lo, lt, li, lj, start); 1241 } 1242 1243 Cell_handle locate(const Point& p, Offset& lo, 1244 Cell_handle start = Cell_handle()) const 1245 { 1246 Locate_type lt; 1247 int li, lj; 1248 return locate( p, lo, lt, li, lj, start); 1249 } 1250 1251 Cell_handle locate(const Point& p, Offset& lo, 1252 Locate_type& lt, int& li, int& lj, 1253 Cell_handle start = Cell_handle()) const 1254 { 1255 Cell_handle ch = periodic_locate(p, Offset(), lo, lt, li, lj, start); 1256 for(unsigned i = 0; i < 3; ++i) 1257 if(lo[i] >= 1) 1258 lo[i] = -1; 1259 return ch; 1260 } 1261 1262 Bounded_side side_of_cell(const Point& p, Cell_handle c, 1263 Locate_type& lt, int& i, int& j) const 1264 { 1265 if(number_of_vertices() == 0) { 1266 lt = EMPTY; 1267 return ON_UNBOUNDED_SIDE; 1268 } 1269 return side_of_cell(p,Offset(),c,lt,i,j); 1270 } 1271 1272 private: 1273 /** @name Insertion helpers */ 1274 template < class Conflict_tester, class Point_hider, class CoverManager > 1275 Vertex_handle periodic_insert(const Point& p, const Offset& o, Locate_type lt, 1276 Cell_handle c, const Conflict_tester& tester, 1277 Point_hider& hider, CoverManager& cover_manager, 1278 Vertex_handle vh = Vertex_handle()); 1279 1280 template <class Point_iterator, class Offset_iterator> 1281 void periodic_sort(Point_iterator /*p_begin*/, Point_iterator /*p_end*/, 1282 Offset_iterator /*o_begin*/, Offset_iterator /*o_end*/) const { 1283 std::cout << "Periodic_sort not yet implemented" << std::endl; 1284 } 1285 1286 Vertex_handle create_initial_triangulation(const Point& p); 1287 1288 public: 1289 std::vector<Vertex_handle> insert_dummy_points(); 1290 1291 protected: 1292 // this is needed for compatibility reasons 1293 template <class Conflict_test, class OutputIteratorBoundaryFacets, 1294 class OutputIteratorCells, class OutputIteratorInternalFacets> 1295 Triple<OutputIteratorBoundaryFacets, OutputIteratorCells, 1296 OutputIteratorInternalFacets> 1297 find_conflicts(Cell_handle c, 1298 const Conflict_test& tester, 1299 Triple<OutputIteratorBoundaryFacets, OutputIteratorCells, 1300 OutputIteratorInternalFacets> it) const 1301 { 1302 bool b = false; 1303 Offset off = get_location_offset(tester, c, b); 1304 if(b) 1305 return find_conflicts(c,off,tester,it); 1306 return it; 1307 } 1308 1309 template <class Conflict_test, class OutputIteratorBoundaryFacets, 1310 class OutputIteratorCells, class OutputIteratorInternalFacets> 1311 Triple<OutputIteratorBoundaryFacets, OutputIteratorCells, 1312 OutputIteratorInternalFacets> 1313 find_conflicts(Cell_handle c, const Offset& current_off, 1314 const Conflict_test& tester, 1315 Triple<OutputIteratorBoundaryFacets, 1316 OutputIteratorCells, 1317 OutputIteratorInternalFacets> it) const; 1318 1319 protected: 1320 1321 // COMMON INSERTION for DELAUNAY and REGULAR TRIANGULATION 1322 template < class Conflict_tester, class Point_hider, class CoverManager > 1323 Vertex_handle insert_in_conflict(const Point& p, Cell_handle start, 1324 const Conflict_tester& tester, 1325 Point_hider& hider, 1326 CoverManager& cover_manager) 1327 { 1328 Locate_type lt = Locate_type(); 1329 int li=0, lj=0; 1330 Offset lo; 1331 Cell_handle c = periodic_locate(p, Offset(), lo, lt, li, lj, start); 1332 return insert_in_conflict(p,lt,c,li,lj,tester,hider, cover_manager); 1333 } 1334 1335 template < class Conflict_tester, class Point_hider, class CoverManager > 1336 Vertex_handle insert_in_conflict(const Point& p, Locate_type lt, 1337 Cell_handle c, int li, int lj, 1338 const Conflict_tester& tester, 1339 Point_hider& hider, 1340 CoverManager& cover_manager); 1341 1342 template < class InputIterator, class Conflict_tester, 1343 class Point_hider, class CoverManager> 1344 std::vector<Vertex_handle> insert_in_conflict( 1345 InputIterator begin, InputIterator end, Cell_handle start, 1346 Conflict_tester& tester, Point_hider& hider, CoverManager& cover_manager) 1347 { 1348 Vertex_handle new_vertex; 1349 std::vector<Vertex_handle> double_vertices; 1350 Locate_type lt = Locate_type(); 1351 int li=0, lj=0; 1352 CGAL_triangulation_assertion_code( Locate_type lta = Locate_type(); ) 1353 CGAL_triangulation_assertion_code( int ia = 0; ) 1354 CGAL_triangulation_assertion_code( int ja = 0; ) 1355 Cell_handle hint; 1356 while(begin!=end) { 1357 tester.set_point(*begin); 1358 Offset lo; 1359 hint = periodic_locate(*begin, Offset(), lo, lt, li, lj, start); 1360 CGAL_triangulation_assertion_code( if(number_of_vertices() != 0) { ); 1361 CGAL_triangulation_assertion(side_of_cell( 1362 *begin,Offset(), hint, lta, ia, ja) != ON_UNBOUNDED_SIDE); 1363 CGAL_triangulation_assertion(lta == lt); 1364 CGAL_triangulation_assertion(ia == li); 1365 CGAL_triangulation_assertion(ja == lj); 1366 CGAL_triangulation_assertion_code( } 1367 ); 1368 1369 new_vertex = insert_in_conflict(*begin,lt,hint,li,lj,tester,hider, cover_manager); 1370 if(lt == VERTEX) 1371 double_vertices.push_back(new_vertex); 1372 if(new_vertex != Vertex_handle()) 1373 start = new_vertex->cell(); 1374 begin++; 1375 } 1376 return double_vertices; 1377 } 1378 1379 private: 1380 /** @name Removal helpers */ 1381 Vertex_triple make_vertex_triple(const Facet& f) const 1382 { 1383 Cell_handle ch = f.first; 1384 int i = f.second; 1385 return Vertex_triple(ch->vertex(vertex_triple_index(i,0)), 1386 ch->vertex(vertex_triple_index(i,1)), 1387 ch->vertex(vertex_triple_index(i,2))); 1388 } 1389 1390 void make_canonical(Vertex_triple& t) const; 1391 1392 void make_hole(Vertex_handle v, std::map<Vertex_triple,Facet>& outer_map, 1393 std::vector<Cell_handle>& hole); 1394 1395 protected: 1396 /** @name Removal */ 1397 template < class PointRemover, class CoverManager > 1398 bool periodic_remove(Vertex_handle v, PointRemover& remover, CoverManager& cover_manager, 1399 const bool abort_if_cover_change = false); 1400 1401 template < class PointRemover, class CT, class CoverManager > 1402 void remove(Vertex_handle v, PointRemover& remover, CT& ct, CoverManager& cover_manager); 1403 1404 void delete_vertex(Vertex_handle vertex_handle) 1405 { 1406 tds().delete_vertex(vertex_handle); 1407 1408 if(!is_1_cover()) 1409 { 1410 typename Virtual_vertex_map::iterator iter = this->virtual_vertices.find(vertex_handle); 1411 if(iter != this->virtual_vertices.end()) 1412 { 1413 Vertex_handle vh = iter->second.first; 1414 this->virtual_vertices.erase(iter); 1415 1416 typename Virtual_vertex_reverse_map::iterator origin_it = this->virtual_vertices_reverse.find(vh); 1417 std::vector<Vertex_handle>& copies = origin_it->second; 1418 typename std::vector<Vertex_handle>::iterator copy_iter = std::find(copies.begin(), copies.end(), vertex_handle); 1419 CGAL_triangulation_assertion(copy_iter != copies.end()); 1420 copies.erase(copy_iter); 1421 if(copies.empty()) 1422 virtual_vertices_reverse.erase(origin_it); 1423 } 1424 return; 1425 } 1426 1427 CGAL_triangulation_assertion(this->virtual_vertices.find(vertex_handle) == this->virtual_vertices.end()); 1428 CGAL_triangulation_assertion(this->virtual_vertices_reverse.find(vertex_handle) == this->virtual_vertices_reverse.end()); 1429 } 1430 1431 public: 1432 /** @name Traversal */ 1433 Cell_iterator cells_begin() const { 1434 return _tds.cells_begin(); 1435 } 1436 Cell_iterator cells_end() const { 1437 return _tds.cells_end(); 1438 } 1439 1440 Vertex_iterator vertices_begin() const { 1441 return _tds.vertices_begin(); 1442 } 1443 Vertex_iterator vertices_end() const { 1444 return _tds.vertices_end(); 1445 } 1446 1447 Edge_iterator edges_begin() const { 1448 return _tds.edges_begin(); 1449 } 1450 Edge_iterator edges_end() const { 1451 return _tds.edges_end(); 1452 } 1453 1454 Facet_iterator facets_begin() const { 1455 return _tds.facets_begin(); 1456 } 1457 Facet_iterator facets_end() const { 1458 return _tds.facets_end(); 1459 } 1460 1461 Cell_iterator finite_cells_begin() const { 1462 return _tds.cells_begin(); 1463 } 1464 Cell_iterator finite_cells_end() const { 1465 return _tds.cells_end(); 1466 } 1467 1468 Vertex_iterator finite_vertices_begin() const { 1469 return _tds.vertices_begin(); 1470 } 1471 Vertex_iterator finite_vertices_end() const { 1472 return _tds.vertices_end(); 1473 } 1474 1475 Edge_iterator finite_edges_begin() const { 1476 return _tds.edges_begin(); 1477 } 1478 Edge_iterator finite_edges_end() const { 1479 return _tds.edges_end(); 1480 } 1481 1482 Facet_iterator finite_facets_begin() const { 1483 return _tds.facets_begin(); 1484 } 1485 Facet_iterator finite_facets_end() const { 1486 return _tds.facets_end(); 1487 } 1488 1489 All_cells_iterator all_cells_begin() const { 1490 return _tds.cells_begin(); 1491 } 1492 All_cells_iterator all_cells_end() const { 1493 return _tds.cells_end(); 1494 } 1495 1496 All_vertices_iterator all_vertices_begin() const { 1497 return _tds.vertices_begin(); 1498 } 1499 All_vertices_iterator all_vertices_end() const { 1500 return _tds.vertices_end(); 1501 } 1502 1503 All_edges_iterator all_edges_begin() const { 1504 return _tds.edges_begin(); 1505 } 1506 All_edges_iterator all_edges_end() const { 1507 return _tds.edges_end(); 1508 } 1509 1510 All_facets_iterator all_facets_begin() const { 1511 return _tds.facets_begin(); 1512 } 1513 All_facets_iterator all_facets_end() const { 1514 return _tds.facets_end(); 1515 } 1516 1517 Unique_vertex_iterator unique_vertices_begin() const { 1518 return CGAL::filter_iterator(vertices_end(), Domain_tester<Self>(this), 1519 vertices_begin()); 1520 } 1521 Unique_vertex_iterator unique_vertices_end() const { 1522 return CGAL::filter_iterator(vertices_end(), Domain_tester<Self>(this)); 1523 } 1524 1525 // Geometric iterators 1526 Periodic_tetrahedron_iterator periodic_tetrahedra_begin( 1527 Iterator_type it = STORED) const { 1528 return Periodic_tetrahedron_iterator(this, it); 1529 } 1530 Periodic_tetrahedron_iterator periodic_tetrahedra_end( 1531 Iterator_type it = STORED) const { 1532 return Periodic_tetrahedron_iterator(this, 1, it); 1533 } 1534 1535 Periodic_triangle_iterator periodic_triangles_begin( 1536 Iterator_type it = STORED) const { 1537 return Periodic_triangle_iterator(this, it); 1538 } 1539 Periodic_triangle_iterator periodic_triangles_end( 1540 Iterator_type it = STORED) const { 1541 return Periodic_triangle_iterator(this, 1, it); 1542 } 1543 1544 Periodic_segment_iterator periodic_segments_begin( 1545 Iterator_type it = STORED) const { 1546 return Periodic_segment_iterator(this, it); 1547 } 1548 Periodic_segment_iterator periodic_segments_end( 1549 Iterator_type it = STORED) const { 1550 return Periodic_segment_iterator(this, 1, it); 1551 } 1552 1553 Periodic_point_iterator periodic_points_begin( 1554 Iterator_type it = STORED) const { 1555 return Periodic_point_iterator(this, it); 1556 } 1557 Periodic_point_iterator periodic_points_end( 1558 Iterator_type it = STORED) const { 1559 return Periodic_point_iterator(this, 1, it); 1560 } 1561 1562 // Circulators 1563 Cell_circulator incident_cells(const Edge& e) const { 1564 return _tds.incident_cells(e); 1565 } 1566 Cell_circulator incident_cells(Cell_handle c, int i, int j) const { 1567 return _tds.incident_cells(c, i, j); 1568 } 1569 Cell_circulator incident_cells(const Edge& e, Cell_handle start) const { 1570 return _tds.incident_cells(e, start); 1571 } 1572 Cell_circulator incident_cells(Cell_handle c, int i, int j, 1573 Cell_handle start) const { 1574 return _tds.incident_cells(c, i, j, start); 1575 } 1576 1577 Facet_circulator incident_facets(const Edge& e) const { 1578 return _tds.incident_facets(e); 1579 } 1580 Facet_circulator incident_facets(Cell_handle c, int i, int j) const { 1581 return _tds.incident_facets(c, i, j); 1582 } 1583 Facet_circulator incident_facets(const Edge& e, const Facet& start) const { 1584 return _tds.incident_facets(e, start); 1585 } 1586 Facet_circulator incident_facets(Cell_handle c, int i, int j, 1587 const Facet& start) const { 1588 return _tds.incident_facets(c, i, j, start); 1589 } 1590 Facet_circulator incident_facets(const Edge& e, 1591 Cell_handle start, int f) const { 1592 return _tds.incident_facets(e, start, f); 1593 } 1594 Facet_circulator incident_facets(Cell_handle c, int i, int j, 1595 Cell_handle start, int f) const { 1596 return _tds.incident_facets(c, i, j, start, f); 1597 } 1598 1599 // around a vertex 1600 template <class OutputIterator> 1601 OutputIterator incident_cells(Vertex_handle v, OutputIterator cells) const { 1602 return _tds.incident_cells(v, cells); 1603 } 1604 1605 template <class OutputIterator> 1606 OutputIterator incident_facets(Vertex_handle v, OutputIterator facets) const { 1607 return _tds.incident_facets(v, facets); 1608 } 1609 1610 template <class OutputIterator> 1611 OutputIterator incident_edges(Vertex_handle v, OutputIterator edges) const { 1612 return _tds.incident_edges(v, edges); 1613 } 1614 1615 template <class OutputIterator> 1616 OutputIterator adjacent_vertices(Vertex_handle v, OutputIterator vertices) const { 1617 return _tds.adjacent_vertices(v, vertices); 1618 } 1619 1620 //deprecated, don't use anymore 1621 template <class OutputIterator> 1622 OutputIterator incident_vertices(Vertex_handle v, OutputIterator vertices) const { 1623 return _tds.adjacent_vertices(v, vertices); 1624 } 1625 1626 size_type degree(Vertex_handle v) const { 1627 return _tds.degree(v); 1628 } 1629 1630 // Functions forwarded from TDS. 1631 int mirror_index(Cell_handle c, int i) const { 1632 return _tds.mirror_index(c, i); 1633 } 1634 1635 Vertex_handle mirror_vertex(Cell_handle c, int i) const { 1636 return _tds.mirror_vertex(c, i); 1637 } 1638 1639 Facet mirror_facet(Facet f) const { 1640 return _tds.mirror_facet(f); 1641 } 1642 1643 private: 1644 /** @name Checking helpers */ 1645 /// calls has_self_edges for every cell of the triangulation 1646 bool has_self_edges() const 1647 { 1648 Cell_iterator it; 1649 for( it = all_cells_begin(); it != all_cells_end(); ++it ) 1650 if(has_self_edges(it)) return true; 1651 return false; 1652 } 1653 1654 bool has_self_edges(Cell_handle c) const; 1655 1656 public: 1657 /** @name Checking */ 1658 bool is_valid(bool verbose = false, int level = 0) const; 1659 bool is_valid(Cell_handle c, bool verbose = false, int level = 0) const; 1660 1661 protected: 1662 template <class ConflictTester> 1663 bool is_valid_conflict(ConflictTester& tester, bool verbose = false, 1664 int level = 0) const; 1665 1666 public: 1667 /** @name Functors */ 1668 class Perturbation_order 1669 { 1670 const Self *t; 1671 1672 public: 1673 Perturbation_order(const Self *tr) 1674 : t(tr) {} 1675 1676 template<typename P> 1677 bool operator()(const P* p, const P* q) const { 1678 return t->compare_xyz(*p, *q) == SMALLER; 1679 } 1680 }; 1681 1682 public: 1683 // undocumented access functions 1684 Offset get_offset(Cell_handle ch, int i) const 1685 { 1686 if(is_1_cover()) 1687 return int_to_off(ch->offset(i)); 1688 1689 Virtual_vertex_map_it it = virtual_vertices.find(ch->vertex(i)); 1690 if(it != virtual_vertices.end()) 1691 return combine_offsets(it->second.second, int_to_off(ch->offset(i))); 1692 else 1693 return combine_offsets(Offset(), int_to_off(ch->offset(i))); 1694 } 1695 1696 Offset get_offset(Vertex_handle vh) const 1697 { 1698 if(is_1_cover()) 1699 return Offset(); 1700 1701 Virtual_vertex_map_it it = virtual_vertices.find(vh); 1702 if(it != virtual_vertices.end()) 1703 return it->second.second; 1704 else 1705 return Offset(); 1706 } 1707 1708 Vertex_handle get_original_vertex(Vertex_handle vh) const 1709 { 1710 if(is_1_cover()) 1711 return vh; 1712 1713 Virtual_vertex_map_it it = virtual_vertices.find(vh); 1714 if(it != virtual_vertices.end()) 1715 return it->second.first; 1716 else 1717 return vh; 1718 } 1719 1720 Offset combine_offsets(const Offset& o_c, const Offset& o_t) const 1721 { 1722 Offset o_ct(_cover[0]*o_t.x(), _cover[1]*o_t.y(), _cover[2]*o_t.z()); 1723 return o_c + o_ct; 1724 } 1725 1726 // These functions give the pair (vertex, offset) that corresponds to the 1727 // i-th vertex of cell ch or vertex vh, respectively. 1728 void get_vertex(Cell_handle ch, int i, Vertex_handle& vh, Offset& off) const; 1729 void get_vertex(Vertex_handle vh_i, Vertex_handle& vh, Offset& off) const; 1730 1731 protected: 1732 // Auxiliary functions 1733 Cell_handle get_cell(const Vertex_handle* vh) const; 1734 1735 template<class Conflict_tester> 1736 Offset get_location_offset(const Conflict_tester& tester, 1737 Cell_handle c) const; 1738 1739 template<class Conflict_tester> 1740 Offset get_location_offset(const Conflict_tester& tester, 1741 Cell_handle c, bool& found) const; 1742 1743 Offset neighbor_offset(Cell_handle ch, int i, Cell_handle nb) const; 1744 1745 public: 1746 Offset neighbor_offset(Cell_handle ch, int i) const 1747 { 1748 return neighbor_offset(ch, i, ch->neighbor(i)); 1749 } 1750 1751 protected: 1752 /** @name Friends */ 1753 friend std::istream& operator>> <> 1754 (std::istream& is, Periodic_3_triangulation_3<GT,TDS>& tr); 1755 friend std::ostream& operator<< <> 1756 (std::ostream& os, const Periodic_3_triangulation_3<GT,TDS>& tr); 1757 1758 protected: 1759 template <class ConstructCircumcenter> 1760 Periodic_point_3 periodic_circumcenter(Cell_handle c, 1761 ConstructCircumcenter construct_circumcenter) const 1762 { 1763 CGAL_triangulation_precondition(c != Cell_handle()); 1764 1765 Point_3 p = construct_circumcenter(c->vertex(0)->point(), c->vertex(1)->point(), 1766 c->vertex(2)->point(), c->vertex(3)->point(), 1767 get_offset(c, 0), get_offset(c, 1), 1768 get_offset(c, 2), get_offset(c, 3)); 1769 1770 return construct_periodic_point(p); 1771 } 1772 1773 public: 1774 bool is_canonical(const Facet& f) const 1775 { 1776 if(number_of_sheets() == CGAL::make_array(1,1,1)) 1777 return true; 1778 1779 Offset cell_off0 = int_to_off(f.first->offset((f.second+1)&3)); 1780 Offset cell_off1 = int_to_off(f.first->offset((f.second+2)&3)); 1781 Offset cell_off2 = int_to_off(f.first->offset((f.second+3)&3)); 1782 Offset diff_off((cell_off0.x() == 1 1783 && cell_off1.x() == 1 1784 && cell_off2.x() == 1) ? -1 : 0, 1785 (cell_off0.y() == 1 1786 && cell_off1.y() == 1 1787 && cell_off2.y() == 1) ? -1 : 0, 1788 (cell_off0.z() == 1 1789 && cell_off1.z() == 1 1790 && cell_off2.z() == 1) ? -1 : 0); 1791 Offset off0 = combine_offsets(get_offset(f.first, (f.second+1)&3), diff_off); 1792 Offset off1 = combine_offsets(get_offset(f.first, (f.second+2)&3), diff_off); 1793 Offset off2 = combine_offsets(get_offset(f.first, (f.second+3)&3), diff_off); 1794 1795 // If there is one offset with entries larger than 1 then we are 1796 // talking about a vertex that is too far away from the original 1797 // domain to belong to a canonical triangle. 1798 if(off0.x() > 1) return false; 1799 if(off0.y() > 1) return false; 1800 if(off0.z() > 1) return false; 1801 if(off1.x() > 1) return false; 1802 if(off1.y() > 1) return false; 1803 if(off1.z() > 1) return false; 1804 if(off2.x() > 1) return false; 1805 if(off2.y() > 1) return false; 1806 if(off2.z() > 1) return false; 1807 1808 // If there is one direction of space for which all offsets are 1809 // non-zero then the edge is not canonical because we can 1810 // take the copy closer towards the origin in that direction. 1811 int offx = off0.x() & off1.x() & off2.x(); 1812 int offy = off0.y() & off1.y() & off2.y(); 1813 int offz = off0.z() & off1.z() & off2.z(); 1814 1815 return (offx == 0 && offy == 0 && offz == 0); 1816 } 1817 1818 protected: 1819 template <class ConstructCircumcenter> 1820 bool canonical_dual_segment(Cell_handle c, int i, Periodic_segment_3& ps, 1821 ConstructCircumcenter construct_circumcenter) const 1822 { 1823 CGAL_triangulation_precondition(c != Cell_handle()); 1824 Offset off = neighbor_offset(c,i,c->neighbor(i)); 1825 Periodic_point_3 p1 = periodic_circumcenter(c, construct_circumcenter); 1826 Periodic_point_3 p2 = periodic_circumcenter(c->neighbor(i), construct_circumcenter); 1827 Offset o1 = -p1.second; 1828 Offset o2 = combine_offsets(-p2.second,-off); 1829 Offset cumm_off((std::min)(o1.x(),o2.x()), 1830 (std::min)(o1.y(),o2.y()), 1831 (std::min)(o1.z(),o2.z())); 1832 const Periodic_point_3 pp1 = std::make_pair(construct_point(p1), o1-cumm_off); 1833 const Periodic_point_3 pp2 = std::make_pair(construct_point(p2), o2-cumm_off); 1834 ps = CGAL::make_array(pp1, pp2); 1835 return (cumm_off == Offset(0,0,0)); 1836 } 1837 1838 template <class OutputIterator, class ConstructCircumcenter> 1839 OutputIterator dual(Cell_handle c, int i, int j, 1840 OutputIterator points, 1841 ConstructCircumcenter construct_circumcenter) const 1842 { 1843 Cell_circulator cstart = incident_cells(c, i, j); 1844 1845 Offset offv = periodic_point(c,i).second; 1846 Vertex_handle v = c->vertex(i); 1847 1848 Cell_circulator ccit = cstart; 1849 do { 1850 Point_3 dual_orig = periodic_circumcenter(ccit, construct_circumcenter).first; 1851 int idx = ccit->index(v); 1852 Offset off = periodic_point(ccit,idx).second; 1853 Point_3 dual = construct_point(std::make_pair(dual_orig, -off+offv)); 1854 *points++ = dual; 1855 ++ccit; 1856 } while(ccit != cstart); 1857 1858 return points; 1859 } 1860 1861 template <class OutputIterator, class ConstructCircumcenter> 1862 OutputIterator dual(Vertex_handle v, OutputIterator points, 1863 ConstructCircumcenter construct_circumcenter) const 1864 { 1865 std::vector<Cell_handle> cells; 1866 incident_cells(v,std::back_inserter(cells)); 1867 1868 for(unsigned int i=0; i<cells.size(); i++) { 1869 Point_3 dual_orig = periodic_circumcenter(cells[i], construct_circumcenter).first; 1870 int idx = cells[i]->index(v); 1871 Offset off = periodic_point(cells[i],idx).second; 1872 Point_3 dual = construct_point(std::make_pair(dual_orig, -off)); 1873 *points++ = dual; 1874 } 1875 1876 return points; 1877 } 1878 1879 template <class Stream, class ConstructCircumcenter> 1880 Stream& draw_dual(Stream& os, ConstructCircumcenter construct_circumcenter) const 1881 { 1882 CGAL_triangulation_assertion_code( unsigned int i = 0; ) 1883 for(Facet_iterator fit = facets_begin(), end = facets_end(); 1884 fit != end; ++fit) { 1885 if(!is_canonical(*fit)) 1886 continue; 1887 1888 Periodic_segment_3 pso; 1889 canonical_dual_segment(fit->first, fit->second, pso, construct_circumcenter); 1890 Segment so = construct_segment(pso); 1891 CGAL_triangulation_assertion_code ( ++i; ) 1892 os << so.source() << '\n'; 1893 os << so.target() << '\n'; 1894 } 1895 CGAL_triangulation_assertion( i == number_of_facets() ); 1896 1897 return os; 1898 } 1899 1900 /// Volume computations 1901 1902 // Note: Polygon area computation requires to evaluate square roots 1903 // and thus cannot be done without changing the Traits concept. 1904 1905 template <class ConstructCircumcenter> 1906 FT dual_volume(Vertex_handle v, ConstructCircumcenter construct_circumcenter) const 1907 { 1908 std::list<Edge> edges; 1909 incident_edges(v, std::back_inserter(edges)); 1910 1911 FT vol(0); 1912 for(typename std::list<Edge>::iterator eit = edges.begin(); 1913 eit != edges.end(); ++eit) { 1914 1915 // compute the dual of the edge *eit but handle the translations 1916 // with respect to the dual of v. That is why we cannot use one 1917 // of the existing dual functions here. 1918 Facet_circulator fstart = incident_facets(*eit); 1919 Facet_circulator fcit = fstart; 1920 std::vector<Point_3> pts; 1921 do { 1922 // TODO: possible speed-up by caching the circumcenters 1923 Point_3 dual_orig = periodic_circumcenter(fcit->first, construct_circumcenter).first; 1924 int idx = fcit->first->index(v); 1925 Offset off = periodic_point(fcit->first,idx).second; 1926 pts.push_back(construct_point(std::make_pair(dual_orig,-off))); 1927 ++fcit; 1928 } while(fcit != fstart); 1929 1930 Point_3 orig(0,0,0); 1931 for(unsigned int i=1; i<pts.size()-1; i++) 1932 vol += Tetrahedron(orig,pts[0],pts[i],pts[i+1]).volume(); 1933 } 1934 return vol; 1935 } 1936 1937 /// Centroid computations 1938 1939 // Note: Centroid computation for polygons requires to evaluate 1940 // square roots and thus cannot be done without changing the 1941 // Traits concept. 1942 1943 template <class ConstructCircumcenter> 1944 Point_3 dual_centroid(Vertex_handle v, ConstructCircumcenter construct_circumcenter) const 1945 { 1946 std::list<Edge> edges; 1947 incident_edges(v, std::back_inserter(edges)); 1948 1949 FT vol(0); 1950 FT x(0), y(0), z(0); 1951 for(typename std::list<Edge>::iterator eit = edges.begin(); 1952 eit != edges.end(); ++eit) { 1953 1954 // compute the dual of the edge *eit but handle the translations 1955 // with respect to the dual of v. That is why we cannot use one 1956 // of the existing dual functions here. 1957 Facet_circulator fstart = incident_facets(*eit); 1958 Facet_circulator fcit = fstart; 1959 std::vector<Point_3> pts; 1960 do { 1961 // TODO: possible speed-up by caching the circumcenters 1962 Point_3 dual_orig = periodic_circumcenter(fcit->first, construct_circumcenter).first; 1963 int idx = fcit->first->index(v); 1964 Offset off = periodic_point(fcit->first,idx).second; 1965 pts.push_back(construct_point(std::make_pair(dual_orig,-off))); 1966 ++fcit; 1967 } while(fcit != fstart); 1968 1969 Point_3 orig(0,0,0); 1970 FT tetvol; 1971 for(unsigned int i=1; i<pts.size()-1; i++) { 1972 tetvol = Tetrahedron(orig,pts[0],pts[i],pts[i+1]).volume(); 1973 x += (pts[0].x() + pts[i].x() + pts[i+1].x()) * tetvol; 1974 y += (pts[0].y() + pts[i].y() + pts[i+1].y()) * tetvol; 1975 z += (pts[0].z() + pts[i].z() + pts[i+1].z()) * tetvol; 1976 vol += tetvol; 1977 } 1978 } 1979 x /= ( 4 * vol ); 1980 y /= ( 4 * vol ); 1981 z /= ( 4 * vol ); 1982 1983 Iso_cuboid d = domain(); 1984 x = (x < d.xmin() ? x+d.xmax()-d.xmin() 1985 : (x >= d.xmax() ? x-d.xmax()+d.xmin() : x)); 1986 y = (y < d.ymin() ? y+d.ymax()-d.ymin() 1987 : (y >= d.ymax() ? y-d.ymax()+d.ymin() : y)); 1988 z = (z < d.zmin() ? z+d.zmax()-d.zmin() 1989 : (z >= d.zmax() ? z-d.zmax()+d.zmin() : z)); 1990 1991 CGAL_triangulation_postcondition((d.xmin()<=x)&&(x<d.xmax())); 1992 CGAL_triangulation_postcondition((d.ymin()<=y)&&(y<d.ymax())); 1993 CGAL_triangulation_postcondition((d.zmin()<=z)&&(z<d.zmax())); 1994 1995 return Point_3(x,y,z); 1996 } 1997 }; 1998 1999 template < class GT, class TDS > 2000 inline bool 2001 Periodic_3_triangulation_3<GT,TDS>:: 2002 is_triangulation_in_1_sheet() const 2003 { 2004 if(is_1_cover()) 2005 return true; 2006 2007 for(Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) { 2008 if(virtual_vertices.find(vit) == virtual_vertices.end()) 2009 continue; 2010 2011 std::vector<Vertex_handle> nb_v; 2012 std::set<Vertex_handle> nb_v_odom; 2013 Vertex_handle vh; 2014 Offset off; 2015 adjacent_vertices(vit, std::back_inserter(nb_v)); 2016 for(unsigned int i=0; i<nb_v.size(); i++) { 2017 get_vertex(nb_v[i],vh,off); 2018 nb_v_odom.insert(vh); 2019 } 2020 if(nb_v.size() != nb_v_odom.size()) 2021 return false; 2022 } 2023 return true; 2024 } 2025 2026 template < class GT, class TDS > 2027 inline bool 2028 Periodic_3_triangulation_3<GT,TDS>:: 2029 is_vertex(const Point& p, Vertex_handle& v) const 2030 { 2031 Locate_type lt; 2032 int li, lj; 2033 Cell_handle c = locate( p, lt, li, lj ); 2034 if( lt != VERTEX ) 2035 return false; 2036 v = c->vertex(li); 2037 return true; 2038 } 2039 2040 template < class GT, class TDS > 2041 inline void 2042 Periodic_3_triangulation_3<GT,TDS>:: 2043 make_canonical(Vertex_triple& t) const 2044 { 2045 int i = (t.first < t.second) ? 0 : 1; 2046 if(i==0) { 2047 i = (t.first < t.third) ? 0 : 2; 2048 } else { 2049 i = (t.second < t.third) ? 1 : 2; 2050 } 2051 Vertex_handle tmp; 2052 switch(i) { 2053 case 0: return; 2054 case 1: 2055 tmp = t.first; 2056 t.first = t.second; 2057 t.second = t.third; 2058 t.third = tmp; 2059 return; 2060 default: 2061 tmp = t.first; 2062 t.first = t.third; 2063 t.third = t.second; 2064 t.second = tmp; 2065 } 2066 } 2067 2068 /** Assumes a point, an offset, and a cell to start from. 2069 * Gives the locate type and the simplex (cell and indices) containing p. 2070 * 2071 * Performs a remembering stochastic walk if the triangulation is not empty. 2072 * After the walk the type of the simplex containing p is determined. 2073 * 2074 * returns the cell p lies in 2075 * starts at cell "start" 2076 * returns a cell Cell_handel if lt == CELL 2077 * returns a facet (Cell_handle,li) if lt == FACET 2078 * returns an edge (Cell_handle,li,lj) if lt == EDGE 2079 * returns a vertex (Cell_handle,li) if lt == VERTEX 2080 */ 2081 template < class GT, class TDS > 2082 inline typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle 2083 Periodic_3_triangulation_3<GT,TDS>:: 2084 #ifdef CGAL_NO_STRUCTURAL_FILTERING 2085 periodic_locate 2086 #else 2087 exact_periodic_locate 2088 #endif 2089 (const Point& p, const Offset& o_p, Offset& lo, 2090 Locate_type& lt, int& li, int& lj, Cell_handle start) const 2091 { 2092 int cumm_off = 0; 2093 Offset off_query = o_p; 2094 if(number_of_vertices() == 0) { 2095 lo = Offset(); 2096 lt = EMPTY; 2097 return Cell_handle(); 2098 } 2099 CGAL_triangulation_assertion(number_of_vertices() != 0); 2100 2101 if(start == Cell_handle()) { 2102 start = cells_begin(); 2103 } 2104 2105 cumm_off = start->offset(0) | start->offset(1) 2106 | start->offset(2) | start->offset(3); 2107 if(is_1_cover() && cumm_off != 0) { 2108 if(((cumm_off & 4) == 4) && (FT(2)*p.x()<(domain().xmax()+domain().xmin()))) 2109 off_query += Offset(1,0,0); 2110 if(((cumm_off & 2) == 2) && (FT(2)*p.y()<(domain().ymax()+domain().ymin()))) 2111 off_query += Offset(0,1,0); 2112 if(((cumm_off & 1) == 1) && (FT(2)*p.z()<(domain().zmax()+domain().zmin()))) 2113 off_query += Offset(0,0,1); 2114 } 2115 2116 CGAL_triangulation_postcondition(start!=Cell_handle()); 2117 CGAL_triangulation_assertion(start->neighbor(0)->neighbor( 2118 start->neighbor(0)->index(start))==start); 2119 CGAL_triangulation_assertion(start->neighbor(1)->neighbor( 2120 start->neighbor(1)->index(start))==start); 2121 CGAL_triangulation_assertion(start->neighbor(2)->neighbor( 2122 start->neighbor(2)->index(start))==start); 2123 CGAL_triangulation_assertion(start->neighbor(3)->neighbor( 2124 start->neighbor(3)->index(start))==start); 2125 2126 // We implement the remembering visibility/stochastic walk. 2127 2128 // Remembers the previous cell to avoid useless orientation tests. 2129 Cell_handle previous = Cell_handle(); 2130 Cell_handle c = start; 2131 2132 // Stores the results of the 4 orientation tests. It will be used 2133 // at the end to decide if p lies on a face/edge/vertex/interior. 2134 Orientation o[4]; 2135 2136 boost::rand48 rng; 2137 boost::uniform_smallint<> four(0, 3); 2138 boost::variate_generator<boost::rand48&, boost::uniform_smallint<> > die4(rng, four); 2139 2140 // Now treat the cell c. 2141 try_next_cell: 2142 // For the remembering stochastic walk, 2143 // we need to start trying with a random index : 2144 int i = die4(); 2145 // For the remembering visibility walk (Delaunay only), we don't : 2146 // int i = 0; 2147 2148 cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); 2149 2150 bool simplicity_criterion = (cumm_off == 0) && (off_query.is_null()); 2151 2152 // We know that the 4 vertices of c are positively oriented. 2153 // So, in order to test if p is seen outside from one of c's facets, 2154 // we just replace the corresponding point by p in the orientation 2155 // test. We do this using the arrays below. 2156 2157 Offset off[4]; 2158 const Point* pts[4] = { &(c->vertex(0)->point()), 2159 &(c->vertex(1)->point()), 2160 &(c->vertex(2)->point()), 2161 &(c->vertex(3)->point()) }; 2162 2163 if(!simplicity_criterion && is_1_cover() ) { 2164 for(int i=0; i<4; i++) { 2165 off[i] = int_to_off(c->offset(i)); 2166 } 2167 } 2168 2169 if(!is_1_cover()) { 2170 // Just fetch the vertices of c as points with offsets 2171 for(int i=0; i<4; i++) { 2172 pts[i] = &(c->vertex(i)->point()); 2173 off[i] = get_offset(c,i); 2174 } 2175 } 2176 2177 for(int j=0; j != 4; ++j, i = (i+1)&3) { 2178 Cell_handle next = c->neighbor(i); 2179 if(previous == next) { 2180 o[i] = POSITIVE; 2181 continue; 2182 } 2183 2184 CGAL_triangulation_assertion(next->neighbor(next->index(c)) == c); 2185 2186 // We temporarily put p at i's place in pts. 2187 const Point* backup = pts[i]; 2188 pts[i] = &p; 2189 2190 if(simplicity_criterion && is_1_cover() ) { 2191 o[i] = orientation(*pts[0], *pts[1], *pts[2], *pts[3]); 2192 2193 if( o[i] != NEGATIVE ) { 2194 pts[i] = backup; 2195 continue; 2196 } 2197 } 2198 else { 2199 Offset backup_off; 2200 2201 backup_off = off[i]; 2202 off[i] = off_query; 2203 o[i] = orientation(*pts[0], *pts[1], *pts[2], *pts[3], 2204 off[0], off[1], off[2], off[3]); 2205 2206 if( o[i] != NEGATIVE ) { 2207 pts[i] = backup; 2208 off[i] = backup_off; 2209 continue; 2210 } 2211 } 2212 2213 // Test whether we need to adapt the offset of the query point. 2214 // This means, if we get out of the current cover. 2215 off_query = combine_offsets(off_query, neighbor_offset(c,i,next)); 2216 previous = c; 2217 c = next; 2218 goto try_next_cell; 2219 } 2220 2221 2222 // Ok, now we have found the cell. It remains to find the dimension of the 2223 // intersected simplex. 2224 // now p is in c or on its boundary 2225 int sum = ( o[0] == COPLANAR ) 2226 + ( o[1] == COPLANAR ) 2227 + ( o[2] == COPLANAR ) 2228 + ( o[3] == COPLANAR ); 2229 switch(sum) { 2230 case 0: 2231 lt = CELL; 2232 break; 2233 case 1: 2234 lt = FACET; 2235 li = ( o[0] == COPLANAR ) ? 0 : 2236 ( o[1] == COPLANAR ) ? 1 : 2237 ( o[2] == COPLANAR ) ? 2 : 3; 2238 break; 2239 case 2: 2240 lt = EDGE; 2241 li = ( o[0] != COPLANAR ) ? 0 : 2242 ( o[1] != COPLANAR ) ? 1 : 2; 2243 lj = ( o[li+1] != COPLANAR ) ? li+1 : 2244 ( o[li+2] != COPLANAR ) ? li+2 : li+3; 2245 break; 2246 case 3: 2247 lt = VERTEX; 2248 li = ( o[0] != COPLANAR ) ? 0 : 2249 ( o[1] != COPLANAR ) ? 1 : 2250 ( o[2] != COPLANAR ) ? 2 : 3; 2251 break; 2252 default: 2253 // the point cannot lie on four facets 2254 CGAL_triangulation_assertion(false); 2255 } 2256 lo = off_query; 2257 2258 return c; 2259 } 2260 2261 2262 #ifndef CGAL_NO_STRUCTURAL_FILTERING 2263 template < class GT, class TDS > 2264 typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle 2265 Periodic_3_triangulation_3<GT,TDS>:: 2266 inexact_periodic_locate(const Point& p, const Offset& o_p, 2267 Cell_handle start, 2268 int n_of_turns) const 2269 { 2270 int cumm_off = 0; 2271 Offset off_query = o_p; 2272 if(number_of_vertices() == 0) { 2273 return Cell_handle(); 2274 } 2275 CGAL_triangulation_assertion(number_of_vertices() != 0); 2276 2277 if(start == Cell_handle()) { 2278 start = cells_begin(); 2279 } 2280 2281 cumm_off = start->offset(0) | start->offset(1) 2282 | start->offset(2) | start->offset(3); 2283 if(is_1_cover() && cumm_off != 0) { 2284 if(((cumm_off & 4) == 4) && (FT(2)*p.x()<(domain().xmax()+domain().xmin()))) 2285 off_query += Offset(1,0,0); 2286 if(((cumm_off & 2) == 2) && (FT(2)*p.y()<(domain().ymax()+domain().ymin()))) 2287 off_query += Offset(0,1,0); 2288 if(((cumm_off & 1) == 1) && (FT(2)*p.z()<(domain().zmax()+domain().zmin()))) 2289 off_query += Offset(0,0,1); 2290 } 2291 2292 CGAL_triangulation_postcondition(start!=Cell_handle()); 2293 CGAL_triangulation_assertion(start->neighbor(0)->neighbor( 2294 start->neighbor(0)->index(start))==start); 2295 CGAL_triangulation_assertion(start->neighbor(1)->neighbor( 2296 start->neighbor(1)->index(start))==start); 2297 CGAL_triangulation_assertion(start->neighbor(2)->neighbor( 2298 start->neighbor(2)->index(start))==start); 2299 CGAL_triangulation_assertion(start->neighbor(3)->neighbor( 2300 start->neighbor(3)->index(start))==start); 2301 2302 // We implement the remembering visibility/stochastic walk. 2303 2304 // Remembers the previous cell to avoid useless orientation tests. 2305 Cell_handle previous = Cell_handle(); 2306 Cell_handle c = start; 2307 2308 // Now treat the cell c. 2309 try_next_cell: 2310 --n_of_turns; 2311 cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); 2312 2313 bool simplicity_criterion = (cumm_off == 0) && (off_query.is_null()); 2314 2315 // We know that the 4 vertices of c are positively oriented. 2316 // So, in order to test if p is seen outside from one of c's facets, 2317 // we just replace the corresponding point by p in the orientation 2318 // test. We do this using the arrays below. 2319 2320 Offset off[4]; 2321 const Point* pts[4] = { &(c->vertex(0)->point()), 2322 &(c->vertex(1)->point()), 2323 &(c->vertex(2)->point()), 2324 &(c->vertex(3)->point()) }; 2325 2326 if(!simplicity_criterion && is_1_cover() ) { 2327 for(int i=0; i<4; i++) { 2328 off[i] = int_to_off(c->offset(i)); 2329 } 2330 } 2331 2332 if(!is_1_cover()) { 2333 // Just fetch the vertices of c as points with offsets 2334 for(int i=0; i<4; i++) { 2335 pts[i] = &(c->vertex(i)->point()); 2336 off[i] = get_offset(c,i); 2337 } 2338 } 2339 2340 for(int i=0; i != 4; ++i) { 2341 Cell_handle next = c->neighbor(i); 2342 if(previous == next) { 2343 continue; 2344 } 2345 2346 // We temporarily put p at i's place in pts. 2347 const Point* backup = pts[i]; 2348 pts[i] = &p; 2349 2350 if(simplicity_criterion && is_1_cover() ) { 2351 if( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3]) != NEGATIVE ) { 2352 pts[i] = backup; 2353 continue; 2354 } 2355 } 2356 else { 2357 Offset backup_off; 2358 2359 backup_off = off[i]; 2360 off[i] = off_query; 2361 2362 if( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3], 2363 off[0], off[1], off[2], off[3]) != NEGATIVE ) { 2364 pts[i] = backup; 2365 off[i] = backup_off; 2366 continue; 2367 } 2368 } 2369 2370 // Test whether we need to adapt the offset of the query point. 2371 // This means, if we get out of the current cover. 2372 off_query = combine_offsets(off_query, neighbor_offset(c,i,next)); 2373 previous = c; 2374 c = next; 2375 if(n_of_turns) 2376 goto try_next_cell; 2377 } 2378 2379 return c; 2380 } 2381 #endif 2382 2383 2384 /** 2385 * returns 2386 * ON_BOUNDED_SIDE if (q, off) inside the cell 2387 * ON_BOUNDARY if (q, off) on the boundary of the cell 2388 * ON_UNBOUNDED_SIDE if (q, off) lies outside the cell 2389 * 2390 * lt has a meaning only when ON_BOUNDED_SIDE or ON_BOUNDARY 2391 */ 2392 template < class GT, class TDS > 2393 inline Bounded_side Periodic_3_triangulation_3<GT,TDS>::side_of_cell( 2394 const Point& q, const Offset& off, Cell_handle c, 2395 Locate_type& lt, int& i, int& j) const 2396 { 2397 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 2398 2399 Orientation o0,o1,o2,o3; 2400 o0 = o1 = o2 = o3 = ZERO; 2401 2402 int cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); 2403 if((cumm_off == 0) && (is_1_cover())) { 2404 CGAL_triangulation_assertion(off == Offset()); 2405 const Point& p0 = c->vertex(0)->point(); 2406 const Point& p1 = c->vertex(1)->point(); 2407 const Point& p2 = c->vertex(2)->point(); 2408 const Point& p3 = c->vertex(3)->point(); 2409 2410 if(((o0 = orientation(q ,p1,p2,p3)) == NEGATIVE) || 2411 ((o1 = orientation(p0,q ,p2,p3)) == NEGATIVE) || 2412 ((o2 = orientation(p0,p1,q ,p3)) == NEGATIVE) || 2413 ((o3 = orientation(p0,p1,p2,q )) == NEGATIVE) ) { 2414 return ON_UNBOUNDED_SIDE; 2415 } 2416 } else { // Special case for the periodic space. 2417 Offset off_q; 2418 Offset offs[4]; 2419 const Point *p[4]; 2420 for(int i=0; i<4; i++) { 2421 p[i] = &(c->vertex(i)->point()); 2422 offs[i] = get_offset(c,i); 2423 } 2424 CGAL_triangulation_assertion(orientation(*p[0], *p[1], *p[2], *p[3], 2425 offs[0], offs[1], offs[2], offs[3]) == POSITIVE); 2426 bool found=false; 2427 for(int i=0; (i<8)&&(!found); i++) { 2428 if((cumm_off | ((~i)&7)) == 7) { 2429 o0 = o1 = o2 = o3 = NEGATIVE; 2430 off_q = combine_offsets(off, int_to_off(i)); 2431 2432 if(((o0 = orientation(q, *p[1], *p[2], *p[3], 2433 off_q ,offs[1],offs[2],offs[3])) != NEGATIVE)&& 2434 ((o1 = orientation(*p[0], q, *p[2], *p[3], 2435 offs[0], off_q,offs[2],offs[3])) != NEGATIVE)&& 2436 ((o2 = orientation(*p[0], *p[1], q, *p[3], 2437 offs[0],offs[1], off_q,offs[3])) != NEGATIVE)&& 2438 ((o3 = orientation(*p[0], *p[1], *p[2], q, 2439 offs[0],offs[1],offs[2], off_q)) != NEGATIVE)) { 2440 found = true; 2441 } 2442 } 2443 } 2444 if(!found) return ON_UNBOUNDED_SIDE; 2445 } 2446 2447 // now all the oi's are >=0 2448 // sum gives the number of facets p lies on 2449 int sum = ( (o0 == ZERO) ? 1 : 0 ) 2450 + ( (o1 == ZERO) ? 1 : 0 ) 2451 + ( (o2 == ZERO) ? 1 : 0 ) 2452 + ( (o3 == ZERO) ? 1 : 0 ); 2453 2454 switch(sum) { 2455 case 0: 2456 { 2457 lt = CELL; 2458 return ON_BOUNDED_SIDE; 2459 } 2460 case 1: 2461 { 2462 lt = FACET; 2463 // i = index such that p lies on facet(i) 2464 i = ( o0 == ZERO ) ? 0 : 2465 ( o1 == ZERO ) ? 1 : 2466 ( o2 == ZERO ) ? 2 : 3; 2467 return ON_BOUNDARY; 2468 } 2469 case 2: 2470 { 2471 lt = EDGE; 2472 // i = smallest index such that p does not lie on facet(i) 2473 // i must be < 3 since p lies on 2 facets 2474 i = ( o0 == POSITIVE ) ? 0 : 2475 ( o1 == POSITIVE ) ? 1 : 2; 2476 // j = larger index such that p not on facet(j) 2477 // j must be > 0 since p lies on 2 facets 2478 j = ( o3 == POSITIVE ) ? 3 : 2479 ( o2 == POSITIVE ) ? 2 : 1; 2480 return ON_BOUNDARY; 2481 } 2482 case 3: 2483 { 2484 lt = VERTEX; 2485 // i = index such that p does not lie on facet(i) 2486 i = ( o0 == POSITIVE ) ? 0 : 2487 ( o1 == POSITIVE ) ? 1 : 2488 ( o2 == POSITIVE ) ? 2 : 3; 2489 return ON_BOUNDARY; 2490 } 2491 default: 2492 { 2493 // impossible : cannot be on 4 facets for a real tetrahedron 2494 CGAL_triangulation_assertion(false); 2495 return ON_BOUNDARY; 2496 } 2497 } 2498 } // side_of_cell 2499 2500 /*! \brief inserts point. 2501 * 2502 * Inserts the point p into the triangulation. It assumes that 2503 * the cell c containing p is already known. 2504 * 2505 * Implementation: 2506 * - some precondition checking 2507 * - find and mark conflicting cells --> find_conflicts 2508 * Conflicting cells are stored in the vector cells. 2509 * - backup hidden points 2510 * - Delete the edges of the marked cells from too_long_edges 2511 * - Insert the new vertex in the hole obtained by removing the 2512 * conflicting cells (star-approach) --> _tds._insert_in_hole 2513 * - find out about offsets 2514 * - Insert the newly added edges that are "too long" 2515 * to too_long_edges 2516 * - reinsert hidden points 2517 */ 2518 template < class GT, class TDS > 2519 template < class Conflict_tester, class Point_hider, class CoverManager > 2520 inline typename Periodic_3_triangulation_3<GT,TDS>::Vertex_handle 2521 Periodic_3_triangulation_3<GT,TDS>::periodic_insert( 2522 const Point& p, const Offset& o, 2523 Locate_type /*lt*/, Cell_handle c, const Conflict_tester& tester, 2524 Point_hider& hider, CoverManager& cover_manager, Vertex_handle vh) 2525 { 2526 Vertex_handle v; 2527 CGAL_triangulation_precondition(number_of_vertices() != 0); 2528 CGAL_triangulation_assertion_code( 2529 Locate_type lt_assert; int i_assert; int j_assert;); 2530 CGAL_triangulation_precondition(side_of_cell(tester.point(),o, c, 2531 lt_assert, i_assert, j_assert) != ON_UNBOUNDED_SIDE); 2532 2533 tester.set_offset(o); 2534 2535 // Choose the periodic copy of tester.point() that is in conflict with c. 2536 bool found = false; 2537 Offset current_off = get_location_offset(tester, c, found); 2538 2539 CGAL_triangulation_assertion(side_of_cell(tester.point(), 2540 combine_offsets(o,current_off),c,lt_assert,i_assert,j_assert) 2541 != ON_UNBOUNDED_SIDE); 2542 2543 // If the new point is not in conflict with its cell, it is hidden. 2544 if(!found || !tester.test_initial_cell(c, current_off)) { 2545 hider.hide_point(c,p); 2546 return Vertex_handle(); 2547 } 2548 2549 // Ok, we really insert the point now. 2550 // First, find the conflict region. 2551 std::vector<Cell_handle> cells; 2552 cells.reserve(32); 2553 2554 Facet facet; 2555 2556 find_conflicts(c, current_off, tester, 2557 make_triple(Oneset_iterator<Facet>(facet), 2558 std::back_inserter(cells), 2559 Emptyset_iterator())); 2560 2561 // Remember the points that are hidden by the conflicting cells, 2562 // as they will be deleted during the insertion. 2563 hider.set_vertices(cells.begin(), cells.end()); 2564 2565 if(!is_1_cover()) 2566 cover_manager.delete_unsatisfying_elements(cells.begin(), cells.end()); 2567 2568 // Insertion. Warning: facets[0].first MUST be in conflict! 2569 // Compute the star and put it into the data structure. 2570 // Store the new cells from the star in nbs. 2571 v = _tds._insert_in_hole(cells.begin(), cells.end(), facet.first, facet.second); 2572 v->set_point(p); 2573 2574 //TODO: this could be done within the _insert_in_hole without losing any 2575 //time because each cell is visited in any case. 2576 //- Do timings to argue to modify _insert_in_conflicts if need be 2577 //- Find the modified _insert_in_hole in the branch svn history of TDS 2578 std::vector<Cell_handle> nbs; 2579 incident_cells(v, std::back_inserter(nbs)); 2580 // For all neighbors of the newly added vertex v: fetch their offsets from 2581 // the tester and reset them in the triangulation data structure. 2582 for(typename std::vector<Cell_handle>::iterator cit = nbs.begin(); 2583 cit != nbs.end(); cit++) { 2584 Offset off[4]; 2585 for(int i=0; i<4; i++) { 2586 off[i] = (*cit)->vertex(i)->offset(); 2587 } 2588 set_offsets(*cit, off[0], off[1], off[2], off[3]); 2589 } 2590 2591 for(typename std::vector<Vertex_handle>::iterator voit = v_offsets.begin(); 2592 voit != v_offsets.end(); ++voit) { 2593 (*voit)->clear_offset(); 2594 } 2595 v_offsets.clear(); 2596 2597 if(vh != Vertex_handle()) { 2598 // CGAL_triangulation_assertion(virtual_vertices.find(v) == virtual_vertices.end()); 2599 virtual_vertices[v] = Virtual_vertex(vh,o); 2600 virtual_vertices_reverse[vh].push_back(v); 2601 } 2602 2603 if(!is_1_cover()) 2604 cover_manager.insert_unsatisfying_elements(v, nbs.begin(), nbs.end()); 2605 2606 // Store the hidden points in their new cells. 2607 hider.reinsert_vertices(v); 2608 return v; 2609 } 2610 2611 /** Inserts the first point to a triangulation. 2612 * 2613 * With inserting the first point the 3-sheeted covering is constructed. 2614 * So first, the 27 vertices are inserted and are added to virtual_vertices 2615 * Then 6*27 cells are created. 2616 * Then all links are set. 2617 */ 2618 template < class GT, class TDS > 2619 inline typename Periodic_3_triangulation_3<GT,TDS>::Vertex_handle 2620 Periodic_3_triangulation_3<GT,TDS>::create_initial_triangulation(const Point& p) 2621 { 2622 /// Virtual vertices, one per periodic instance 2623 Vertex_handle vir_vertices[3][3][3]; 2624 /// Virtual cells, 6 per periodic instance 2625 Cell_handle cells[3][3][3][6]; 2626 2627 // Initialise vertices: 2628 vir_vertices[0][0][0] = _tds.create_vertex(); 2629 vir_vertices[0][0][0]->set_point(p); 2630 virtual_vertices_reverse[vir_vertices[0][0][0]] = std::vector<Vertex_handle>(); 2631 for(int i=0; i<_cover[0]; i++) { 2632 for(int j=0; j<_cover[1]; j++) { 2633 for(int k=0; k<_cover[2]; k++) { 2634 if((i!=0)||(j!=0)||(k!=0)) { 2635 // Initialise virtual vertices out of the domain for debugging 2636 vir_vertices[i][j][k] = 2637 _tds.create_vertex(); 2638 vir_vertices[i][j][k]->set_point(p); //+Offset(i,j,k)); 2639 virtual_vertices[vir_vertices[i][j][k]] = 2640 Virtual_vertex(vir_vertices[0][0][0], Offset(i,j,k)); 2641 virtual_vertices_reverse[vir_vertices[0][0][0]].push_back( 2642 vir_vertices[i][j][k]); 2643 } 2644 CGAL_triangulation_assertion(vir_vertices[i][j][k] != Vertex_handle()); 2645 CGAL_triangulation_assertion(vir_vertices[0][0][0]->point() == p); 2646 } 2647 } 2648 } 2649 2650 // Create cells: 2651 for(int i=0; i<_cover[0]; i++) { 2652 for(int j=0; j<_cover[1]; j++) { 2653 for(int k=0; k<_cover[2]; k++) { 2654 for(int l=0; l<6; l++) { 2655 // 6 cells per 'cube' 2656 cells[i][j][k][l] = _tds.create_cell(); 2657 for(int n=0; n<4; n++) 2658 CGAL_triangulation_assertion(cells[i][j][k][l] != Cell_handle()); 2659 } 2660 } 2661 } 2662 } 2663 // set vertex and neighbor information 2664 // index to the right vertex: [number of cells][vertex][offset] 2665 int vertex_ind[6][4][3] = { 2666 { {0, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 0, 0} }, 2667 { {1, 1, 0}, {0, 1, 1}, {1, 0, 1}, {1, 1, 1} }, 2668 { {1, 0, 0}, {0, 1, 1}, {0, 1, 0}, {0, 0, 1} }, 2669 { {1, 0, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1} }, 2670 { {1, 0, 0}, {0, 1, 1}, {1, 0, 1}, {1, 1, 0} }, 2671 { {1, 0, 0}, {0, 1, 1}, {1, 1, 0}, {0, 1, 0} } 2672 }; 2673 int neighb_ind[6][4][4] = { 2674 { { 0, 0, 0, 2}, { 0,-1, 0, 5}, { 0, 0,-1, 3}, {-1, 0, 0, 4} }, 2675 { { 0, 0, 1, 5}, { 1, 0, 0, 2}, { 0, 1, 0, 3}, { 0, 0, 0, 4} }, 2676 { {-1, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 3}, { 0, 0, 0, 5} }, 2677 { { 0, 0, 1, 0}, { 0,-1, 0, 1}, { 0, 0, 0, 4}, { 0, 0, 0, 2} }, 2678 { { 0, 0, 0, 1}, { 1, 0, 0, 0}, { 0, 0, 0, 5}, { 0, 0, 0, 3} }, 2679 { { 0, 1, 0, 0}, { 0, 0,-1, 1}, { 0, 0, 0, 2}, { 0, 0, 0, 4} } 2680 }; 2681 for(int i=0; i<_cover[0]; i++) { 2682 for(int j=0; j<_cover[1]; j++) { 2683 for(int k=0; k<_cover[2]; k++) { 2684 int offset = (i==_cover[0]-1 ? 4 : 0) | 2685 (j==_cover[1]-1 ? 2 : 0) | 2686 (k==_cover[2]-1 ? 1 : 0); 2687 for(int l=0; l<6; l++) { 2688 // cell 0: 2689 cells[i][j][k][l]->set_vertices( 2690 vir_vertices 2691 [(i+vertex_ind[l][0][0])%_cover[0]] 2692 [(j+vertex_ind[l][0][1])%_cover[1]] 2693 [(k+vertex_ind[l][0][2])%_cover[2]], 2694 vir_vertices 2695 [(i+vertex_ind[l][1][0])%_cover[0]] 2696 [(j+vertex_ind[l][1][1])%_cover[1]] 2697 [(k+vertex_ind[l][1][2])%_cover[2]], 2698 vir_vertices 2699 [(i+vertex_ind[l][2][0])%_cover[0]] 2700 [(j+vertex_ind[l][2][1])%_cover[1]] 2701 [(k+vertex_ind[l][2][2])%_cover[2]], 2702 vir_vertices 2703 [(i+vertex_ind[l][3][0])%_cover[0]] 2704 [(j+vertex_ind[l][3][1])%_cover[1]] 2705 [(k+vertex_ind[l][3][2])%_cover[2]]); 2706 set_offsets(cells[i][j][k][l], 2707 offset & (vertex_ind[l][0][0]*4 + 2708 vertex_ind[l][0][1]*2 + 2709 vertex_ind[l][0][2]*1), 2710 offset & (vertex_ind[l][1][0]*4 + 2711 vertex_ind[l][1][1]*2 + 2712 vertex_ind[l][1][2]*1), 2713 offset & (vertex_ind[l][2][0]*4 + 2714 vertex_ind[l][2][1]*2 + 2715 vertex_ind[l][2][2]*1), 2716 offset & (vertex_ind[l][3][0]*4 + 2717 vertex_ind[l][3][1]*2 + 2718 vertex_ind[l][3][2]*1)); 2719 cells[i][j][k][l]->set_neighbors( 2720 cells [(i+_cover[0]+neighb_ind[l][0][0])%_cover[0]] 2721 [(j+_cover[1]+neighb_ind[l][0][1])%_cover[1]] 2722 [(k+_cover[2]+neighb_ind[l][0][2])%_cover[2]] 2723 [ neighb_ind[l][0][3] ], 2724 cells [(i+_cover[0]+neighb_ind[l][1][0])%_cover[0]] 2725 [(j+_cover[1]+neighb_ind[l][1][1])%_cover[1]] 2726 [(k+_cover[2]+neighb_ind[l][1][2])%_cover[2]] 2727 [ neighb_ind[l][1][3] ], 2728 cells [(i+_cover[0]+neighb_ind[l][2][0])%_cover[0]] 2729 [(j+_cover[1]+neighb_ind[l][2][1])%_cover[1]] 2730 [(k+_cover[2]+neighb_ind[l][2][2])%_cover[2]] 2731 [ neighb_ind[l][2][3] ], 2732 cells [(i+_cover[0]+neighb_ind[l][3][0])%_cover[0]] 2733 [(j+_cover[1]+neighb_ind[l][3][1])%_cover[1]] 2734 [(k+_cover[2]+neighb_ind[l][3][2])%_cover[2]] 2735 [ neighb_ind[l][3][3] ] 2736 ); 2737 } 2738 } 2739 } 2740 } 2741 // set pointers from the vertices to incident cells. 2742 for(int i=0; i<_cover[0]; i++) { 2743 for(int j=0; j<_cover[1]; j++) { 2744 for(int k=0; k<_cover[2]; k++) { 2745 vir_vertices[i][j][k]->set_cell(cells[i][j][k][0]); 2746 } 2747 } 2748 } 2749 2750 _tds.set_dimension(3); 2751 2752 return vir_vertices[0][0][0]; 2753 } 2754 #define CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H 2755 #include <CGAL/internal/Periodic_3_triangulation_dummy_36.h> 2756 #undef CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H 2757 2758 /** finds all cells that are in conflict with the currently added point 2759 * (stored in tester). 2760 * 2761 * The result will be a hole of which the following data is returned: 2762 * - boundary facets 2763 * - cells 2764 * - internal facets. 2765 * 2766 * c is the current cell, which must be in conflict. 2767 * tester is the function object that tests if a cell is in conflict. 2768 */ 2769 template <class GT, class TDS> 2770 template <class Conflict_test, 2771 class OutputIteratorBoundaryFacets, 2772 class OutputIteratorCells, 2773 class OutputIteratorInternalFacets> 2774 Triple<OutputIteratorBoundaryFacets, 2775 OutputIteratorCells, 2776 OutputIteratorInternalFacets> 2777 Periodic_3_triangulation_3<GT,TDS>:: 2778 find_conflicts(Cell_handle d, const Offset& current_off, 2779 const Conflict_test& tester, 2780 Triple<OutputIteratorBoundaryFacets, 2781 OutputIteratorCells, 2782 OutputIteratorInternalFacets> it) const 2783 { 2784 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 2785 CGAL_triangulation_precondition( tester(d, current_off) ); 2786 2787 std::stack<std::pair<Cell_handle, Offset> > cell_stack; 2788 cell_stack.push(std::make_pair(d,current_off)); 2789 d->tds_data().mark_in_conflict(); 2790 *it.second++ = d; 2791 2792 do { 2793 Cell_handle c = cell_stack.top().first; 2794 Offset current_off2 = cell_stack.top().second; 2795 cell_stack.pop(); 2796 2797 for(int i=0; i< 4; ++i) { 2798 Cell_handle test = c->neighbor(i); 2799 if(test->tds_data().is_in_conflict()) { 2800 if(c < test) { 2801 *it.third++ = Facet(c, i); // Internal facet. 2802 } 2803 continue; // test was already in conflict. 2804 } 2805 if(test->tds_data().is_clear()) { 2806 Offset o_test = current_off2 + neighbor_offset(c, i, test); 2807 if(tester(test,o_test)) { 2808 if(c < test) 2809 *it.third++ = Facet(c, i); // Internal facet. 2810 2811 cell_stack.push(std::make_pair(test,o_test)); 2812 test->tds_data().mark_in_conflict(); 2813 *it.second++ = test; 2814 continue; 2815 } 2816 test->tds_data().mark_on_boundary(); // test is on the boundary. 2817 } 2818 *it.first++ = Facet(c, i); 2819 for(int j = 0; j<4; j++) { 2820 if(j==i) continue; 2821 if(!c->vertex(j)->get_offset_flag()) { 2822 c->vertex(j)->set_offset(int_to_off(c->offset(j))-current_off2); 2823 v_offsets.push_back(c->vertex(j)); 2824 } 2825 } 2826 } 2827 } while(!cell_stack.empty()); 2828 return it; 2829 } 2830 2831 /*! \brief inserts point into triangulation. 2832 * 2833 * Inserts the point p into the triangulation. It expects 2834 * - a cell to start the point location 2835 * - a testing function to determine cells in conflict 2836 * - a testing function to determine if a vertex is hidden. 2837 * 2838 * Implementation: 2839 * - If the triangulation is empty call a special function 2840 * (create_initial_triangulation) to construct the basic 2841 * triangulation. 2842 * - Run point location to get the cell c containing point p. 2843 * - Call periodic_insert to insert p into the 3-cover. 2844 * - Also insert the eight periodic copies of p. 2845 */ 2846 template < class GT, class TDS > 2847 template < class Conflict_tester, class Point_hider, class CoverManager > 2848 inline typename Periodic_3_triangulation_3<GT,TDS>::Vertex_handle 2849 Periodic_3_triangulation_3<GT,TDS>::insert_in_conflict(const Point& p, 2850 Locate_type lt, Cell_handle c, int li, int lj, 2851 const Conflict_tester& tester, Point_hider& hider, CoverManager& cover_manager) { 2852 CGAL_triangulation_assertion((domain().xmin() <= p.x()) 2853 && (p.x() < domain().xmax())); 2854 CGAL_triangulation_assertion((domain().ymin() <= p.y()) 2855 && (p.y() < domain().ymax())); 2856 CGAL_triangulation_assertion((domain().zmin() <= p.z()) 2857 && (p.z() < domain().zmax())); 2858 2859 if(number_of_vertices() == 0) { 2860 Vertex_handle vh = create_initial_triangulation(p); 2861 cover_manager.create_initial_triangulation(); 2862 return vh; 2863 } 2864 2865 if((lt == VERTEX) && 2866 (tester.compare_weight(c->vertex(li)->point(),p) == 0) ) { 2867 return c->vertex(li); 2868 } 2869 2870 Vertex_handle vstart; 2871 if(!is_1_cover()) { 2872 Virtual_vertex_map_it vvmit = virtual_vertices.find(c->vertex(0)); 2873 if(vvmit == virtual_vertices.end()) 2874 vstart = c->vertex(0); 2875 else 2876 vstart = vvmit->second.first; 2877 CGAL_triangulation_assertion(virtual_vertices.find(vstart) 2878 == virtual_vertices.end()); 2879 CGAL_triangulation_assertion(virtual_vertices_reverse.find(vstart) 2880 != virtual_vertices_reverse.end()); 2881 } 2882 2883 CGAL_triangulation_assertion( number_of_vertices() != 0 ); 2884 CGAL_triangulation_expensive_assertion(is_valid()); 2885 hider.set_original_cube(true); 2886 Vertex_handle vh = periodic_insert(p, Offset(), lt, c, tester, hider, cover_manager); 2887 if(is_1_cover()) { 2888 return vh; 2889 } 2890 2891 hider.set_original_cube(false); 2892 2893 virtual_vertices_reverse[vh] = std::vector<Vertex_handle>(); 2894 Offset lo; 2895 // insert 26 periodic copies 2896 for(int i=0; i<_cover[0]; i++) { 2897 for(int j=0; j<_cover[1]; j++) { 2898 for(int k=0; k<_cover[2]; k++) { 2899 if((i!=0)||(j!=0)||(k!=0)) { 2900 c = periodic_locate(p, Offset(i,j,k), lo, lt, li, lj, Cell_handle()); 2901 periodic_insert(p, Offset(i,j,k), lt, c, tester, hider,cover_manager,vh); 2902 } 2903 } 2904 } 2905 } 2906 CGAL_triangulation_expensive_assertion(is_valid()); 2907 2908 // Fall back to 1-cover if the criterion that the longest edge is shorter 2909 // than sqrt(0.166) is fulfilled. 2910 if( cover_manager.can_be_converted_to_1_sheet() ) { 2911 CGAL_triangulation_expensive_assertion(is_valid()); 2912 convert_to_1_sheeted_covering(); 2913 CGAL_triangulation_expensive_assertion( is_valid() ); 2914 } 2915 return vh; 2916 } 2917 2918 /// tests if two vertices of one cell are just periodic copies of each other 2919 template < class GT, class TDS > 2920 inline bool Periodic_3_triangulation_3<GT,TDS>::has_self_edges(Cell_handle c) const { 2921 CGAL_triangulation_assertion((c->vertex(0) != c->vertex(1)) || 2922 (c->offset(0) != c->offset(1))); 2923 CGAL_triangulation_assertion((c->vertex(0) != c->vertex(2)) || 2924 (c->offset(0) != c->offset(2))); 2925 CGAL_triangulation_assertion((c->vertex(0) != c->vertex(3)) || 2926 (c->offset(0) != c->offset(3))); 2927 CGAL_triangulation_assertion((c->vertex(1) != c->vertex(2)) || 2928 (c->offset(1) != c->offset(2))); 2929 CGAL_triangulation_assertion((c->vertex(1) != c->vertex(3)) || 2930 (c->offset(1) != c->offset(3))); 2931 CGAL_triangulation_assertion((c->vertex(2) != c->vertex(3)) || 2932 (c->offset(2) != c->offset(3))); 2933 return ((c->vertex(0) == c->vertex(1)) || 2934 (c->vertex(0) == c->vertex(2)) || 2935 (c->vertex(0) == c->vertex(3)) || 2936 (c->vertex(1) == c->vertex(2)) || 2937 (c->vertex(1) == c->vertex(3)) || 2938 (c->vertex(2) == c->vertex(3))); 2939 } 2940 2941 /*! \brief tests if the triangulation is valid. 2942 * 2943 * A triangulation is valid if 2944 * - A cell is not its own neighbor. 2945 * - A cell has no two equal neighbors 2946 * - A cell has no two equal vertex-offset pairs 2947 * - A cell is positively oriented. 2948 * - The point of a neighbor of cell c that does not belong to c is not inside 2949 * the circumcircle of c. 2950 */ 2951 template < class GT, class TDS > 2952 bool 2953 Periodic_3_triangulation_3<GT,TDS>:: 2954 is_valid(bool verbose, int level) const 2955 { 2956 if(!is_1_cover()) 2957 { 2958 for(Virtual_vertex_reverse_map_it iter = virtual_vertices_reverse.begin(), end_iter = virtual_vertices_reverse.end(); 2959 iter != end_iter; 2960 ++iter) 2961 { 2962 for(typename Virtual_vertex_reverse_map::mapped_type::const_iterator iter_2 = iter->second.begin(), 2963 end_iter_2 = iter->second.end(); 2964 iter_2 != end_iter_2; 2965 ++iter_2) 2966 { 2967 CGAL_triangulation_assertion(virtual_vertices.find(*iter_2) != virtual_vertices.end()); 2968 CGAL_triangulation_assertion(virtual_vertices.at(*iter_2).first == iter->first); 2969 } 2970 } 2971 } 2972 2973 bool error = false; 2974 for(Cell_iterator cit = cells_begin(); 2975 cit != cells_end(); ++cit) { 2976 for(int i=0; i<4; i++) { 2977 CGAL_triangulation_assertion(cit != cit->neighbor(i)); 2978 for(int j=i+1; j<4; j++) { 2979 CGAL_triangulation_assertion(cit->neighbor(i) != cit->neighbor(j)); 2980 CGAL_triangulation_assertion(cit->vertex(i) != cit->vertex(j)); 2981 } 2982 } 2983 // Check positive orientation: 2984 const Point *p[4]; Offset off[4]; 2985 for(int i=0; i<4; i++) { 2986 p[i] = &cit->vertex(i)->point(); 2987 off[i] = get_offset(cit,i); 2988 } 2989 if(orientation(*p[0], *p[1], *p[2], *p[3], 2990 off[0], off[1], off[2], off[3]) != POSITIVE) { 2991 if(verbose) { 2992 std::cerr<<"Periodic_3_triangulation_3: wrong orientation:"<<std::endl; 2993 std::cerr<<off[0]<<'\t'<<*p[0]<<'\n' 2994 <<off[1]<<'\t'<<*p[1]<<'\n' 2995 <<off[2]<<'\t'<<*p[2]<<'\n' 2996 <<off[3]<<'\t'<<*p[3]<<std::endl; 2997 } 2998 error = true; 2999 } 3000 } 3001 3002 if(!has_self_edges()) { 3003 if(! _tds.is_valid(verbose, level) ) { 3004 return false; 3005 } 3006 } 3007 3008 return !error; 3009 } 3010 3011 template < class GT, class TDS > 3012 bool Periodic_3_triangulation_3<GT,TDS>::is_valid(Cell_handle ch, 3013 bool verbose, int level) const 3014 { 3015 if( ! _tds.is_valid(ch,verbose,level) ) 3016 return false; 3017 3018 bool error = false; 3019 const Point *p[4]; Offset off[4]; 3020 for(int i=0; i<4; i++) { 3021 p[i] = &ch->vertex(i)->point(); 3022 off[i] = get_offset(ch,i); 3023 } 3024 3025 if(orientation(*p[0], *p[1], *p[2], *p[3], 3026 off[0], off[1], off[2], off[3]) != POSITIVE) { 3027 error = true; 3028 } 3029 3030 return !error; 3031 } 3032 3033 template < class GT, class TDS > 3034 template < class ConflictTester > 3035 bool Periodic_3_triangulation_3<GT,TDS>:: 3036 is_valid_conflict(ConflictTester& tester, bool verbose, int level) const 3037 { 3038 Cell_iterator it; 3039 for( it = cells_begin(); it != cells_end(); ++it ) { 3040 is_valid(it, verbose, level); 3041 for(int i=0; i<4; i++ ) { 3042 Offset o_nb = neighbor_offset(it, i, it->neighbor(i)); 3043 Offset o_vt = get_offset(it->neighbor(i), it->neighbor(i)->index(it)); 3044 if(tester(it, 3045 it->neighbor(i)->vertex(it->neighbor(i)->index(it))->point(), 3046 o_vt - o_nb)) { 3047 if(verbose) { 3048 std::cerr << "non-empty sphere:\n" 3049 << "Point[0]: " << it->vertex(0)->point() 3050 << " Off: " << int_to_off(it->offset(0)) << "\n" 3051 << "Point[1]: " << it->vertex(1)->point() 3052 << " Off: " << int_to_off(it->offset(1)) << "\n" 3053 << "Point[2]: " << it->vertex(2)->point() 3054 << " Off: " << int_to_off(it->offset(2)) << "\n" 3055 << "Point[3]: " << it->vertex(3)->point() 3056 << " Off: " << int_to_off(it->offset(3)) << "\n" 3057 << "Foreign point: " 3058 << it->neighbor(i)->vertex(it->neighbor(i)->index(it))->point() 3059 << "\t Off: " << o_vt - o_nb << std::endl; 3060 } 3061 return false; 3062 } 3063 } 3064 } 3065 return true; 3066 } 3067 3068 template < class GT, class TDS > 3069 inline void Periodic_3_triangulation_3<GT,TDS>::make_hole(Vertex_handle v, 3070 std::map<Vertex_triple,Facet>& outer_map, 3071 std::vector<Cell_handle>& hole) 3072 { 3073 incident_cells(v, std::back_inserter(hole)); 3074 3075 for(typename std::vector<Cell_handle>::iterator cit = hole.begin(); 3076 cit != hole.end(); ++cit) { 3077 int indv = (*cit)->index(v); 3078 Cell_handle opp_cit = (*cit)->neighbor( indv ); 3079 Facet f(opp_cit, opp_cit->index(*cit)); 3080 Vertex_triple vt = make_vertex_triple(f); 3081 make_canonical(vt); 3082 outer_map[vt] = f; 3083 for(int i=0; i<4; i++) 3084 if( i != indv ) 3085 (*cit)->vertex(i)->set_cell(opp_cit); 3086 } 3087 } 3088 3089 /*! \brief removes a vertex from the triangulation. 3090 * 3091 * Removes vertex v from the triangulation. 3092 */ 3093 template < class GT, class TDS > 3094 template < class PointRemover, class Conflict_tester, class CoverManager> 3095 inline void 3096 Periodic_3_triangulation_3<GT,TDS>:: 3097 remove(Vertex_handle v, PointRemover& r, Conflict_tester& t, CoverManager& cover_manager) 3098 { 3099 CGAL_expensive_precondition(is_vertex(v)); 3100 3101 if(!is_1_cover()) 3102 { 3103 std::vector<Vertex_handle> vhrem; 3104 if(number_of_vertices() == 1) { 3105 clear(); 3106 return; 3107 } 3108 Virtual_vertex_map_it vvmit = virtual_vertices.find(v); 3109 if(vvmit != virtual_vertices.end()) v = vvmit->second.first; 3110 CGAL_triangulation_assertion(virtual_vertices_reverse.find(v) 3111 != virtual_vertices_reverse.end()); 3112 vhrem = virtual_vertices_reverse.find(v)->second; 3113 virtual_vertices_reverse.erase(v); 3114 CGAL_triangulation_assertion(vhrem.size()==26); 3115 for(int i=0; i<26; i++) { 3116 periodic_remove(vhrem[i],r, cover_manager); 3117 virtual_vertices.erase(vhrem[i]); 3118 CGAL_triangulation_expensive_assertion(is_valid()); 3119 } 3120 periodic_remove(v,r, cover_manager); 3121 } 3122 else 3123 { 3124 periodic_remove(v, r, cover_manager); 3125 if(!is_1_cover()) 3126 remove(v, r, t, cover_manager); 3127 } 3128 } 3129 3130 /*! \brief removes a vertex from the triangulation. 3131 * 3132 * Removes vertex v from the triangulation. 3133 * It expects a reference to an instance of a PointRemover. 3134 * 3135 * Implementation: 3136 * - Compute the hole, that is, all cells incident to v. Cells outside of 3137 * this hole are not affected by the deletion of v. 3138 * - Triangulate the hole. This is done computing the triangulation 3139 * in Euclidean space for the points on the border of the hole. 3140 * - Sew this triangulation into the hole. 3141 * - Test for all newly added edges, whether they are shorter than the 3142 * edge_length_threshold. If not, convert to 3-cover. 3143 */ 3144 template < class GT, class TDS > 3145 template < class PointRemover, class CoverManager > 3146 inline bool 3147 Periodic_3_triangulation_3<GT,TDS>:: 3148 periodic_remove(Vertex_handle v, PointRemover& remover, CoverManager& cover_manager, 3149 const bool abort_if_cover_change) 3150 { 3151 // Construct the set of vertex triples on the boundary 3152 // with the facet just behind 3153 typedef std::map<Vertex_triple,Facet> Vertex_triple_Facet_map; 3154 typedef PointRemover Point_remover; 3155 typedef typename Point_remover::CellE_handle CellE_handle; 3156 typedef typename Point_remover::VertexE_handle VertexE_handle; 3157 typedef typename Point_remover::FacetE FacetE; 3158 typedef typename Point_remover::VertexE_triple VertexE_triple; 3159 typedef typename Point_remover::Finite_cellsE_iterator 3160 Finite_cellsE_iterator; 3161 typedef typename Point_remover::Vertex_triple_FacetE_map 3162 Vertex_triple_FacetE_map; 3163 3164 // First compute the hole and its boundary vertices. 3165 std::vector<Cell_handle> hole; 3166 hole.reserve(64); 3167 Vertex_triple_Facet_map outer_map; 3168 Vertex_triple_FacetE_map inner_map; 3169 3170 make_hole(v, outer_map, hole); 3171 3172 CGAL_triangulation_assertion(outer_map.size()==hole.size()); 3173 3174 if(!is_1_cover()) { 3175 cover_manager.delete_unsatisfying_elements(hole.begin(), hole.end()); 3176 } 3177 3178 // Build up the map between Vertices on the boundary and offsets 3179 // collect all vertices on the boundary 3180 std::vector<Vertex_handle> vertices; 3181 vertices.reserve(64); 3182 3183 // The set is needed to ensure that each vertex is inserted only once. 3184 std::set<Vertex_handle> tmp_vertices; 3185 // The map connects vertices to offsets in the hole 3186 std::map<Vertex_handle, Offset> vh_off_map; 3187 3188 for(typename std::vector<Cell_handle>::iterator cit = hole.begin(); 3189 cit != hole.end(); ++cit) 3190 { 3191 // Put all incident vertices in tmp_vertices. 3192 for(int j=0; j<4; ++j) { 3193 if((*cit)->vertex(j) != v) { 3194 tmp_vertices.insert((*cit)->vertex(j)); 3195 vh_off_map[(*cit)->vertex(j)] = int_to_off((*cit)->offset(j)) 3196 - int_to_off((*cit)->offset((*cit)->index(v))); 3197 } 3198 } 3199 } 3200 3201 // Now output the vertices. 3202 std::copy(tmp_vertices.begin(), tmp_vertices.end(), std::back_inserter(vertices)); 3203 3204 // create a Delaunay/regular triangulation of the points on the boundary 3205 // in Euclidean space and make a map from the vertices in remover.tmp 3206 // towards the vertices in *this 3207 3208 Unique_hash_map<VertexE_handle,Vertex_handle> vmap; 3209 CellE_handle ch; 3210 remover.tmp.clear(); 3211 3212 for(unsigned int i=0; i < vertices.size(); i++) { 3213 typedef typename Point_remover::Triangulation_R3::Point TRPoint; 3214 CGAL_triangulation_assertion( 3215 get_offset(vertices[i]) + combine_offsets(Offset(), vh_off_map[vertices[i]]) 3216 == combine_offsets(get_offset(vertices[i]),vh_off_map[vertices[i]])); 3217 TRPoint trp = std::make_pair(vertices[i]->point(), 3218 combine_offsets( get_offset(vertices[i]), vh_off_map[vertices[i]]) ); 3219 VertexE_handle vh = remover.tmp.insert(trp, ch); 3220 vmap[vh] = vertices[i]; 3221 CGAL_triangulation_assertion(vmap.is_defined(vh)); 3222 } 3223 CGAL_triangulation_assertion(remover.tmp.number_of_vertices() != 0); 3224 3225 // Construct the set of vertex triples of tmp 3226 // We reorient the vertex triple so that it matches those from outer_map 3227 // Also note that we use the vertices of *this, not of tmp 3228 for(Finite_cellsE_iterator it = remover.tmp.finite_cells_begin(); 3229 it != remover.tmp.finite_cells_end(); 3230 ++it) { 3231 VertexE_triple vt_aux; 3232 for(int i=0; i < 4; i++) { 3233 FacetE f = std::pair<CellE_handle,int>(it,i); 3234 vt_aux = VertexE_triple( 3235 f.first->vertex(vertex_triple_index(f.second,0)), 3236 f.first->vertex(vertex_triple_index(f.second,1)), 3237 f.first->vertex(vertex_triple_index(f.second,2))); 3238 if(vmap.is_defined(vt_aux.first) 3239 && vmap.is_defined(vt_aux.second) 3240 && vmap.is_defined(vt_aux.third) ) { 3241 Vertex_triple vt(vmap[vt_aux.first],vmap[vt_aux.third], 3242 vmap[vt_aux.second]); 3243 make_canonical(vt); 3244 inner_map[vt]= f; 3245 } 3246 } 3247 } 3248 3249 // A structure for storing the new neighboring relations 3250 typedef boost::tuple<Cell_handle, int, Cell_handle> Neighbor_relation; 3251 std::vector<Neighbor_relation> nr_vec; 3252 std::vector<Cell_handle> new_cells; 3253 3254 // Grow inside the hole, by extending the surface 3255 while(! outer_map.empty()) { 3256 typename Vertex_triple_Facet_map::iterator oit = outer_map.begin(); 3257 3258 typename Vertex_triple_Facet_map::value_type o_vt_f_pair = *oit; 3259 Cell_handle o_ch = o_vt_f_pair.second.first; 3260 unsigned int o_i = o_vt_f_pair.second.second; 3261 3262 typename Vertex_triple_FacetE_map::iterator iit = 3263 inner_map.find(o_vt_f_pair.first); 3264 3265 CGAL_triangulation_assertion(iit != inner_map.end()); 3266 typename Vertex_triple_FacetE_map::value_type i_vt_f_pair = *iit; 3267 CellE_handle i_ch = i_vt_f_pair.second.first; 3268 unsigned int i_i = i_vt_f_pair.second.second; 3269 3270 // create a new cell to glue to the outer surface 3271 Cell_handle new_ch = _tds.create_cell(); 3272 new_cells.push_back(new_ch); 3273 new_ch->set_vertices(vmap[i_ch->vertex(0)], vmap[i_ch->vertex(1)], 3274 vmap[i_ch->vertex(2)], vmap[i_ch->vertex(3)]); 3275 set_offsets(new_ch, vh_off_map[vmap[i_ch->vertex(0)]], 3276 vh_off_map[vmap[i_ch->vertex(1)]], 3277 vh_off_map[vmap[i_ch->vertex(2)]], 3278 vh_off_map[vmap[i_ch->vertex(3)]]); 3279 3280 // Update the edge length management 3281 if(cover_manager.update_cover_data_during_management(new_ch, new_cells, 3282 abort_if_cover_change)) 3283 { 3284 CGAL_triangulation_expensive_postcondition(_tds.is_valid()); 3285 return false; // removing would cause / has caused a cover change 3286 } 3287 3288 // The neighboring relation needs to be stored temporarily in 3289 // nr_vec. It cannot be applied directly because then we could not 3290 // easily cancel the removing process if a cell is encountered 3291 // that does not obey the edge-length criterion. 3292 nr_vec.push_back(boost::make_tuple(o_ch,o_i,new_ch)); 3293 nr_vec.push_back(boost::make_tuple(new_ch,i_i,o_ch)); 3294 3295 // for the other faces check, if they can also be glued 3296 for(unsigned int i = 0; i < 4; i++) { 3297 if(i != i_i) { 3298 Facet f = std::pair<Cell_handle,int>(new_ch,i); 3299 Vertex_triple vt = make_vertex_triple(f); 3300 make_canonical(vt); 3301 std::swap(vt.second,vt.third); 3302 typename Vertex_triple_Facet_map::iterator oit2 = outer_map.find(vt); 3303 if(oit2 == outer_map.end()) { 3304 std::swap(vt.second,vt.third); 3305 outer_map[vt]= f; 3306 } else { 3307 // glue the faces 3308 typename Vertex_triple_Facet_map::value_type o_vt_f_pair2 = *oit2; 3309 Cell_handle o_ch2 = o_vt_f_pair2.second.first; 3310 int o_i2 = o_vt_f_pair2.second.second; 3311 nr_vec.push_back(boost::make_tuple(o_ch2,o_i2,new_ch)); 3312 nr_vec.push_back(boost::make_tuple(new_ch,i,o_ch2)); 3313 outer_map.erase(oit2); 3314 } 3315 } 3316 } 3317 outer_map.erase(oit); 3318 } 3319 3320 // finally set the neighboring relations 3321 for(unsigned int i=0; i<nr_vec.size(); i++) { 3322 nr_vec[i].template get<0>()->set_neighbor(nr_vec[i].template get<1>(),nr_vec[i].template get<2>()); 3323 } 3324 3325 // Output the hidden points. 3326 for(typename std::vector<Cell_handle>::iterator 3327 hi = hole.begin(), hend = hole.end(); hi != hend; ++hi) 3328 { 3329 remover.add_hidden_points(*hi); 3330 } 3331 3332 _tds.delete_vertex(v); 3333 _tds.delete_cells(hole.begin(), hole.end()); 3334 CGAL_triangulation_expensive_assertion(is_valid()); 3335 return true; // sucessfully removed the vertex 3336 } 3337 3338 // ############################################################################ 3339 // ############################################################################ 3340 // ############################################################################ 3341 3342 /** \brief deletes each redundant cell and the not anymore needed data 3343 * structures. 3344 * 3345 * This function consists of four iterations over all cells and one 3346 * iteration over all vertices: 3347 * 1. cell iteration: mark all cells that are to delete 3348 * 2. cell iteration: redirect neighbors of remaining cells 3349 * 3. cell iteration: redirect vertices of remaining cells 3350 * 4. cell iteration: delete all cells marked in the 1. iteration 3351 * Vertex iteration: delete all vertices outside the original domain. 3352 */ 3353 template < class GT, class TDS > 3354 inline void 3355 Periodic_3_triangulation_3<GT,TDS>::convert_to_1_sheeted_covering() 3356 { 3357 // ################################################################### 3358 // ### First cell iteration ########################################## 3359 // ################################################################### 3360 { 3361 if(is_1_cover()) 3362 return; 3363 3364 bool to_delete, has_simplifiable_offset; 3365 Virtual_vertex_map_it vvmit; 3366 // First iteration over all cells: Mark the cells that are to be deleted. 3367 // Cells will be deleted if they cannot be translated anymore in the 3368 // direction of one of the axes without yielding negative offsets. 3369 for( Cell_iterator it = all_cells_begin(); 3370 it != all_cells_end(); ++it ) { 3371 to_delete = false; 3372 // for all directions in 3D Space 3373 for( int j=0; j<3; j++ ) { 3374 has_simplifiable_offset = true; 3375 // for all vertices of cell it 3376 for( int i=0; i<4; i++ ) { 3377 vvmit = virtual_vertices.find(it->vertex(i)); 3378 // if it->vertex(i) lies inside the original domain: 3379 if(vvmit == virtual_vertices.end()) { 3380 // the cell cannot be moved any more because if we did, then 3381 // it->vertex(i) will get at least one negative offset. 3382 has_simplifiable_offset = false; 3383 // if it->vertex(i) lies outside the original domain: 3384 } else { 3385 // The cell can certainly be deleted if the offset contains a 2 3386 to_delete = to_delete || (vvmit->second.second[j] == 2); 3387 // The cell can be moved into one direction only if the offset of 3388 // all for vertices is >=1 for this direction. Since we already 3389 // tested for 2 it is sufficient to test here for 1. 3390 has_simplifiable_offset = has_simplifiable_offset 3391 && (vvmit->second.second[j] == 1); 3392 } 3393 } 3394 // if the offset can be simplified, i.e. the cell can be moved, then 3395 // it can be deleted. 3396 if(has_simplifiable_offset) 3397 to_delete = true; 3398 } 3399 // Mark all cells that are to delete. They cannot be deleted yet, 3400 // because neighboring information still needs to be extracted. 3401 if(to_delete) { 3402 it->set_additional_flag(1); 3403 } 3404 } 3405 } 3406 3407 // ################################################################### 3408 // ### Second cell iteration ######################################### 3409 // ################################################################### 3410 { 3411 Vertex_handle vert[4], nbv[4]; 3412 Offset off[4]; 3413 Cell_handle nb, new_neighbor; 3414 std::vector<Triple<Cell_handle, int, Cell_handle> > new_neighbor_relations; 3415 3416 // Second iteration over all cells: redirect neighbors where necessary 3417 for(Cell_iterator it = all_cells_begin(); 3418 it != all_cells_end(); ++it) { 3419 // Skip all cells that are to delete. 3420 if(it->get_additional_flag() == 1) 3421 continue; 3422 3423 // Redirect neighbors: Only neighbors that are marked by the 3424 // additional_flag have to be substituted by one of their periodic 3425 // copies. The unmarked neighbors stay the same. 3426 for( int i = 0; i < 4; i++ ) { 3427 if( it->neighbor(i)->get_additional_flag() != 1 ) 3428 continue; 3429 3430 nb = it->neighbor(i); 3431 3432 for( int j = 0; j < 4; j++ ) { 3433 off[j] = Offset(); 3434 get_vertex( nb, j, vert[j], off[j]); 3435 } 3436 int x,y,z; 3437 x = (std::min) ( (std::min) ( off[0][0], off[1][0] ), 3438 (std::min) ( off[2][0], off[3][0] ) ); 3439 y = (std::min) ( (std::min) ( off[0][1], off[1][1] ), 3440 (std::min) ( off[2][1], off[3][1] ) ); 3441 z = (std::min) ( (std::min) ( off[0][2], off[1][2] ), 3442 (std::min) ( off[2][2], off[3][2] ) ); 3443 3444 // The vector from nb to the "original" periodic copy of nb, that is 3445 // the copy that will not be deleted. 3446 Offset difference_offset(x,y,z); 3447 CGAL_triangulation_assertion( !difference_offset.is_null() ); 3448 3449 // We now have to find the "original" periodic copy of nb from 3450 // its vertices. Therefore, we first have to find the vertices. 3451 for( int j = 0; j < 4; j++ ) { 3452 CGAL_triangulation_assertion( (off[j]-difference_offset)[0] >= 0); 3453 CGAL_triangulation_assertion( (off[j]-difference_offset)[1] >= 0); 3454 CGAL_triangulation_assertion( (off[j]-difference_offset)[2] >= 0); 3455 CGAL_triangulation_assertion( (off[j]-difference_offset)[0] < 3); 3456 CGAL_triangulation_assertion( (off[j]-difference_offset)[1] < 3); 3457 CGAL_triangulation_assertion( (off[j]-difference_offset)[2] < 3); 3458 3459 // find the Vertex_handles of the vertices of the "original" 3460 // periodic copy of nb. If the vertex is inside the original 3461 // domain, there is nothing to do 3462 if( (off[j]-difference_offset).is_null() ) { 3463 nbv[j] = vert[j]; 3464 // If the vertex is outside the original domain, we have to search 3465 // in virtual_vertices in the "wrong" direction. That means, we 3466 // cannot use virtual_vertices.find but have to use 3467 // virtual_vertices_reverse. 3468 } else { 3469 Offset nbo = off[j]-difference_offset; 3470 nbv[j] = virtual_vertices_reverse.find(vert[j]) 3471 ->second[nbo[0]*9+nbo[1]*3+nbo[2]-1]; 3472 } 3473 } 3474 // Find the new neighbor by its 4 vertices 3475 new_neighbor = get_cell( nbv ); 3476 3477 // Store the new neighbor relation. This cannot be applied yet because 3478 // it would disturb the functioning of get_cell( ... ) 3479 new_neighbor_relations.push_back(make_triple(it, i, new_neighbor)); 3480 } 3481 } 3482 // Apply the new neighbor relations now. 3483 for(unsigned int i=0; i<new_neighbor_relations.size(); i++) { 3484 new_neighbor_relations[i].first->set_neighbor( 3485 new_neighbor_relations[i].second, 3486 new_neighbor_relations[i].third); 3487 } 3488 } 3489 3490 // ################################################################### 3491 // ### Third cell iteration ########################################## 3492 // ################################################################### 3493 { 3494 Vertex_handle vert[4]; 3495 Offset off[4]; 3496 // Third iteration over all cells: redirect vertices where necessary 3497 for(Cell_iterator it = all_cells_begin(); 3498 it != all_cells_end(); ++it) { 3499 // Skip all cells that are marked to delete 3500 if(it->get_additional_flag() == 1) 3501 continue; 3502 3503 // Find the corresponding vertices of it in the original domain 3504 // and set them as new vertices of it. 3505 for( int i = 0; i < 4; i++ ) { 3506 off[i] = Offset(); 3507 get_vertex( it, i, vert[i], off[i]); 3508 it->set_vertex( i, vert[i]); 3509 CGAL_triangulation_assertion(vert[i]->point()[0] < domain().xmax()); 3510 CGAL_triangulation_assertion(vert[i]->point()[1] < domain().ymax()); 3511 CGAL_triangulation_assertion(vert[i]->point()[2] < domain().zmax()); 3512 CGAL_triangulation_assertion(vert[i]->point()[0] >= domain().xmin()); 3513 CGAL_triangulation_assertion(vert[i]->point()[1] >= domain().ymin()); 3514 CGAL_triangulation_assertion(vert[i]->point()[2] >= domain().zmin()); 3515 3516 // redirect also the cell pointer of the vertex. 3517 it->vertex(i)->set_cell(it); 3518 } 3519 // Set the offsets. 3520 set_offsets(it, off[0], off[1], off[2], off[3] ); 3521 CGAL_triangulation_assertion( int_to_off(it->offset(0)) == off[0] ); 3522 CGAL_triangulation_assertion( int_to_off(it->offset(1)) == off[1] ); 3523 CGAL_triangulation_assertion( int_to_off(it->offset(2)) == off[2] ); 3524 CGAL_triangulation_assertion( int_to_off(it->offset(3)) == off[3] ); 3525 } 3526 } 3527 3528 // ################################################################### 3529 // ### Fourth cell iteration ######################################### 3530 // ################################################################### 3531 std::vector<Point> hidden_points; 3532 { 3533 // Delete the marked cells. 3534 std::vector<Cell_handle> cells_to_delete; 3535 for( Cell_iterator cit = all_cells_begin(); 3536 cit != all_cells_end(); ++cit ) { 3537 if( cit->get_additional_flag() == 1 ) 3538 { 3539 gather_cell_hidden_points(cit, hidden_points); 3540 cells_to_delete.push_back( cit ); 3541 } 3542 } 3543 _tds.delete_cells(cells_to_delete.begin(), cells_to_delete.end()); 3544 } 3545 3546 // ################################################################### 3547 // ### Vertex iteration ############################################## 3548 // ################################################################### 3549 { 3550 // Delete all the vertices in virtual_vertices, that is all vertices 3551 // outside the original domain. 3552 std::vector<Vertex_handle> vertices_to_delete; 3553 for( Vertex_iterator vit = all_vertices_begin(); 3554 vit != all_vertices_end(); ++vit ) { 3555 if( virtual_vertices.count( vit ) != 0 ) { 3556 CGAL_triangulation_assertion( virtual_vertices.count( vit ) == 1 ); 3557 vertices_to_delete.push_back( vit ); 3558 } 3559 } 3560 _tds.delete_vertices(vertices_to_delete.begin(), vertices_to_delete.end()); 3561 } 3562 _cover = CGAL::make_array(1,1,1); 3563 virtual_vertices.clear(); 3564 virtual_vertices_reverse.clear(); 3565 reinsert_hidden_points_after_converting_to_1_sheeted(hidden_points); 3566 } 3567 3568 template < class GT, class TDS > 3569 inline void 3570 Periodic_3_triangulation_3<GT,TDS>::convert_to_27_sheeted_covering() 3571 { 3572 if(_cover == CGAL::make_array(3,3,3)) 3573 return; 3574 3575 CGAL_triangulation_precondition(is_1_cover()); 3576 3577 // Create 27 copies of each vertex and write virtual_vertices and 3578 // virtual_vertices_reverse 3579 std::list<Vertex_handle> original_vertices; 3580 // try to use std::copy instead of the following loop. 3581 for(Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) 3582 original_vertices.push_back(vit); 3583 for(typename std::list<Vertex_handle>::iterator vit 3584 = original_vertices.begin(); vit != original_vertices.end(); ++vit) { 3585 Vertex_handle v_cp; 3586 std::vector<Vertex_handle> copies; 3587 for(int i=0; i<3; i++) 3588 for(int j=0; j<3; j++) 3589 for(int k=0; k<3; k++) { 3590 if(i==0 && j==0 && k==0) 3591 continue; 3592 v_cp = _tds.create_vertex(*vit); 3593 copies.push_back(v_cp); 3594 virtual_vertices.insert(std::make_pair(v_cp, 3595 std::make_pair(*vit,Offset(i,j,k)))); 3596 } 3597 virtual_vertices_reverse.insert(std::make_pair(*vit,copies)); 3598 } 3599 3600 // Create 27 copies of each cell from the respective copies of the 3601 // vertices and write virtual_cells and virtual_cells_reverse. 3602 typedef std::map<Cell_handle, std::vector<Cell_handle > > 3603 Virtual_cell_reverse_map; 3604 typedef typename Virtual_cell_reverse_map::const_iterator VCRMIT; 3605 3606 Virtual_cell_reverse_map virtual_cells_reverse; 3607 3608 std::list<Cell_handle> original_cells; 3609 for(Cell_iterator cit = cells_begin(); cit != cells_end(); ++cit) 3610 original_cells.push_back(cit); 3611 3612 // Store vertex offsets in a separate data structure 3613 std::list< Offset > off_v; 3614 for(typename std::list<Vertex_handle>::iterator vit 3615 = original_vertices.begin(); vit != original_vertices.end(); ++vit) { 3616 Cell_handle ccc = (*vit)->cell(); 3617 int v_index = ccc->index(*vit); 3618 off_v.push_back(int_to_off(ccc->offset(v_index))); 3619 } 3620 3621 // Store neighboring offsets in a separate data structure 3622 std::list<std::array<Offset,4> > off_nb; 3623 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3624 cit != original_cells.end(); ++cit) { 3625 std::array<Offset,4> off_nb_c; 3626 for(int i=0; i<4; i++) { 3627 Cell_handle ccc = *cit; 3628 Cell_handle nnn = ccc->neighbor(i); 3629 off_nb_c[i] = neighbor_offset(ccc,i,nnn); 3630 } 3631 off_nb.push_back(off_nb_c); 3632 } 3633 3634 // Create copies of cells 3635 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3636 cit != original_cells.end(); ++cit) { 3637 Cell_handle c_cp; 3638 std::vector<Cell_handle> copies; 3639 Virtual_vertex_reverse_map_it vvrmit[4]; 3640 Offset vvoff[4]; 3641 for(int i=0; i<4; i++) { 3642 vvrmit[i] = virtual_vertices_reverse.find((*cit)->vertex(i)); 3643 CGAL_triangulation_assertion( 3644 vvrmit[i] != virtual_vertices_reverse.end()); 3645 vvoff[i] = int_to_off((*cit)->offset(i)); 3646 } 3647 Vertex_handle vvh[4]; 3648 for(int n=0; n<26; n++) { 3649 for(int i=0; i<4; i++) { 3650 // Decomposition of n into an offset (nx,ny,nz): 3651 // nx = (n+1)/9, ny = ((n+1)/3)%3, nz = (n+1)%3 3652 int o_i = ((n+1)/9+vvoff[i].x()+3)%3; 3653 int o_j = ((n+1)/3+vvoff[i].y()+3)%3; 3654 int o_k = ((n+1)+vvoff[i].z()+3)%3; 3655 int n_c = 9*o_i+3*o_j+o_k-1; 3656 CGAL_triangulation_assertion(n_c >= -1); 3657 if(n_c == -1) vvh[i] = (*cit)->vertex(i); 3658 else vvh[i] = vvrmit[i]->second[n_c]; 3659 } 3660 c_cp = _tds.create_cell(vvh[0], vvh[1], vvh[2], vvh[3]); 3661 copies.push_back(c_cp); 3662 } 3663 virtual_cells_reverse.insert(std::make_pair(*cit,copies)); 3664 } 3665 3666 // Set new vertices of boundary cells of the original domain. 3667 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3668 cit != original_cells.end(); ++cit) { 3669 for(int i=0; i<4; i++) { 3670 Virtual_vertex_reverse_map_it vvrmit 3671 = virtual_vertices_reverse.find((*cit)->vertex(i)); 3672 CGAL_triangulation_assertion(vvrmit != virtual_vertices_reverse.end()); 3673 Offset vvoff = int_to_off((*cit)->offset(i)); 3674 if(!vvoff.is_null()) { 3675 int n_c = 9*vvoff.x()+3*vvoff.y()+vvoff.z()-1; 3676 CGAL_triangulation_assertion(n_c >= 0); 3677 CGAL_triangulation_assertion(static_cast<unsigned int>(n_c) 3678 < vvrmit->second.size()); 3679 (*cit)->set_vertex(i,vvrmit->second[n_c]); 3680 } 3681 } 3682 } 3683 3684 // Set neighboring relations of cell copies 3685 typename std::list<std::array<Offset,4> >::iterator oit = off_nb.begin(); 3686 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3687 cit != original_cells.end(); ++cit, ++oit) { 3688 CGAL_triangulation_assertion( oit != off_nb.end() ); 3689 VCRMIT c_cp = virtual_cells_reverse.find(*cit); 3690 CGAL_triangulation_assertion(c_cp != virtual_cells_reverse.end()); 3691 for(int i=0; i<4; i++) { 3692 Cell_handle cit_nb = (*cit)->neighbor(i); 3693 VCRMIT c_cp_nb = virtual_cells_reverse.find(cit_nb); 3694 CGAL_triangulation_assertion(c_cp_nb != virtual_cells_reverse.end()); 3695 Offset nboff = (*oit)[i]; 3696 for(int n=0; n<26; n++) { 3697 int n_nb; 3698 if(nboff.is_null()) n_nb = n; 3699 else { 3700 int o_i = ((n+1)/9-nboff.x()+3)%3; 3701 int o_j = ((n+1)/3-nboff.y()+3)%3; 3702 int o_k = (n+1-nboff.z()+3)%3; 3703 n_nb = 9*o_i+3*o_j+o_k-1; 3704 } 3705 if(n_nb == -1) { 3706 CGAL_triangulation_assertion(cit_nb->has_vertex( 3707 c_cp->second[n]->vertex((i+1)%4)) ); 3708 CGAL_triangulation_assertion(cit_nb->has_vertex( 3709 c_cp->second[n]->vertex((i+2)%4)) ); 3710 CGAL_triangulation_assertion(cit_nb->has_vertex( 3711 c_cp->second[n]->vertex((i+3)%4)) ); 3712 c_cp->second[n]->set_neighbor(i,cit_nb); 3713 } 3714 else { 3715 CGAL_triangulation_assertion(n_nb >= 0); 3716 CGAL_triangulation_assertion(static_cast<unsigned int>(n_nb) 3717 <= c_cp_nb->second.size()); 3718 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3719 ->has_vertex(c_cp->second[n]->vertex((i+1)%4)) ); 3720 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3721 ->has_vertex(c_cp->second[n]->vertex((i+2)%4)) ); 3722 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3723 ->has_vertex(c_cp->second[n]->vertex((i+3)%4)) ); 3724 c_cp->second[n]->set_neighbor(i,c_cp_nb->second[n_nb]); 3725 } 3726 } 3727 } 3728 } 3729 3730 // Set neighboring relations of original cells 3731 oit = off_nb.begin(); 3732 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3733 cit != original_cells.end(); ++cit, ++oit) { 3734 CGAL_triangulation_assertion( oit != off_nb.end() ); 3735 for(int i=0; i<4; i++) { 3736 Offset nboff = (*oit)[i]; 3737 if(!nboff.is_null()) { 3738 Cell_handle cit_nb = (*cit)->neighbor(i); 3739 VCRMIT c_cp_nb = virtual_cells_reverse.find(cit_nb); 3740 CGAL_triangulation_assertion(c_cp_nb != virtual_cells_reverse.end()); 3741 int o_i = (3-nboff.x())%3; 3742 int o_j = (3-nboff.y())%3; 3743 int o_k = (3-nboff.z())%3; 3744 int n_nb = 9*o_i+3*o_j+o_k-1; 3745 CGAL_triangulation_assertion(n_nb >= 0); 3746 CGAL_triangulation_assertion(static_cast<unsigned int>(n_nb) 3747 <= c_cp_nb->second.size()); 3748 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3749 ->has_vertex((*cit)->vertex((i+1)%4)) ); 3750 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3751 ->has_vertex((*cit)->vertex((i+2)%4)) ); 3752 CGAL_triangulation_assertion(c_cp_nb->second[n_nb] 3753 ->has_vertex((*cit)->vertex((i+3)%4)) ); 3754 (*cit)->set_neighbor(i,c_cp_nb->second[n_nb]); 3755 } 3756 } 3757 } 3758 3759 // Set incident cells 3760 for(Cell_iterator cit = cells_begin(); cit != cells_end(); ++cit) { 3761 for(int i=0; i<4; i++) { 3762 cit->vertex(i)->set_cell(cit); 3763 } 3764 } 3765 3766 // Set offsets where necessary 3767 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3768 cit != original_cells.end(); ++cit) { 3769 VCRMIT c_cp = virtual_cells_reverse.find(*cit); 3770 CGAL_triangulation_assertion( c_cp != virtual_cells_reverse.end()); 3771 Offset off[4]; 3772 for(int i=0; i<4; i++) 3773 off[i] = int_to_off((*cit)->offset(i)); 3774 3775 if(off[0].is_null() && off[1].is_null() && 3776 off[2].is_null() && off[3].is_null()) 3777 continue; 3778 3779 for(int n=0; n<26; n++) { 3780 Offset off_cp[4]; 3781 int o_i = (n+1)/9; 3782 int o_j = ((n+1)/3)%3; 3783 int o_k = (n+1)%3; 3784 if(o_i!=2 && o_j!=2 && o_k !=2) 3785 continue; 3786 3787 for(int i=0; i<4; i++) { 3788 off_cp[i] = Offset((o_i==2)?off[i].x():0, 3789 (o_j==2)?off[i].y():0, 3790 (o_k==2)?off[i].z():0); 3791 CGAL_triangulation_assertion(off_cp[i].x() == 0 || off_cp[i].x() == 1); 3792 CGAL_triangulation_assertion(off_cp[i].y() == 0 || off_cp[i].y() == 1); 3793 CGAL_triangulation_assertion(off_cp[i].z() == 0 || off_cp[i].z() == 1); 3794 } 3795 set_offsets(c_cp->second[n],off_cp[0],off_cp[1],off_cp[2],off_cp[3]); 3796 } 3797 } 3798 3799 // Iterate over all original cells and reset offsets. 3800 for(typename std::list<Cell_handle>::iterator cit = original_cells.begin(); 3801 cit != original_cells.end(); ++cit) { 3802 //This statement does not seem to have any effect 3803 set_offsets(*cit, 0,0,0,0); 3804 CGAL_triangulation_assertion((*cit)->offset(0) == 0); 3805 CGAL_triangulation_assertion((*cit)->offset(1) == 0); 3806 CGAL_triangulation_assertion((*cit)->offset(2) == 0); 3807 CGAL_triangulation_assertion((*cit)->offset(3) == 0); 3808 } 3809 3810 _cover = CGAL::make_array(3,3,3); 3811 CGAL_triangulation_expensive_assertion(is_valid()); 3812 3813 update_cover_data_after_converting_to_27_sheeted_covering(); 3814 } 3815 3816 template < class GT, class TDS > 3817 class Periodic_3_triangulation_3<GT,TDS>::Finder 3818 { 3819 const Self* _t; 3820 const Point& _p; 3821 3822 public: 3823 Finder(const Self* t, const Point &p) : _t(t), _p(p) {} 3824 bool operator()(const Vertex_handle v) { 3825 return _t->equal(v->point(), _p); 3826 } 3827 }; 3828 3829 /** Find the cell that consists of the four given vertices 3830 * 3831 * Iterates over all cells and compare the four vertices of each cell 3832 * with the four vertices in vh. 3833 */ 3834 template < class GT, class TDS > 3835 inline typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle 3836 Periodic_3_triangulation_3<GT,TDS>::get_cell(const Vertex_handle* vh) const 3837 { 3838 bool contains_v[4]; 3839 std::vector<Cell_handle> cells; 3840 incident_cells(vh[3],std::back_inserter(cells)); 3841 for( typename std::vector<Cell_handle>::iterator it = cells.begin(); 3842 it != cells.end(); it++ ) { 3843 CGAL_triangulation_assertion( 3844 (*it)->vertex(0) == vh[3] || (*it)->vertex(1) == vh[3] 3845 ||(*it)->vertex(2) == vh[3] || (*it)->vertex(3) == vh[3]); 3846 for( int j=0; j<3; j++ ) { 3847 contains_v[j] = false; 3848 contains_v[j] = ( (*it)->vertex(0) == vh[j] ) 3849 || ( (*it)->vertex(1) == vh[j] ) 3850 || ( (*it)->vertex(2) == vh[j] ) 3851 || ( (*it)->vertex(3) == vh[j] ); 3852 } 3853 if(contains_v[0] && contains_v[1] && contains_v[2]) { 3854 return (*it); 3855 } 3856 } 3857 CGAL_triangulation_assertion(false); 3858 return Cell_handle(); 3859 } 3860 3861 /*! \brief gets the offset of tester.point() such that 3862 * this point is in conflict with c w.r.t tester.get_offset(). 3863 * 3864 * Implementation: Just try all eight possibilities. 3865 */ 3866 template < class GT, class TDS > 3867 template < class Conflict_tester > 3868 inline typename Periodic_3_triangulation_3<GT,TDS>::Offset 3869 Periodic_3_triangulation_3<GT,TDS>::get_location_offset( 3870 const Conflict_tester& tester, Cell_handle c) const 3871 { 3872 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 3873 3874 int cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); 3875 if(cumm_off == 0) { 3876 // default case: 3877 CGAL_triangulation_assertion(tester(c, Offset())); 3878 return Offset(); 3879 } else { 3880 // Main idea seems to just test all possibilities. 3881 for(int i=0; i<8; i++) { 3882 if(((cumm_off | (~i))&7) == 7) { 3883 if(tester(c,int_to_off(i))) { 3884 return int_to_off(i); 3885 } 3886 } 3887 } 3888 } 3889 CGAL_triangulation_assertion(false); 3890 return Offset(); 3891 } 3892 3893 template < class GT, class TDS > 3894 template < class Conflict_tester > 3895 inline typename Periodic_3_triangulation_3<GT,TDS>::Offset 3896 Periodic_3_triangulation_3<GT,TDS>::get_location_offset( 3897 const Conflict_tester& tester, Cell_handle c, bool& found) const 3898 { 3899 CGAL_triangulation_precondition( number_of_vertices() != 0 ); 3900 3901 found = false; 3902 3903 int cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); 3904 if(cumm_off == 0 && tester(c, Offset())) { 3905 // default case: 3906 found = true; 3907 return Offset(); 3908 } else { 3909 // Main idea seems to just test all possibilities. 3910 for(int i=0; i<8; i++) { 3911 if(((cumm_off | (~i))&7) == 7) { 3912 if(tester(c,int_to_off(i))) { 3913 found = true; 3914 return int_to_off(i); 3915 } 3916 } 3917 } 3918 } 3919 3920 return Offset(); 3921 } 3922 3923 /** Get the offset between the origins of the internal offset coordinate 3924 * systems of two neighboring cells with respect from ch to nb. 3925 * 3926 * - Find two corresponding vertices from each cell 3927 * - Return the difference of their offsets. 3928 */ 3929 template < class GT, class TDS > 3930 inline typename Periodic_3_triangulation_3<GT,TDS>::Offset 3931 Periodic_3_triangulation_3<GT,TDS>::neighbor_offset( 3932 Cell_handle ch, int i, Cell_handle nb) const 3933 { 3934 // Redundance in the signature! 3935 CGAL_triangulation_precondition(ch->neighbor(i) == nb); 3936 CGAL_triangulation_precondition(nb->neighbor(nb->index(ch)) == ch); 3937 3938 Vertex_handle vertex_ch; 3939 int index_ch, index_nb; 3940 // ensure that vertex_ch \in nb and vertex_nb \in ch 3941 index_ch = (i==0? 1 : 0); 3942 vertex_ch = ch->vertex(index_ch); 3943 index_nb = nb->index(vertex_ch); 3944 3945 return int_to_off(nb->offset(index_nb)) - int_to_off(ch->offset(index_ch)); 3946 } 3947 3948 /** 3949 * - ch->offset(i) is an bit triple encapsulated in an integer. Each bit 3950 * represents the offset in one direction --> 2-cover! 3951 * - it_to_off(int) decodes this again. 3952 * - Finally the offset vector is multiplied by cover. 3953 * So if we are working in 3-cover we translate it to the neighboring 3954 * 3-cover and not only to the neighboring domain. 3955 */ 3956 template < class GT, class TDS > 3957 inline void Periodic_3_triangulation_3<GT, TDS>::get_vertex( 3958 Cell_handle ch, int i, Vertex_handle& vh, Offset& off) const 3959 { 3960 off = combine_offsets(Offset(),int_to_off(ch->offset(i))); 3961 vh = ch->vertex(i); 3962 3963 if(is_1_cover()) 3964 return; 3965 3966 Vertex_handle vh_i = vh; 3967 get_vertex(vh_i, vh, off); 3968 return; 3969 } 3970 3971 template < class GT, class TDS > 3972 inline void Periodic_3_triangulation_3<GT, TDS>::get_vertex( 3973 Vertex_handle vh_i, Vertex_handle& vh, Offset& off) const 3974 { 3975 Virtual_vertex_map_it it = virtual_vertices.find(vh_i); 3976 3977 if(it == virtual_vertices.end()) { 3978 // if 'vh_i' is not contained in virtual_vertices, then it is in 3979 // the original domain. 3980 vh = vh_i; 3981 CGAL_triangulation_assertion(vh != Vertex_handle()); 3982 } else { 3983 // otherwise it has to be looked up as well as its offset. 3984 vh = it->second.first; 3985 off += it->second.second; 3986 CGAL_triangulation_assertion(vh->point().x() < domain().xmax()); 3987 CGAL_triangulation_assertion(vh->point().y() < domain().ymax()); 3988 CGAL_triangulation_assertion(vh->point().z() < domain().zmax()); 3989 CGAL_triangulation_assertion(vh->point().x() >= domain().xmin()); 3990 CGAL_triangulation_assertion(vh->point().y() >= domain().ymin()); 3991 CGAL_triangulation_assertion(vh->point().z() >= domain().zmin()); 3992 } 3993 } 3994 3995 template < class GT, class TDS > 3996 std::istream& 3997 operator>> (std::istream& is, Periodic_3_triangulation_3<GT,TDS>& tr) 3998 // reads 3999 // the current covering that guarantees the triangulation to be a 4000 // simplicial complex 4001 // the number of vertices 4002 // the non combinatorial information on vertices (points in case of 1-sheeted 4003 // covering, point-offset pairs otherwise) 4004 // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY 4005 // the number of cells 4006 // the cells by the indices of their vertices in the preceding list 4007 // of vertices, plus the non combinatorial information on each cell 4008 // the neighbors of each cell by their index in the preceding list of cells 4009 { 4010 CGAL_triangulation_precondition(is.good()); 4011 4012 typedef Periodic_3_triangulation_3<GT,TDS> Triangulation; 4013 typedef typename Triangulation::size_type size_type; 4014 typedef typename Triangulation::Vertex_handle Vertex_handle; 4015 typedef typename Triangulation::Cell_handle Cell_handle; 4016 typedef typename Triangulation::Offset Offset; 4017 typedef typename Triangulation::Iso_cuboid Iso_cuboid; 4018 4019 tr.clear(); 4020 4021 Iso_cuboid domain(0,0,0,1,1,1); 4022 int cx=0, cy=0, cz=0; 4023 size_type n=0; 4024 4025 if(IO::is_ascii(is)) { 4026 is >> domain; 4027 is >> cx >> cy >> cz; 4028 is >> n; 4029 } 4030 else { 4031 is >> domain; 4032 read(is,cx); 4033 read(is,cy); 4034 read(is,cz); 4035 read(is,n); 4036 } 4037 4038 CGAL_triangulation_assertion((n/(cx*cy*cz))*cx*cy*cz == n); 4039 4040 tr.tds().set_dimension((n==0?-2:3)); 4041 tr._gt.set_domain(domain); 4042 tr._cover = CGAL::make_array(cx,cy,cz); 4043 4044 if( n==0 ) return is; 4045 4046 std::vector< Vertex_handle > V(n); 4047 4048 if(cx==1 && cy==1 && cz==1) { 4049 for(std::size_t i=0; i < n; i++) { 4050 V[i] = tr.tds().create_vertex(); 4051 is >> *V[i]; 4052 } 4053 } else { 4054 Vertex_handle v,w; 4055 std::vector<Vertex_handle> vv; 4056 Offset off; 4057 for(std::size_t i=0; i < n; i++) { 4058 v = tr.tds().create_vertex(); 4059 V[i] = v; 4060 is >> *V[i] >> off; 4061 vv.clear(); 4062 for(int j=1; j<cx*cy*cz; j++) { 4063 i++; 4064 w = tr.tds().create_vertex(); 4065 V[i] = w; 4066 is >> *V[i] >> off; 4067 vv.push_back(w); 4068 tr.virtual_vertices[w]=std::make_pair(v,off); 4069 } 4070 tr.virtual_vertices_reverse[v]=vv; 4071 } 4072 } 4073 4074 std::vector< Cell_handle > C; 4075 std::size_t m; 4076 tr.tds().read_cells(is, V, m, C); 4077 4078 // read offsets 4079 int off[4] = {0,0,0,0}; 4080 for(std::size_t j=0; j < m; j++) { 4081 if(IO::is_ascii(is)) 4082 is >> off[0] >> off[1] >> off[2] >> off[3]; 4083 else { 4084 read(is,off[0]); 4085 read(is,off[1]); 4086 read(is,off[2]); 4087 read(is,off[3]); 4088 } 4089 tr.set_offsets(C[j],off[0],off[1],off[2],off[3]); 4090 } 4091 4092 // read potential other information 4093 for(std::size_t j=0; j < m; j++) 4094 is >> *(C[j]); 4095 4096 CGAL_triangulation_expensive_assertion( tr.is_valid() ); 4097 return is; 4098 } 4099 4100 template < class GT, class TDS > 4101 std::ostream& 4102 operator<< (std::ostream& os,const Periodic_3_triangulation_3<GT,TDS>& tr) 4103 // writes : 4104 // the number of vertices 4105 // the domain as six coordinates: xmin ymin zmin xmax ymax zmax 4106 // the current covering that guarantees the triangulation to be a 4107 // simplicial complex 4108 // the non combinatorial information on vertices (points in case of 1-sheeted 4109 // covering, point-offset pairs otherwise) 4110 // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY 4111 // the number of cells 4112 // the cells by the indices of their vertices in the preceding list 4113 // of vertices, plus the non combinatorial information on each cell 4114 // the neighbors of each cell by their index in the preceding list of cells 4115 { 4116 typedef Periodic_3_triangulation_3<GT,TDS> Triangulation; 4117 typedef typename Triangulation::size_type size_type; 4118 typedef typename Triangulation::Vertex_handle Vertex_handle; 4119 typedef typename Triangulation::Vertex_iterator Vertex_iterator; 4120 typedef typename Triangulation::Cell_handle Cell_handle; 4121 typedef typename Triangulation::Cell_iterator Cell_iterator; 4122 typedef typename Triangulation::Covering_sheets Covering_sheets; 4123 typedef typename Triangulation::Offset Offset; 4124 typedef typename Triangulation::Virtual_vertex_map_it Virtual_vertex_map_it; 4125 typedef typename Triangulation::Iso_cuboid Iso_cuboid; 4126 4127 // outputs dimension, domain and number of vertices 4128 Iso_cuboid domain = tr.domain(); 4129 Covering_sheets cover = tr.number_of_sheets(); 4130 size_type n = tr.number_of_vertices(); 4131 4132 if(IO::is_ascii(os)) 4133 os << domain << std::endl 4134 << cover[0] << " " << cover[1] << " " << cover[2] << std::endl 4135 << n*cover[0]*cover[1]*cover[2] << std::endl; 4136 else { 4137 os << domain; 4138 write(os,cover[0]); 4139 write(os,cover[1]); 4140 write(os,cover[2]); 4141 write(os,n*cover[0]*cover[1]*cover[2]); 4142 } 4143 4144 if(n == 0) 4145 return os; 4146 4147 // write the vertices 4148 Unique_hash_map<Vertex_handle, std::size_t > V; 4149 std::size_t i=0; 4150 if(tr.is_1_cover()) { 4151 for(Vertex_iterator it=tr.vertices_begin(); it!=tr.vertices_end(); ++it) { 4152 V[it] = i++; 4153 os << it->point(); 4154 if(IO::is_ascii(os)) 4155 os << std::endl; 4156 } 4157 } else { 4158 Virtual_vertex_map_it vit, vvit; 4159 std::vector<Vertex_handle> vv; 4160 for(Vertex_iterator it=tr.vertices_begin(); it!=tr.vertices_end(); ++it) { 4161 vit = tr.virtual_vertices.find(it); 4162 if(vit != tr.virtual_vertices.end()) 4163 continue; 4164 4165 V[it]=i++; 4166 if(IO::is_ascii(os)) 4167 os << it->point() << std::endl 4168 << Offset(0,0,0) << std::endl; 4169 else os << it->point() << Offset(0,0,0); 4170 CGAL_triangulation_assertion(tr.virtual_vertices_reverse.find(it) 4171 != tr.virtual_vertices_reverse.end()); 4172 vv = tr.virtual_vertices_reverse.find(it)->second; 4173 CGAL_triangulation_assertion(vv.size() == 26); 4174 for(std::size_t j=0; j<vv.size(); j++) { 4175 vvit = tr.virtual_vertices.find(vv[j]); 4176 CGAL_triangulation_assertion(vvit != tr.virtual_vertices.end()); 4177 V[vv[j]] = i++; 4178 if(IO::is_ascii(os)) 4179 os << vv[j]->point() << std::endl 4180 << vvit->second.second << std::endl; 4181 else os << vv[j]->point() << vvit->second.second; 4182 } 4183 } 4184 } 4185 CGAL_triangulation_postcondition(i == tr.number_of_sheets()[0] * 4186 tr.number_of_sheets()[1] * 4187 tr.number_of_sheets()[2] * n); 4188 4189 // asks the tds for the combinatorial information 4190 tr.tds().print_cells(os, V); 4191 4192 // write offsets 4193 //for(unsigned int i=0; i<tr.number_of_cells(); i++) { 4194 for(Cell_iterator it=tr.cells_begin(); it!=tr.cells_end(); ++it) { 4195 //Cell_handle ch = std::find(tr.cells_begin(), tr.cells_end(), i); 4196 Cell_handle ch(it); 4197 for(int j=0; j<4; j++) { 4198 if(IO::is_ascii(os)) { 4199 os << ch->offset(j); 4200 if( j==3 ) 4201 os << std::endl; 4202 else 4203 os << ' '; 4204 } 4205 else write(os,ch->offset(j)); 4206 } 4207 } 4208 4209 // write the non combinatorial information on the cells 4210 // using the << operator of Cell 4211 // works because the iterator of the tds traverses the cells in the 4212 // same order as the iterator of the triangulation 4213 if(tr.number_of_vertices() != 0) { 4214 for(Cell_iterator it=tr.cells_begin(); it != tr.cells_end(); ++it) { 4215 os << *it; // other information 4216 if(IO::is_ascii(os)) 4217 os << std::endl; 4218 } 4219 } 4220 return os; 4221 } 4222 4223 namespace internal { 4224 4225 /// Internal function used by operator==. 4226 // This function tests and registers the 4 neighbors of c1/c2, 4227 // and performs a bfs traversal 4228 // Returns false if an inequality has been found. 4229 //TODO: introduce offsets 4230 template <class GT, class TDS1, class TDS2> 4231 bool 4232 test_next(const Periodic_3_triangulation_3<GT, TDS1>& t1, 4233 const Periodic_3_triangulation_3<GT, TDS2>& t2, 4234 typename Periodic_3_triangulation_3<GT, TDS1>::Cell_handle c1, 4235 typename Periodic_3_triangulation_3<GT, TDS2>::Cell_handle c2, 4236 std::map<typename Periodic_3_triangulation_3<GT, TDS1>::Cell_handle, 4237 typename Periodic_3_triangulation_3<GT, TDS2>::Cell_handle>& Cmap, 4238 std::map<typename Periodic_3_triangulation_3<GT, TDS1>::Vertex_handle, 4239 typename Periodic_3_triangulation_3<GT, TDS2>::Vertex_handle>& Vmap) 4240 { 4241 typedef Periodic_3_triangulation_3<GT, TDS1> Tr1; 4242 typedef Periodic_3_triangulation_3<GT, TDS2> Tr2; 4243 typedef typename Tr1::Vertex_handle Vertex_handle1; 4244 typedef typename Tr1::Cell_handle Cell_handle1; 4245 typedef typename Tr2::Vertex_handle Vertex_handle2; 4246 typedef typename Tr2::Cell_handle Cell_handle2; 4247 typedef typename std::map<Cell_handle1, Cell_handle2>::const_iterator Cit; 4248 typedef typename std::map<Vertex_handle1, Vertex_handle2>::const_iterator Vit; 4249 4250 std::vector<std::pair<Cell_handle1, Cell_handle2> > queue; 4251 queue.push_back(std::make_pair(c1,c2)); 4252 4253 while(! queue.empty()) { 4254 boost::tie(c1,c2) = queue.back(); 4255 queue.pop_back(); 4256 4257 // Precondition: c1, c2 have been registered as well as their 4 vertices. 4258 CGAL_triangulation_precondition(t1.number_of_vertices() != 0); 4259 CGAL_triangulation_precondition(Cmap[c1] == c2); 4260 CGAL_triangulation_precondition(Vmap.find(c1->vertex(0)) != Vmap.end()); 4261 CGAL_triangulation_precondition(Vmap.find(c1->vertex(1)) != Vmap.end()); 4262 CGAL_triangulation_precondition(Vmap.find(c1->vertex(2)) != Vmap.end()); 4263 CGAL_triangulation_precondition(Vmap.find(c1->vertex(3)) != Vmap.end()); 4264 4265 for(int i=0; i <= 3; ++i) { 4266 Cell_handle1 n1 = c1->neighbor(i); 4267 Cit cit = Cmap.find(n1); 4268 Vertex_handle1 v1 = c1->vertex(i); 4269 Vertex_handle2 v2 = Vmap[v1]; 4270 Cell_handle2 n2 = c2->neighbor(c2->index(v2)); 4271 if(cit != Cmap.end()) { 4272 // n1 was already registered. 4273 if(cit->second != n2) 4274 return false; 4275 continue; 4276 } 4277 // n1 has not yet been registered. 4278 // We check that the new vertices match geometrically. 4279 // And we register them. 4280 Vertex_handle1 vn1 = n1->vertex(n1->index(c1)); 4281 Vertex_handle2 vn2 = n2->vertex(n2->index(c2)); 4282 Vit vit = Vmap.find(vn1); 4283 if(vit != Vmap.end()) { 4284 // vn1 already registered 4285 if(vit->second != vn2) 4286 return false; 4287 } 4288 else { 4289 if(t1.geom_traits().compare_xyz_3_object()(t1.construct_point(vn1->point()), 4290 t2.construct_point(vn2->point())) != 0) 4291 return false; 4292 4293 // We register vn1/vn2. 4294 Vmap.insert(std::make_pair(vn1, vn2)); 4295 } 4296 4297 // We register n1/n2. 4298 Cmap.insert(std::make_pair(n1, n2)); 4299 queue.push_back(std::make_pair(n1, n2)); 4300 } 4301 } 4302 return true; 4303 } 4304 4305 } // namespace internal 4306 4307 4308 template < class GT, class TDS1, class TDS2 > 4309 bool 4310 operator==(const Periodic_3_triangulation_3<GT,TDS1>& t1, 4311 const Periodic_3_triangulation_3<GT,TDS2>& t2) 4312 { 4313 typedef typename Periodic_3_triangulation_3<GT,TDS1>::Vertex_handle 4314 Vertex_handle1; 4315 typedef typename Periodic_3_triangulation_3<GT,TDS1>::Cell_handle 4316 Cell_handle1; 4317 typedef typename Periodic_3_triangulation_3<GT,TDS2>::Vertex_handle 4318 Vertex_handle2; 4319 typedef typename Periodic_3_triangulation_3<GT,TDS2>::Vertex_handle 4320 Vertex_iterator2; 4321 typedef typename Periodic_3_triangulation_3<GT,TDS2>::Cell_handle 4322 Cell_handle2; 4323 4324 typedef typename Periodic_3_triangulation_3<GT,TDS1>::Point Point; 4325 typedef typename Periodic_3_triangulation_3<GT,TDS1>::Offset Offset; 4326 4327 // typedef typename Periodic_3_triangulation_3<GT,TDS1> 4328 // ::Geometric_traits::Compare_xyz_3 Compare_xyz_3; 4329 // Compare_xyz_3 cmp1 = t1.geom_traits().compare_xyz_3_object(); 4330 // Compare_xyz_3 cmp2 = t2.geom_traits().compare_xyz_3_object(); 4331 4332 // Some quick checks. 4333 if( t1.domain() != t2.domain() 4334 || t1.number_of_sheets() != t2.number_of_sheets()) 4335 return false; 4336 4337 if( t1.number_of_vertices() != t2.number_of_vertices() 4338 || t1.number_of_cells() != t2.number_of_cells()) 4339 return false; 4340 4341 // Special case for empty triangulations 4342 if(t1.number_of_vertices() == 0) 4343 return true; 4344 4345 // We will store the mapping between the 2 triangulations vertices and 4346 // cells in 2 maps. 4347 std::map<Vertex_handle1, Vertex_handle2> Vmap; 4348 std::map<Cell_handle1, Cell_handle2> Cmap; 4349 4350 // find a common point 4351 Vertex_handle1 v1 = static_cast<Vertex_handle1>(t1.vertices_begin()); 4352 Vertex_handle2 iv2; 4353 for(Vertex_iterator2 vit2 = t2.vertices_begin(); 4354 vit2 != t2.vertices_end(); ++vit2) { 4355 if(!t1.equal(vit2->point(), v1->point(), 4356 t2.get_offset(vit2), t1.get_offset(v1))) 4357 continue; 4358 iv2 = static_cast<Vertex_handle2>(vit2); 4359 break; 4360 } 4361 if(iv2 == Vertex_handle2()) 4362 return false; 4363 Vmap.insert(std::make_pair(v1, iv2)); 4364 4365 // We pick one cell of t1, and try to match it against the 4366 // cells of t2. 4367 Cell_handle1 c = v1->cell(); 4368 Vertex_handle1 v2 = c->vertex((c->index(v1)+1)%4); 4369 Vertex_handle1 v3 = c->vertex((c->index(v1)+2)%4); 4370 Vertex_handle1 v4 = c->vertex((c->index(v1)+3)%4); 4371 Point p2 = v2->point(); 4372 Point p3 = v3->point(); 4373 Point p4 = v4->point(); 4374 Offset o2 = t1.get_offset(v2); 4375 Offset o3 = t1.get_offset(v3); 4376 Offset o4 = t1.get_offset(v4); 4377 4378 std::vector<Cell_handle2> ics; 4379 t2.incident_cells(iv2, std::back_inserter(ics)); 4380 for(typename std::vector<Cell_handle2>::const_iterator cit = ics.begin(); 4381 cit != ics.end(); ++cit) { 4382 int inf = (*cit)->index(iv2); 4383 4384 if(t1.equal(p2, (*cit)->vertex((inf+1)%4)->point(), 4385 o2, t2.get_offset((*cit)->vertex((inf+1)%4)))) 4386 Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+1)%4))); 4387 else if(t1.equal(p2, (*cit)->vertex((inf+2)%4)->point(), 4388 o2, t2.get_offset((*cit)->vertex((inf+2)%4)))) 4389 Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+2)%4))); 4390 else if(t1.equal(p2, (*cit)->vertex((inf+3)%4)->point(), 4391 o2, t2.get_offset((*cit)->vertex((inf+3)%4)))) 4392 Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+3)%4))); 4393 else 4394 continue; // None matched v2. 4395 4396 if(t1.equal(p3, (*cit)->vertex((inf+1)%4)->point(), 4397 o3, t2.get_offset((*cit)->vertex((inf+1)%4)))) 4398 Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+1)%4))); 4399 else if(t1.equal(p3, (*cit)->vertex((inf+2)%4)->point(), 4400 o3, t2.get_offset((*cit)->vertex((inf+2)%4)))) 4401 Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+2)%4))); 4402 else if(t1.equal(p3, (*cit)->vertex((inf+3)%4)->point(), 4403 o3, t2.get_offset((*cit)->vertex((inf+3)%4)))) 4404 Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+3)%4))); 4405 else 4406 continue; // None matched v3. 4407 4408 if(t1.equal(p4, (*cit)->vertex((inf+1)%4)->point(), 4409 o4, t2.get_offset((*cit)->vertex((inf+1)%4)))) 4410 Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+1)%4))); 4411 else if(t1.equal(p4, (*cit)->vertex((inf+2)%4)->point(), 4412 o4, t2.get_offset((*cit)->vertex((inf+2)%4)))) 4413 Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+2)%4))); 4414 else if(t1.equal(p4, (*cit)->vertex((inf+3)%4)->point(), 4415 o4, t2.get_offset((*cit)->vertex((inf+3)%4)))) 4416 Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+3)%4))); 4417 else 4418 continue; // None matched v4. 4419 4420 // Found it ! 4421 Cmap.insert(std::make_pair(c, *cit)); 4422 break; 4423 } 4424 4425 if(Cmap.size() == 0) 4426 return false; 4427 4428 // We now have one cell, we need to compare in a bfs graph traversal 4429 return internal::test_next(t1, t2, Cmap.begin()->first, Cmap.begin()->second, Cmap, Vmap); 4430 } 4431 4432 template < class GT, class TDS1, class TDS2 > 4433 inline 4434 bool 4435 operator!=(const Periodic_3_triangulation_3<GT,TDS1>& t1, 4436 const Periodic_3_triangulation_3<GT,TDS2>& t2) 4437 { 4438 return ! (t1 == t2); 4439 } 4440 4441 } // namespace CGAL 4442 4443 #endif // CGAL_PERIODIC_3_TRIANGULATION_3_H 4444