// Copyright (c) 1999-2003,2006-2009 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // You can redistribute it and/or modify it under the terms of the GNU // General Public License as published by the Free Software Foundation, // either version 3 of the License, or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Monique Teillaud // Sylvain Pion // Nico Kruithof // Manuel Caroli #ifndef CGAL_PERIODIC_3_TRIANGULATION_3_H #define CGAL_PERIODIC_3_TRIANGULATION_3_H #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifndef CGAL_NO_STRUCTURAL_FILTERING #include #include #endif // no CGAL_NO_STRUCTURAL_FILTERING namespace CGAL { template < class GT, class TDS > class Periodic_3_triangulation_3; template < class GT, class TDS > std::istream& operator>> (std::istream& is, Periodic_3_triangulation_3 &tr); template < class GT, class TDS > std::ostream& operator<< (std::ostream& os, const Periodic_3_triangulation_3 &tr); #ifndef CGAL_NO_STRUCTURAL_FILTERING namespace internal { // structural filtering is performed only for EPIC struct Periodic_structural_filtering_3_tag {}; struct No_periodic_structural_filtering_3_tag {}; template struct Periodic_structural_filtering_selector_3 { #ifdef FORCE_STRUCTURAL_FILTERING typedef Periodic_structural_filtering_3_tag Tag; #else typedef No_periodic_structural_filtering_3_tag Tag; #endif }; template <> struct Periodic_structural_filtering_selector_3 { typedef Periodic_structural_filtering_3_tag Tag; }; } #endif // no CGAL_NO_STRUCTURAL_FILTERING /**\class Periodic_3_triangulation_3 * * \brief Implements functionality for computing in periodic space. * * There are several things that are special to computing in $\mathbb{T}^3$ * such as * - periodicity --> offsets * - no infinite vertex * - no degenerate dimensions * All functions that are affected can be found in this class. In case it is * necessary to provide different implementations for Delaunay and regular * triangulation, we work with visitors. */ template < class GT, class TDS = Triangulation_data_structure_3 < Triangulation_vertex_base_3< GT, Periodic_3_triangulation_ds_vertex_base_3<> >, Triangulation_cell_base_3< GT, Periodic_3_triangulation_ds_cell_base_3<> > > > class Periodic_3_triangulation_3 : public Triangulation_utils_3 { friend std::istream& operator>> <> (std::istream& is, Periodic_3_triangulation_3 &tr); typedef Periodic_3_triangulation_3 Self; public: typedef GT Geometric_traits; typedef TDS Triangulation_data_structure; typedef typename GT::Periodic_3_offset_3 Offset; typedef typename GT::Iso_cuboid_3 Iso_cuboid; typedef array Covering_sheets; typedef typename GT::Point_3 Point; typedef typename GT::Segment_3 Segment; typedef typename GT::Triangle_3 Triangle; typedef typename GT::Tetrahedron_3 Tetrahedron; typedef std::pair Periodic_point; typedef array< std::pair, 2> Periodic_segment; typedef array< std::pair, 3> Periodic_triangle; typedef array< std::pair, 4> Periodic_tetrahedron; typedef typename TDS::Vertex Vertex; typedef typename TDS::Cell Cell; typedef typename TDS::Facet Facet; typedef typename TDS::Edge Edge; typedef typename TDS::Vertex_handle Vertex_handle; typedef typename TDS::Cell_handle Cell_handle; typedef typename TDS::size_type size_type; typedef typename TDS::difference_type difference_type; typedef typename TDS::Cell_iterator Cell_iterator; typedef typename TDS::Facet_iterator Facet_iterator; typedef typename TDS::Edge_iterator Edge_iterator; typedef typename TDS::Vertex_iterator Vertex_iterator; typedef typename TDS::Cell_circulator Cell_circulator; typedef typename TDS::Facet_circulator Facet_circulator; typedef Cell_iterator All_cells_iterator; typedef Facet_iterator All_facets_iterator; typedef Edge_iterator All_edges_iterator; typedef Vertex_iterator All_vertices_iterator; typedef Periodic_3_triangulation_unique_vertex_iterator_3 Unique_vertex_iterator; private: typedef typename GT::FT FT; typedef std::pair< Vertex_handle, Offset > Virtual_vertex; typedef std::map Virtual_vertex_map; typedef typename Virtual_vertex_map::const_iterator Virtual_vertex_map_it; typedef std::map > Virtual_vertex_reverse_map; typedef typename Virtual_vertex_reverse_map::const_iterator Virtual_vertex_reverse_map_it; typedef Triple< Vertex_handle, Vertex_handle, Vertex_handle > Vertex_triple; public: typedef Periodic_3_triangulation_tetrahedron_iterator_3 Periodic_tetrahedron_iterator; typedef Periodic_3_triangulation_triangle_iterator_3 Periodic_triangle_iterator; typedef Periodic_3_triangulation_segment_iterator_3 Periodic_segment_iterator; typedef Periodic_3_triangulation_point_iterator_3 Periodic_point_iterator; typedef Point value_type; typedef const value_type& const_reference; typedef Tag_false Weighted_tag; public: enum Iterator_type { STORED=0, UNIQUE, //1 STORED_COVER_DOMAIN, //2 UNIQUE_COVER_DOMAIN };//3 enum Locate_type { VERTEX=0, EDGE, //1 FACET, //2 CELL, //3 EMPTY , //4 OUTSIDE_CONVEX_HULL, // unused, for compatibility with Alpha_shape_3 OUTSIDE_AFFINE_HULL }; // unused, for compatibility with Alpha_shape_3 private: Geometric_traits _gt; Triangulation_data_structure _tds; Iso_cuboid _domain; /// This threshold should be chosen such that if all edges are shorter, /// we can be sure that there are no self-edges anymore. FT edge_length_threshold; /// This adjacency list stores all edges that are longer than /// edge_length_threshold. std::map< Vertex_handle, std::list > too_long_edges; unsigned int too_long_edge_counter; /// map of offsets for periodic copies of vertices Virtual_vertex_map virtual_vertices; Virtual_vertex_reverse_map virtual_vertices_reverse; protected: /// v_offsets temporarily stores all the vertices on the border of a /// conflict region. mutable std::vector v_offsets; private: /// Determines if we currently compute in 3-cover or 1-cover. Covering_sheets _cover; public: /** @name Creation */ //@{ Periodic_3_triangulation_3( const Iso_cuboid & domain = Iso_cuboid(0,0,0,1,1,1), const Geometric_traits & gt = Geometric_traits()) : _gt(gt), _tds(), _domain(domain), too_long_edge_counter(0) { _gt.set_domain(_domain); typedef typename internal::Exact_field_selector::Type EFT; typedef NT_converter NTC; CGAL_USE_TYPE(NTC); CGAL_triangulation_precondition_code( NTC ntc; ) CGAL_triangulation_precondition(ntc(_domain.xmax())-ntc(_domain.xmin()) == ntc(_domain.ymax())-ntc(_domain.ymin())); CGAL_triangulation_precondition(ntc(_domain.ymax())-ntc(_domain.ymin()) == ntc(_domain.zmax())-ntc(_domain.zmin())); CGAL_triangulation_precondition(ntc(_domain.zmax())-ntc(_domain.zmin()) == ntc(_domain.xmax())-ntc(_domain.xmin())); _cover = make_array(3,3,3); init_tds(); edge_length_threshold = FT(0.166) * (_domain.xmax()-_domain.xmin()) * (_domain.xmax()-_domain.xmin()); } private: // Copy constructor helpers class Finder; void copy_multiple_covering(const Periodic_3_triangulation_3 & tr); public: // Copy constructor duplicates vertices and cells Periodic_3_triangulation_3(const Periodic_3_triangulation_3 & tr) : _gt(tr.geom_traits()), _domain(tr._domain), edge_length_threshold(tr.edge_length_threshold), _cover(tr._cover) { if (is_1_cover()) { _tds = tr.tds(); } else { copy_multiple_covering(tr); } CGAL_triangulation_expensive_postcondition(*this == tr); } /** @name Assignment */ //@{ Periodic_3_triangulation_3 & operator=(Periodic_3_triangulation_3 tr) { swap(tr); return *this; } void swap(Periodic_3_triangulation_3 &tr) { std::swap(tr._gt, _gt); _tds.swap(tr._tds); std::swap(_domain,tr._domain); std::swap(edge_length_threshold,tr.edge_length_threshold); std::swap(too_long_edges,tr.too_long_edges); std::swap(too_long_edge_counter,tr.too_long_edge_counter); std::swap(virtual_vertices,tr.virtual_vertices); std::swap(virtual_vertices_reverse,tr.virtual_vertices_reverse); std::swap(_cover, tr._cover); } /// Clears the triangulation and initializes it again. void clear() { _tds.clear(); init_tds(); too_long_edges.clear(); too_long_edge_counter = 0; virtual_vertices.clear(); virtual_vertices_reverse.clear(); _cover = make_array(3,3,3); v_offsets.clear(); } //@} private: /// Initializes the triangulation data structure void init_tds() { _tds.set_dimension(-2); v_offsets.reserve(48); } public: /** @name Access functions */ //@{ const Geometric_traits& geom_traits() const { return _gt; } const TDS & tds() const { return _tds; } TDS & tds() { return _tds; } const Iso_cuboid & domain() const { return _domain; } // TODO: Documentation and tests void set_domain(const Iso_cuboid & domain) { clear(); _domain = domain; _gt.set_domain(domain); edge_length_threshold = FT(0.166) * (_domain.xmax()-_domain.xmin()) * (_domain.xmax()-_domain.xmin()); } const Covering_sheets & number_of_sheets() const { return _cover; } const std::pair original_vertex( const Vertex_handle v) const { return (virtual_vertices.find(v) == virtual_vertices.end()) ? std::make_pair(v,Offset()) : virtual_vertices.find(v)->second; } const std::vector& periodic_copies( const Vertex_handle v) const { CGAL_triangulation_precondition(number_of_sheets() != make_array(1,1,1) ); CGAL_triangulation_precondition( virtual_vertices.find(v) == virtual_vertices.end()); CGAL_triangulation_assertion( virtual_vertices_reverse.find(v) != virtual_vertices_reverse.end()); return virtual_vertices_reverse.find(v)->second; } bool is_extensible_triangulation_in_1_sheet_h1() const; bool is_extensible_triangulation_in_1_sheet_h2() const; bool is_triangulation_in_1_sheet() const; void convert_to_1_sheeted_covering(); void convert_to_27_sheeted_covering(); size_type number_of_cells() const { if (is_1_cover()) return _tds.number_of_cells(); else return _tds.number_of_cells()/27; } size_type number_of_facets() const { if (is_1_cover()) return _tds.number_of_facets(); else return _tds.number_of_facets()/27; } size_type number_of_edges() const { if (is_1_cover()) return _tds.number_of_edges(); else return _tds.number_of_edges()/27; } size_type number_of_vertices() const { if (is_1_cover()) return _tds.number_of_vertices(); else return _tds.number_of_vertices()/27; } size_type number_of_stored_cells() const { return _tds.number_of_cells(); } size_type number_of_stored_facets() const { return _tds.number_of_facets(); } size_type number_of_stored_edges() const { return _tds.number_of_edges(); } size_type number_of_stored_vertices() const { return _tds.number_of_vertices(); } protected: bool is_1_cover() const { bool flag; flag = ((_cover[0] == 1) && (_cover[1] == 1) && (_cover[2] == 1)); return flag; } public: bool is_virtual(Vertex_handle v) { if (is_1_cover()) return false; return (virtual_vertices.find(v) != virtual_vertices.end()); } public: // Offset converters int off_to_int(const Offset & off) const { CGAL_triangulation_assertion( off.x()==0 || off.x() ==1 ); CGAL_triangulation_assertion( off.y()==0 || off.y() ==1 ); CGAL_triangulation_assertion( off.z()==0 || off.z() ==1 ); int i = ((off.x()&1)<<2) + ((off.y()&1)<<1) + ((off.z()&1)); return i; } Offset int_to_off(int i) const { return Offset((i>>2)&1,(i>>1)&1,i&1); } void set_offsets(Cell_handle c, int o0,int o1,int o2,int o3) { int off0[3] = {(o0>>2)&1,(o0>>1)&1,(o0&1)}; int off1[3] = {(o1>>2)&1,(o1>>1)&1,(o1&1)}; int off2[3] = {(o2>>2)&1,(o2>>1)&1,(o2&1)}; int off3[3] = {(o3>>2)&1,(o3>>1)&1,(o3&1)}; for (int i=0; i<3; i++) { int min_off = (std::min)((std::min)(off0[i],off1[i]), (std::min)(off2[i],off3[i])); if (min_off != 0) { off0[i] -= min_off; off1[i] -= min_off; off2[i] -= min_off; off3[i] -= min_off; } } o0 = ((off0[0]&1)<<2)+((off0[1]&1)<<1)+(off0[2]&1); o1 = ((off1[0]&1)<<2)+((off1[1]&1)<<1)+(off1[2]&1); o2 = ((off2[0]&1)<<2)+((off2[1]&1)<<1)+(off2[2]&1); o3 = ((off3[0]&1)<<2)+((off3[1]&1)<<1)+(off3[2]&1); c->set_offsets(o0,o1,o2,o3); } template void set_offsets(Cell_handle c, Offset o0,Offset o1,Offset o2,Offset o3) { int off0[3] = {o0.x(),o0.y(),o0.z()}; int off1[3] = {o1.x(),o1.y(),o1.z()}; int off2[3] = {o2.x(),o2.y(),o2.z()}; int off3[3] = {o3.x(),o3.y(),o3.z()}; for (int i=0; i<3; i++) { int min_off = (std::min)((std::min)(off0[i],off1[i]), (std::min)(off2[i],off3[i])); if (min_off != 0) { off0[i] -= min_off; off1[i] -= min_off; off2[i] -= min_off; off3[i] -= min_off; } } CGAL_triangulation_assertion((std::min)((std::min)(off0[0],off1[0]), (std::min)(off2[0],off3[0])) == 0); CGAL_triangulation_assertion((std::min)((std::min)(off0[1],off1[1]), (std::min)(off2[1],off3[1])) == 0); CGAL_triangulation_assertion((std::min)((std::min)(off0[2],off1[2]), (std::min)(off2[2],off3[2])) == 0); CGAL_triangulation_assertion((0 <= off0[0]) && (off0[0] < 2)); CGAL_triangulation_assertion((0 <= off1[0]) && (off1[0] < 2)); CGAL_triangulation_assertion((0 <= off2[0]) && (off2[0] < 2)); CGAL_triangulation_assertion((0 <= off3[0]) && (off3[0] < 2)); CGAL_triangulation_assertion((0 <= off0[1]) && (off0[1] < 2)); CGAL_triangulation_assertion((0 <= off1[1]) && (off1[1] < 2)); CGAL_triangulation_assertion((0 <= off2[1]) && (off2[1] < 2)); CGAL_triangulation_assertion((0 <= off3[1]) && (off3[1] < 2)); CGAL_triangulation_assertion((0 <= off0[2]) && (off0[2] < 2)); CGAL_triangulation_assertion((0 <= off1[2]) && (off1[2] < 2)); CGAL_triangulation_assertion((0 <= off2[2]) && (off2[2] < 2)); CGAL_triangulation_assertion((0 <= off3[2]) && (off3[2] < 2)); int o0i = ((off0[0]&1)<<2)+((off0[1]&1)<<1)+(off0[2]&1); int o1i = ((off1[0]&1)<<2)+((off1[1]&1)<<1)+(off1[2]&1); int o2i = ((off2[0]&1)<<2)+((off2[1]&1)<<1)+(off2[2]&1); int o3i = ((off3[0]&1)<<2)+((off3[1]&1)<<1)+(off3[2]&1); c->set_offsets(o0i,o1i,o2i,o3i); } public: /** @name Wrapping the traits */ //@{ Comparison_result compare_xyz(const Point &p1, const Point &p2) const { return geom_traits().compare_xyz_3_object()(p1,p2); } Comparison_result compare_xyz(const Point &p1, const Point&p2, const Offset &o1, const Offset &o2) const { return geom_traits().compare_xyz_3_object()(p1,p2,o1,o2); } Orientation orientation( const Point &p1, const Point &p2, const Point &p3, const Point &p4) const { return geom_traits().orientation_3_object()(p1,p2,p3,p4); } Orientation orientation( const Point &p1, const Point &p2, const Point &p3, const Point &p4, const Offset &o1, const Offset &o2, const Offset &o3, const Offset &o4) const { return geom_traits().orientation_3_object()(p1,p2,p3,p4,o1,o2,o3,o4); } bool equal(const Point &p1, const Point &p2) const { return compare_xyz(p1,p2) == EQUAL; } bool equal(const Point &p1, const Point &p2, const Offset &o1, const Offset &o2) const { return compare_xyz(p1,p2,o1,o2) == EQUAL; } bool coplanar( const Point &p1, const Point &p2, const Point &p3, const Point &p4) const { return orientation(p1,p2,p3,p4) == COPLANAR; } bool coplanar( const Point &p1, const Point &p2, const Point &p3, const Point &p4, const Offset &o1, const Offset &o2, const Offset &o3, const Offset &o4) const { return orientation(p1,p2,p3,p4,o1,o2,o3,o4) == COPLANAR; } Periodic_tetrahedron construct_periodic_3_tetrahedron( const Point &p1, const Point &p2, const Point &p3, const Point &p4) const { return make_array(std::make_pair(p1,Offset()), std::make_pair(p2,Offset()), std::make_pair(p3,Offset()), std::make_pair(p4,Offset())); } Periodic_tetrahedron construct_periodic_3_tetrahedron( const Point &p1, const Point &p2, const Point &p3, const Point &p4, const Offset &o1, const Offset &o2, const Offset &o3, const Offset &o4) const { return make_array(std::make_pair(p1,o1), std::make_pair(p2,o2), std::make_pair(p3,o3), std::make_pair(p4,o4)); } Periodic_triangle construct_periodic_3_triangle( const Point &p1, const Point &p2, const Point &p3) const { return make_array(std::make_pair(p1,Offset()), std::make_pair(p2,Offset()), std::make_pair(p3,Offset())); } Periodic_triangle construct_periodic_3_triangle( const Point &p1, const Point &p2, const Point &p3, const Offset &o1, const Offset &o2, const Offset &o3) const { return make_array(std::make_pair(p1,o1), std::make_pair(p2,o2), std::make_pair(p3,o3)); } Periodic_segment construct_periodic_3_segment( const Point &p1, const Point &p2) const { return make_array(std::make_pair(p1,Offset()), std::make_pair(p2,Offset())); } Periodic_segment construct_periodic_3_segment( const Point &p1, const Point &p2, const Offset &o1, const Offset &o2) const { return make_array(std::make_pair(p1,o1), std::make_pair(p2,o2)); } Tetrahedron construct_tetrahedron( const Point &p1, const Point &p2, const Point &p3, const Point &p4) const { return geom_traits().construct_tetrahedron_3_object()(p1,p2,p3,p4); } Tetrahedron construct_tetrahedron( const Point &p1, const Point &p2, const Point &p3, const Point &p4, const Offset &o1, const Offset &o2, const Offset &o3, const Offset &o4) const { return geom_traits().construct_tetrahedron_3_object()(p1,p2,p3,p4, o1,o2,o3,o4); } Tetrahedron construct_tetrahedron(const Periodic_tetrahedron& tet) { return construct_tetrahedron( tet[0].first, tet[1].first, tet[2].first, tet[3].first, tet[0].second, tet[1].second, tet[2].second, tet[3].second); } Triangle construct_triangle( const Point &p1, const Point &p2, const Point &p3) const { return geom_traits().construct_triangle_3_object()(p1,p2,p3); } Triangle construct_triangle( const Point &p1, const Point &p2, const Point &p3, const Offset &o1, const Offset &o2, const Offset &o3) const { return geom_traits().construct_triangle_3_object()(p1,p2,p3,o1,o2,o3); } Triangle construct_triangle(const Periodic_triangle& tri) { return construct_triangle(tri[0].first, tri[1].first, tri[2].first, tri[0].second, tri[1].second, tri[2].second); } Segment construct_segment(const Point &p1, const Point &p2) const { return geom_traits().construct_segment_3_object()(p1, p2); } Segment construct_segment(const Point &p1, const Point &p2, const Offset &o1, const Offset &o2) const { return geom_traits().construct_segment_3_object()(p1,p2,o1,o2); } Segment construct_segment(const Periodic_segment& seg) const { return construct_segment(seg[0].first, seg[1].first, seg[0].second, seg[1].second); } Point construct_point(const Point& p, const Offset &o) const { return geom_traits().construct_point_3_object()(p,o); } Point construct_point(const Periodic_point& pp) const { return construct_point(pp.first, pp.second); } //@} public: /** @name Geometric access functions */ //@{ Periodic_point periodic_point( const Vertex_handle v ) const { if (is_1_cover()) return std::make_pair(v->point(), Offset(0,0,0)); Virtual_vertex_map_it it = virtual_vertices.find(v); if (it == virtual_vertices.end()) { // if v is not contained in virtual_vertices, then it is in the // original domain. return std::make_pair(v->point(), Offset(0,0,0)); } else { // otherwise it has to be looked up as well as its offset. return std::make_pair(it->second.first->point(), it->second.second); } } Periodic_point periodic_point( const Cell_handle c, int i) const { if (is_1_cover()) return std::make_pair(c->vertex(i)->point(), int_to_off(c->offset(i))); Virtual_vertex_map_it it = virtual_vertices.find(c->vertex(i)); if (it == virtual_vertices.end()) { // if c->vertex(i) is not contained in virtual_vertices, then it // is in the original domain. return std::make_pair(c->vertex(i)->point(), combine_offsets(Offset(),int_to_off(c->offset(i))) ); } else { // otherwise it has to be looked up as well as its offset. return std::make_pair(it->second.first->point(), combine_offsets(it->second.second, int_to_off(c->offset(i))) ); } } Periodic_segment periodic_segment(const Cell_handle c, int i, int j) const { CGAL_triangulation_precondition( i != j ); CGAL_triangulation_precondition( number_of_vertices() != 0 ); CGAL_triangulation_precondition( i >= 0 && i <= 3 && j >= 0 && j <= 3 ); return make_array( std::make_pair(c->vertex(i)->point(), get_offset(c,i)), std::make_pair(c->vertex(j)->point(), get_offset(c,j)) ); } Periodic_segment periodic_segment(const Edge & e) const { return periodic_segment(e.first,e.second,e.third); } Periodic_triangle periodic_triangle(const Cell_handle c, int i) const; Periodic_triangle periodic_triangle(const Facet & f) const { return periodic_triangle(f.first, f.second); } Periodic_tetrahedron periodic_tetrahedron(const Cell_handle c) const { CGAL_triangulation_precondition( number_of_vertices() != 0 ); return make_array( std::make_pair(c->vertex(0)->point(), get_offset(c,0)), std::make_pair(c->vertex(1)->point(), get_offset(c,1)), std::make_pair(c->vertex(2)->point(), get_offset(c,2)), std::make_pair(c->vertex(3)->point(), get_offset(c,3)) ); } Point point(const Periodic_point & pp) const { return construct_point(pp.first, pp.second); } Segment segment(const Periodic_segment & ps) const { return construct_segment(ps[0].first,ps[1].first,ps[0].second,ps[1].second); } Triangle triangle(const Periodic_triangle & pt) const { return construct_triangle(pt[0].first, pt[1].first, pt[2].first, pt[0].second,pt[1].second,pt[2].second); } Tetrahedron tetrahedron(const Periodic_tetrahedron & pt) const { return construct_tetrahedron(pt[0].first, pt[1].first, pt[2].first, pt[3].first, pt[0].second,pt[1].second, pt[2].second,pt[3].second); } // @} /** @name Queries */ //@{ bool is_vertex(const Point & p, Vertex_handle & v) const; bool is_vertex(Vertex_handle v) const { return _tds.is_vertex(v); } bool is_edge(Vertex_handle u, Vertex_handle v, Cell_handle & c, int & i, int & j) const { return _tds.is_edge(u, v, c, i, j); } bool is_edge(Vertex_handle u, const Offset & off_u, Vertex_handle v, const Offset & off_v, Cell_handle & c, int & i, int & j) const { if (!_tds.is_edge(u,v,c,i,j)) return false; if ((get_offset(c,i) == off_u) && (get_offset(c,j) == off_v)) return true; // it might be that different cells containing (u,v) yield // different offsets, which forces us to test for all possibilities. else { Cell_circulator ccirc = incident_cells(c,i,j,c); while (++ccirc != c) { i = ccirc->index(u); j = ccirc->index(v); if ((get_offset(ccirc,i) == off_u) && (get_offset(ccirc,j) == off_v)) { c = ccirc; return true; } } return false; } } bool is_facet(Vertex_handle u, Vertex_handle v, Vertex_handle w, Cell_handle & c, int & i, int & j, int & k) const { return _tds.is_facet(u, v, w, c, i, j, k); } bool is_facet(Vertex_handle u, const Offset & off_u, Vertex_handle v, const Offset & off_v, Vertex_handle w, const Offset & off_w, Cell_handle & c, int & i, int & j, int & k) const { if (!_tds.is_facet(u,v,w,c,i,j,k)) return false; if ((get_offset(c,i) == off_u) && (get_offset(c,j) == off_v) && (get_offset(c,k) == off_w) ) return true; // it might be that c and c->neighbor(l) yield different offsets // which forces us to test for both possibilities. else { int l = 6-i-j-k; c = c->neighbor(l); i = c->index(u); j = c->index(v); k = c->index(w); return ((get_offset(c,i) == off_u) && (get_offset(c,j) == off_v) && (get_offset(c,k) == off_w) ); } } bool is_cell(Cell_handle c) const { return _tds.is_cell(c); } bool is_cell(Vertex_handle u, Vertex_handle v, Vertex_handle w, Vertex_handle t, Cell_handle & c, int & i, int & j, int & k, int & l) const { return _tds.is_cell(u, v, w, t, c, i, j, k, l); } bool is_cell(Vertex_handle u, Vertex_handle v, Vertex_handle w, Vertex_handle t, Cell_handle & c) const { int i,j,k,l; return _tds.is_cell(u, v, w, t, c, i, j, k, l); } bool is_cell(Vertex_handle u, const Offset & off_u, Vertex_handle v, const Offset & off_v, Vertex_handle w, const Offset & off_w, Vertex_handle t, const Offset & off_t, Cell_handle & c, int & i, int & j, int & k, int & l) const { if (!_tds.is_cell(u,v,w,t,c,i,j,k,l)) return false; return ((get_offset(c,i) == off_u) && (get_offset(c,j) == off_v) && (get_offset(c,k) == off_w) && (get_offset(c,l) == off_t) ); return false; } bool is_cell(Vertex_handle u, const Offset & off_u, Vertex_handle v, const Offset & off_v, Vertex_handle w, const Offset & off_w, Vertex_handle t, const Offset & off_t, Cell_handle & c) const { int i, j, k, l; return is_cell(u,off_u,v,off_v,w,off_w,t,off_t,c,i,j,k,l); } bool has_vertex(const Facet & f, Vertex_handle v, int & j) const { return _tds.has_vertex(f.first, f.second, v, j); } bool has_vertex(Cell_handle c, int i, Vertex_handle v, int & j) const { return _tds.has_vertex(c, i, v, j); } bool has_vertex(const Facet & f, Vertex_handle v) const { return _tds.has_vertex(f.first, f.second, v); } bool has_vertex(Cell_handle c, int i, Vertex_handle v) const { return _tds.has_vertex(c, i, v); } bool are_equal(Cell_handle c, int i, Cell_handle n, int j) const { return _tds.are_equal(c, i, n, j); } bool are_equal(const Facet & f, const Facet & g) const { return _tds.are_equal(f.first, f.second, g.first, g.second); } bool are_equal(const Facet & f, Cell_handle n, int j) const { return _tds.are_equal(f.first, f.second, n, j); } //@} #ifdef CGAL_NO_STRUCTURAL_FILTERING Cell_handle periodic_locate(const Point & p, const Offset &o_p, Locate_type & lt, int & li, int & lj, Cell_handle start = Cell_handle()) const; #else // no CGAL_NO_STRUCTURAL_FILTERING # ifndef CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS # define CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS 2500 # endif // no CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS public: Cell_handle inexact_periodic_locate(const Point& p, const Offset &o_p, Cell_handle start = Cell_handle(), int max_num_cells = CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const; protected: Cell_handle exact_periodic_locate(const Point& p, const Offset &o_p, Locate_type& lt, int& li, int & lj, Cell_handle start) const; Cell_handle generic_periodic_locate(const Point& p, const Offset &o_p, Locate_type& lt, int& li, int & lj, Cell_handle start, internal::Periodic_structural_filtering_3_tag) const { return exact_periodic_locate(p, o_p, lt, li, lj, inexact_periodic_locate(p, o_p, start)); } Cell_handle generic_periodic_locate(const Point& p, const Offset &o_p, Locate_type& lt, int& li, int & lj, Cell_handle start, internal::No_periodic_structural_filtering_3_tag) const { return exact_periodic_locate(p, o_p, lt, li, lj, start); } Orientation inexact_orientation(const Point &p, const Point &q, const Point &r, const Point &s) const { const double px = to_double(p.x()); const double py = to_double(p.y()); const double pz = to_double(p.z()); const double qx = to_double(q.x()); const double qy = to_double(q.y()); const double qz = to_double(q.z()); const double rx = to_double(r.x()); const double ry = to_double(r.y()); const double rz = to_double(r.z()); const double sx = to_double(s.x()); const double sy = to_double(s.y()); const double sz = to_double(s.z()); const double pqx = qx - px; const double pqy = qy - py; const double pqz = qz - pz; const double prx = rx - px; const double pry = ry - py; const double prz = rz - pz; const double psx = sx - px; const double psy = sy - py; const double psz = sz - pz; const double det = determinant(pqx, pqy, pqz, prx, pry, prz, psx, psy, psz); if (det > 0) return POSITIVE; if (det < 0) return NEGATIVE; return ZERO; } Orientation inexact_orientation(const Point &p, const Point &q, const Point &r, const Point &s, const Offset& o_p, const Offset& o_q, const Offset& o_r, const Offset& o_s) const { return inexact_orientation(construct_point(p, o_p), construct_point(q, o_q), construct_point(r, o_r), construct_point(s, o_s)); } public: Cell_handle periodic_locate(const Point & p, const Offset &o_p, Locate_type & lt, int & li, int & lj, Cell_handle start = Cell_handle()) const { typedef Triangulation_structural_filtering_traits TSFT; typedef typename internal::Periodic_structural_filtering_selector_3< TSFT::Use_structural_filtering_tag::value >::Tag Should_filter_tag; return generic_periodic_locate(p, o_p, lt, li, lj, start, Should_filter_tag()); } Cell_handle inexact_locate(const Point& p, Cell_handle start = Cell_handle(), int max_num_cells = CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const { return inexact_periodic_locate(p, Offset(), start, max_num_cells); } #endif // no CGAL_NO_STRUCTURAL_FILTERING protected: /** @name Location helpers */ //@{ // Cell_handle periodic_locate(const Point & p, const Offset &o_p, // Locate_type & lt, int & li, int & lj, Cell_handle start) const; Bounded_side side_of_cell(const Point & p, const Offset &off, Cell_handle c, Locate_type & lt, int & i, int & j) const; //@} public: /** @name Point Location */ //@{ /** Wrapper function for locate if only the request point is given. */ Cell_handle locate(const Point & p, Cell_handle start=Cell_handle()) const { Locate_type lt; int li, lj; return locate( p, lt, li, lj, start); } /** Wrapper function calling locate with an empty offset if there was no * offset given. */ Cell_handle locate(const Point & p, Locate_type & lt, int & li, int & lj, Cell_handle start = Cell_handle()) const { return periodic_locate(p, Offset(), lt, li, lj, start); } Bounded_side side_of_cell(const Point & p, Cell_handle c, Locate_type & lt, int & i, int & j) const { if (number_of_vertices() == 0) { lt = EMPTY; return ON_UNBOUNDED_SIDE; } return side_of_cell(p,Offset(),c,lt,i,j); } //@} private: /** @name Insertion helpers */ //@{ template void insert_too_long_edges(Vertex_handle v, const CellIt begin, const CellIt end); template void delete_too_long_edges(const CellIt begin, const CellIt end); template < class Conflict_tester, class Point_hider > Vertex_handle periodic_insert(const Point& p, const Offset& o, Locate_type lt, Cell_handle c, const Conflict_tester &tester, Point_hider &hider, Vertex_handle vh = Vertex_handle()); template void periodic_sort(Point_iterator /*p_begin*/, Point_iterator /*p_end*/, Offset_iterator /*o_begin*/, Offset_iterator /*o_end*/) const { std::cout << "Periodic_sort not yet implemented" << std::endl; } Vertex_handle create_initial_triangulation(const Point &p); public: std::vector insert_dummy_points(); protected: // this is needed for compatibility reasons template Triple find_conflicts(Cell_handle c, const Conflict_test &tester, Triple it) const { Offset off = get_location_offset(tester, c); return find_conflicts(c,off,tester,it); } template Triple find_conflicts(Cell_handle c, const Offset ¤t_off, const Conflict_test &tester, Triple it) const; //@} protected: // COMMON INSERTION for DELAUNAY and REGULAR TRIANGULATION template < class Conflict_tester, class Point_hider > Vertex_handle insert_in_conflict(const Point & p, Cell_handle start, const Conflict_tester &tester, Point_hider &hider) { Locate_type lt = Locate_type(); int li=0, lj=0; Cell_handle c = periodic_locate(p, Offset(), lt, li, lj, start); return insert_in_conflict(p,lt,c,li,lj,tester,hider); } template < class Conflict_tester, class Point_hider > Vertex_handle insert_in_conflict(const Point & p, Locate_type lt, Cell_handle c, int li, int lj, const Conflict_tester &tester, Point_hider &hider); template < class InputIterator, class Conflict_tester, class Point_hider> std::vector insert_in_conflict( InputIterator begin, InputIterator end, Cell_handle start, Conflict_tester &tester, Point_hider &hider) { Vertex_handle new_vertex; std::vector double_vertices; Locate_type lt = Locate_type(); int li=0, lj=0; CGAL_triangulation_assertion_code( Locate_type lta = Locate_type(); ) CGAL_triangulation_assertion_code( int ia = 0; ) CGAL_triangulation_assertion_code( int ja = 0; ) Cell_handle hint; while (begin!=end) { tester.set_point(*begin); hint = periodic_locate(*begin, Offset(), lt, li, lj, start); CGAL_triangulation_assertion_code( if (number_of_vertices() != 0) { ); CGAL_triangulation_assertion(side_of_cell( *begin,Offset(), hint, lta, ia, ja) != ON_UNBOUNDED_SIDE); CGAL_triangulation_assertion(lta == lt); CGAL_triangulation_assertion(ia == li); CGAL_triangulation_assertion(ja == lj); CGAL_triangulation_assertion_code( } ); new_vertex = insert_in_conflict(*begin,lt,hint,li,lj,tester,hider); if (lt == VERTEX) double_vertices.push_back(new_vertex); start = new_vertex->cell(); begin++; } return double_vertices; } //@} private: /** @name Removal helpers */ //@{ Vertex_triple make_vertex_triple(const Facet& f) const { Cell_handle ch = f.first; int i = f.second; return Vertex_triple(ch->vertex(vertex_triple_index(i,0)), ch->vertex(vertex_triple_index(i,1)), ch->vertex(vertex_triple_index(i,2))); } void make_canonical(Vertex_triple& t) const; void make_hole(Vertex_handle v, std::map &outer_map, std::vector &hole); template < class PointRemover > void periodic_remove(Vertex_handle v, PointRemover &remover); //@} protected: /** @name Removal */ //@{ template < class PointRemover, class CT > void remove(Vertex_handle v, PointRemover &remover, CT &ct); //@} public: /** @name Traversal */ //@{ Cell_iterator cells_begin() const { return _tds.cells_begin(); } Cell_iterator cells_end() const { return _tds.cells_end(); } Vertex_iterator vertices_begin() const { return _tds.vertices_begin(); } Vertex_iterator vertices_end() const { return _tds.vertices_end(); } Edge_iterator edges_begin() const { return _tds.edges_begin(); } Edge_iterator edges_end() const { return _tds.edges_end(); } Facet_iterator facets_begin() const { return _tds.facets_begin(); } Facet_iterator facets_end() const { return _tds.facets_end(); } Cell_iterator finite_cells_begin() const { return _tds.cells_begin(); } Cell_iterator finite_cells_end() const { return _tds.cells_end(); } Vertex_iterator finite_vertices_begin() const { return _tds.vertices_begin(); } Vertex_iterator finite_vertices_end() const { return _tds.vertices_end(); } Edge_iterator finite_edges_begin() const { return _tds.edges_begin(); } Edge_iterator finite_edges_end() const { return _tds.edges_end(); } Facet_iterator finite_facets_begin() const { return _tds.facets_begin(); } Facet_iterator finite_facets_end() const { return _tds.facets_end(); } All_cells_iterator all_cells_begin() const { return _tds.cells_begin(); } All_cells_iterator all_cells_end() const { return _tds.cells_end(); } All_vertices_iterator all_vertices_begin() const { return _tds.vertices_begin(); } All_vertices_iterator all_vertices_end() const { return _tds.vertices_end(); } All_edges_iterator all_edges_begin() const { return _tds.edges_begin(); } All_edges_iterator all_edges_end() const { return _tds.edges_end(); } All_facets_iterator all_facets_begin() const { return _tds.facets_begin(); } All_facets_iterator all_facets_end() const { return _tds.facets_end(); } Unique_vertex_iterator unique_vertices_begin() const { return CGAL::filter_iterator(vertices_end(), Domain_tester(this), vertices_begin()); } Unique_vertex_iterator unique_vertices_end() const { return CGAL::filter_iterator(vertices_end(), Domain_tester(this)); } // Geometric iterators Periodic_tetrahedron_iterator periodic_tetrahedra_begin( Iterator_type it = STORED) const { return Periodic_tetrahedron_iterator(this, it); } Periodic_tetrahedron_iterator periodic_tetrahedra_end( Iterator_type it = STORED) const { return Periodic_tetrahedron_iterator(this, 1, it); } Periodic_triangle_iterator periodic_triangles_begin( Iterator_type it = STORED) const { return Periodic_triangle_iterator(this, it); } Periodic_triangle_iterator periodic_triangles_end( Iterator_type it = STORED) const { return Periodic_triangle_iterator(this, 1, it); } Periodic_segment_iterator periodic_segments_begin( Iterator_type it = STORED) const { return Periodic_segment_iterator(this, it); } Periodic_segment_iterator periodic_segments_end( Iterator_type it = STORED) const { return Periodic_segment_iterator(this, 1, it); } Periodic_point_iterator periodic_points_begin( Iterator_type it = STORED) const { return Periodic_point_iterator(this, it); } Periodic_point_iterator periodic_points_end( Iterator_type it = STORED) const { return Periodic_point_iterator(this, 1, it); } // Circulators Cell_circulator incident_cells(const Edge & e) const { return _tds.incident_cells(e); } Cell_circulator incident_cells(Cell_handle c, int i, int j) const { return _tds.incident_cells(c, i, j); } Cell_circulator incident_cells(const Edge & e, Cell_handle start) const { return _tds.incident_cells(e, start); } Cell_circulator incident_cells(Cell_handle c, int i, int j, Cell_handle start) const { return _tds.incident_cells(c, i, j, start); } Facet_circulator incident_facets(const Edge & e) const { return _tds.incident_facets(e); } Facet_circulator incident_facets(Cell_handle c, int i, int j) const { return _tds.incident_facets(c, i, j); } Facet_circulator incident_facets(const Edge & e, const Facet & start) const { return _tds.incident_facets(e, start); } Facet_circulator incident_facets(Cell_handle c, int i, int j, const Facet & start) const { return _tds.incident_facets(c, i, j, start); } Facet_circulator incident_facets(const Edge & e, Cell_handle start, int f) const { return _tds.incident_facets(e, start, f); } Facet_circulator incident_facets(Cell_handle c, int i, int j, Cell_handle start, int f) const { return _tds.incident_facets(c, i, j, start, f); } // around a vertex template OutputIterator incident_cells(Vertex_handle v, OutputIterator cells) const { return _tds.incident_cells(v, cells); } template OutputIterator incident_facets(Vertex_handle v, OutputIterator facets) const { return _tds.incident_facets(v, facets); } template OutputIterator incident_edges( Vertex_handle v, OutputIterator edges) const { return _tds.incident_edges(v, edges); } template OutputIterator adjacent_vertices( Vertex_handle v, OutputIterator vertices) const { return _tds.adjacent_vertices(v, vertices); } //deprecated, don't use anymore template OutputIterator incident_vertices( Vertex_handle v, OutputIterator vertices) const { return _tds.adjacent_vertices(v, vertices); } size_type degree(Vertex_handle v) const { return _tds.degree(v); } // Functions forwarded from TDS. int mirror_index(Cell_handle c, int i) const { return _tds.mirror_index(c, i); } Vertex_handle mirror_vertex(Cell_handle c, int i) const { return _tds.mirror_vertex(c, i); } Facet mirror_facet(Facet f) const { return _tds.mirror_facet(f); } //@} private: /** @name Checking helpers */ //@{ /// calls has_self_edges for every cell of the triangulation bool has_self_edges() const { Cell_iterator it; for ( it = all_cells_begin(); it != all_cells_end(); ++it ) if (has_self_edges(it)) return true; return false; } bool has_self_edges(Cell_handle c) const; //@} public: /** @name Checking */ //@{ bool is_valid(bool verbose = false, int level = 0) const; bool is_valid(Cell_handle c, bool verbose = false, int level = 0) const; protected: template bool is_valid_conflict(ConflictTester &tester, bool verbose = false, int level = 0) const; //@} protected: /** @name Functors */ //@{ template < class Cmp > class Perturbation_order; //@} public: // undocumented access functions Offset get_offset(Cell_handle ch, int i) const { if (is_1_cover()) return int_to_off(ch->offset(i)); Virtual_vertex_map_it it = virtual_vertices.find(ch->vertex(i)); if (it != virtual_vertices.end()) return combine_offsets(it->second.second, int_to_off(ch->offset(i))); else return combine_offsets(Offset(), int_to_off(ch->offset(i))); } Offset get_offset(Vertex_handle vh) const { if (is_1_cover()) return Offset(); Virtual_vertex_map_it it = virtual_vertices.find(vh); if (it != virtual_vertices.end()) return it->second.second; else return Offset(); } Vertex_handle get_original_vertex(Vertex_handle vh) const { if (is_1_cover()) return vh; Virtual_vertex_map_it it = virtual_vertices.find(vh); if (it != virtual_vertices.end()) return it->second.first; else return vh; } Offset combine_offsets(const Offset& o_c, const Offset& o_t) const { Offset o_ct(_cover[0]*o_t.x(),_cover[1]*o_t.y(),_cover[2]*o_t.z()); return o_c + o_ct; } // These functions give the pair (vertex, offset) that corresponds to the // i-th vertex of cell ch or vertex vh, respectively. void get_vertex(Cell_handle ch, int i, Vertex_handle &vh, Offset &off) const; void get_vertex(Vertex_handle vh_i, Vertex_handle &vh, Offset &off) const; protected: // Auxiliary functions int find_too_long_edges(std::map >& edges) const; Cell_handle get_cell(const Vertex_handle* vh) const; template Offset get_location_offset(const Conflict_tester& tester, Cell_handle c) const; Offset get_neighbor_offset(Cell_handle ch, int i, Cell_handle nb) const; /** @name Friends */ //@{ friend class Perturbation_order; friend std::istream& operator>> <> (std::istream& is, Periodic_3_triangulation_3 &tr); friend std::ostream& operator<< <> (std::ostream& os, const Periodic_3_triangulation_3 &tr); //@} // unused and undocumented function required to be compatible to // Alpha_shape_3 public: Point point(Cell_handle c, int idx) const { // if (is_1_cover()) return point(periodic_point(c,idx)); Offset vec_off[4]; for (int i=0 ; i<4 ; i++) vec_off[i] = periodic_point(c,i).second; int ox = vec_off[0].x(); int oy = vec_off[0].y(); int oz = vec_off[0].z(); for (int i=1 ; i<4 ; i++) { ox = (std::min)(ox,vec_off[i].x()); oy = (std::min)(oy,vec_off[i].y()); oz = (std::min)(oz,vec_off[i].z()); } Offset diff_off(-ox,-oy,-oz); if (diff_off.is_null()) return point(periodic_point(c,idx)); for (unsigned int i=0 ; i<4 ; i++) vec_off[i] += diff_off; Vertex_handle canonic_vh[4]; for (int i=0 ; i<4 ; i++) { Virtual_vertex_map_it vvmit = virtual_vertices.find(c->vertex(i)); Vertex_handle orig_vh; if (vvmit == virtual_vertices.end()) orig_vh = c->vertex(i); else orig_vh = vvmit->second.first; if (vec_off[i].is_null()) canonic_vh[i] = orig_vh; else { CGAL_assertion(virtual_vertices_reverse.find(orig_vh) != virtual_vertices_reverse.end()); canonic_vh[i] = virtual_vertices_reverse.find(orig_vh) ->second[9*vec_off[i][0]+3*vec_off[i][1]+vec_off[i][2]-1]; } } std::vector cells; incident_cells(canonic_vh[0], std::back_inserter(cells)); for (unsigned int i=0 ; ihas_vertex(canonic_vh[0])); if (cells[i]->has_vertex(canonic_vh[1]) && cells[i]->has_vertex(canonic_vh[2]) && cells[i]->has_vertex(canonic_vh[3]) ) return point(periodic_point(cells[i],cells[i]->index(canonic_vh[idx]))); } CGAL_assertion(false); return Point(); } }; template < class GT, class TDS > inline void Periodic_3_triangulation_3:: copy_multiple_covering(const Periodic_3_triangulation_3 & tr) { // Write the respective offsets in the vertices to make them // automatically copy with the tds. for (Vertex_iterator vit = tr.vertices_begin() ; vit != tr.vertices_end() ; ++vit) { vit->set_offset(tr.get_offset(vit)); } // copy the tds _tds = tr.tds(); // make a list of all vertices that belong to the original // domain and initialize the basic structure of // virtual_vertices_reverse std::list vlist; for (Vertex_iterator vit = vertices_begin() ; vit != vertices_end() ; ++vit) { if (vit->offset() == Offset()) { vlist.push_back(vit); virtual_vertices_reverse.insert( std::make_pair(vit,std::vector(26))); CGAL_triangulation_assertion(virtual_vertices_reverse.find(vit) ->second.size() == 26); } } // Iterate over all vertices that are not in the original domain // and construct the respective entries to virtual_vertices and // virtual_vertices_reverse for (Vertex_iterator vit2 = vertices_begin() ; vit2 != vertices_end() ; ++vit2) { if (vit2->offset() != Offset()) { //TODO: use some binding, maybe boost instead of the Finder. typename std::list::iterator vlist_it = std::find_if(vlist.begin(), vlist.end(), Finder(this,vit2->point())); Offset off = vit2->offset(); virtual_vertices.insert(std::make_pair(vit2, std::make_pair(*vlist_it,off))); virtual_vertices_reverse.find(*vlist_it) ->second[9*off[0]+3*off[1]+off[2]-1]=vit2; CGAL_triangulation_assertion(get_offset(vit2) == off); } } // Cleanup vertex offsets for (Vertex_iterator vit = vertices_begin() ; vit != vertices_end() ; ++vit) vit->clear_offset(); for (Vertex_iterator vit = tr.vertices_begin() ; vit != tr.vertices_end() ; ++vit) vit->clear_offset(); // Build up the too_long_edges container too_long_edge_counter = 0; too_long_edges.clear(); for (Vertex_iterator vit = vertices_begin() ; vit != vertices_end() ; ++vit) too_long_edges[vit] = std::list(); std::pair edge_to_add; Point p1,p2; int i,j; for (Edge_iterator eit = edges_begin() ; eit != edges_end() ; ++eit) { if (&*(eit->first->vertex(eit->second)) < &*(eit->first->vertex(eit->third))) { i = eit->second; j = eit->third; } else { i = eit->third; j = eit->second; } edge_to_add = std::make_pair(eit->first->vertex(i), eit->first->vertex(j)); p1 = construct_point(eit->first->vertex(i)->point(), get_offset(eit->first, i)); p2 = construct_point(eit->first->vertex(j)->point(), get_offset(eit->first, j)); Vertex_handle v_no = eit->first->vertex(i); if (squared_distance(p1,p2) > edge_length_threshold) { CGAL_triangulation_assertion( find(too_long_edges[v_no].begin(), too_long_edges[v_no].end(), edge_to_add.second) == too_long_edges[v_no].end()); too_long_edges[v_no].push_back(edge_to_add.second); too_long_edge_counter++; } } } template < class GT, class TDS > inline bool Periodic_3_triangulation_3:: is_extensible_triangulation_in_1_sheet_h1() const { if (!is_1_cover()) { if (too_long_edge_counter == 0) return true; else return false; } else { typename Geometric_traits::FT longest_edge_squared_length(0); Segment s; for (Periodic_segment_iterator psit = periodic_segments_begin(UNIQUE); psit != periodic_segments_end(UNIQUE) ; ++psit) { s = construct_segment(*psit); longest_edge_squared_length = (std::max)(longest_edge_squared_length, s.squared_length()); } return (longest_edge_squared_length < edge_length_threshold); } } template < class GT, class TDS > inline bool Periodic_3_triangulation_3:: is_extensible_triangulation_in_1_sheet_h2() const { typedef typename Geometric_traits::Construct_circumcenter_3 Construct_circumcenter; typedef typename Geometric_traits::FT FT; Construct_circumcenter construct_circumcenter = _gt.construct_circumcenter_3_object(); for (Periodic_tetrahedron_iterator tit = periodic_tetrahedra_begin(UNIQUE) ; tit != periodic_tetrahedra_end(UNIQUE) ; ++tit) { Point cc = construct_circumcenter( tit->at(0).first, tit->at(1).first, tit->at(2).first, tit->at(3).first, tit->at(0).second, tit->at(1).second, tit->at(2).second, tit->at(3).second); if ( !(FT(16)*squared_distance(cc,point(tit->at(0))) < (_domain.xmax()-_domain.xmin())*(_domain.xmax()-_domain.xmin())) ) return false; } return true; } template < class GT, class TDS > inline bool Periodic_3_triangulation_3:: is_triangulation_in_1_sheet() const { if (is_1_cover()) return true; for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) { if (virtual_vertices.find(vit) == virtual_vertices.end()) continue; std::vector nb_v; std::set nb_v_odom; Vertex_handle vh; Offset off; adjacent_vertices(vit, std::back_inserter(nb_v)); for (unsigned int i=0; i inline bool Periodic_3_triangulation_3:: is_vertex(const Point & p, Vertex_handle & v) const { Locate_type lt; int li, lj; Cell_handle c = locate( p, lt, li, lj ); if ( lt != VERTEX ) return false; v = c->vertex(li); return true; } template < class GT, class TDS > inline void Periodic_3_triangulation_3:: make_canonical(Vertex_triple& t) const { int i = (&*(t.first) < &*(t.second))? 0 : 1; if(i==0) { i = (&*(t.first) < &*(t.third))? 0 : 2; } else { i = (&*(t.second) < &*(t.third))? 1 : 2; } Vertex_handle tmp; switch(i){ case 0: return; case 1: tmp = t.first; t.first = t.second; t.second = t.third; t.third = tmp; return; default: tmp = t.first; t.first = t.third; t.third = t.second; t.second = tmp; } } template < class GT, class TDS > inline typename Periodic_3_triangulation_3::Periodic_triangle Periodic_3_triangulation_3:: periodic_triangle(const Cell_handle c, int i) const { CGAL_triangulation_precondition( number_of_vertices() != 0 ); CGAL_triangulation_precondition( i >= 0 && i <= 3 ); if ( (i&1)==0 ) return make_array(std::make_pair(c->vertex( (i+2)&3 )->point(), get_offset(c,(i+2)&3)), std::make_pair(c->vertex( (i+1)&3 )->point(), get_offset(c,(i+1)&3)), std::make_pair(c->vertex( (i+3)&3 )->point(), get_offset(c,(i+3)&3)) ); return make_array(std::make_pair(c->vertex( (i+1)&3 )->point(), get_offset(c,(i+1)&3)), std::make_pair(c->vertex( (i+2)&3 )->point(), get_offset(c,(i+2)&3)), std::make_pair(c->vertex( (i+3)&3 )->point(), get_offset(c,(i+3)&3)) ); } /** Assumes a point, an offset, and a cell to start from. * Gives the locate type and the simplex (cell and indices) containing p. * * Performs a remembering stochastic walk if the triangulation is not empty. * After the walk the type of the simplex containing p is determined. * * returns the cell p lies in * starts at cell "start" * returns a cell Cell_handel if lt == CELL * returns a facet (Cell_handle,li) if lt == FACET * returns an edge (Cell_handle,li,lj) if lt == EDGE * returns a vertex (Cell_handle,li) if lt == VERTEX */ template < class GT, class TDS > inline typename Periodic_3_triangulation_3::Cell_handle Periodic_3_triangulation_3:: #ifdef CGAL_NO_STRUCTURAL_FILTERING periodic_locate #else exact_periodic_locate #endif (const Point & p, const Offset &o_p, Locate_type & lt, int & li, int & lj, Cell_handle start) const { int cumm_off = 0; Offset off_query = o_p; if (number_of_vertices() == 0) { lt = EMPTY; return Cell_handle(); } CGAL_triangulation_assertion(number_of_vertices() != 0); if (start == Cell_handle()) { start = cells_begin(); } cumm_off = start->offset(0) | start->offset(1) | start->offset(2) | start->offset(3); if (is_1_cover() && cumm_off != 0) { if (((cumm_off & 4) == 4) && (FT(2)*p.x()<(_domain.xmax()+_domain.xmin()))) off_query += Offset(1,0,0); if (((cumm_off & 2) == 2) && (FT(2)*p.y()<(_domain.ymax()+_domain.ymin()))) off_query += Offset(0,1,0); if (((cumm_off & 1) == 1) && (FT(2)*p.z()<(_domain.zmax()+_domain.zmin()))) off_query += Offset(0,0,1); } CGAL_triangulation_postcondition(start!=Cell_handle()); CGAL_triangulation_assertion(start->neighbor(0)->neighbor( start->neighbor(0)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(1)->neighbor( start->neighbor(1)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(2)->neighbor( start->neighbor(2)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(3)->neighbor( start->neighbor(3)->index(start))==start); // We implement the remembering visibility/stochastic walk. // Remembers the previous cell to avoid useless orientation tests. Cell_handle previous = Cell_handle(); Cell_handle c = start; // Stores the results of the 4 orientation tests. It will be used // at the end to decide if p lies on a face/edge/vertex/interior. Orientation o[4]; boost::rand48 rng; boost::uniform_smallint<> four(0, 3); boost::variate_generator > die4(rng, four); // Now treat the cell c. try_next_cell: // For the remembering stochastic walk, // we need to start trying with a random index : int i = die4(); // For the remembering visibility walk (Delaunay only), we don't : // int i = 0; cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); bool simplicity_criterion = (cumm_off == 0) && (off_query.is_null()); // We know that the 4 vertices of c are positively oriented. // So, in order to test if p is seen outside from one of c's facets, // we just replace the corresponding point by p in the orientation // test. We do this using the arrays below. Offset off[4]; const Point* pts[4] = { &(c->vertex(0)->point()), &(c->vertex(1)->point()), &(c->vertex(2)->point()), &(c->vertex(3)->point()) }; if (!simplicity_criterion && is_1_cover() ) { for (int i=0; i<4; i++) { off[i] = int_to_off(c->offset(i)); } } if (!is_1_cover()) { // Just fetch the vertices of c as points with offsets for (int i=0; i<4; i++) { pts[i] = &(c->vertex(i)->point()); off[i] = get_offset(c,i); } } for (int j=0; j != 4; ++j, i = (i+1)&3) { Cell_handle next = c->neighbor(i); if (previous == next) { o[i] = POSITIVE; continue; } CGAL_triangulation_assertion(next->neighbor(next->index(c)) == c); // We temporarily put p at i's place in pts. const Point* backup = pts[i]; pts[i] = &p; if (simplicity_criterion && is_1_cover() ) { o[i] = orientation(*pts[0], *pts[1], *pts[2], *pts[3]); if ( o[i] != NEGATIVE ) { pts[i] = backup; continue; } } else { Offset backup_off; backup_off = off[i]; off[i] = off_query; o[i] = orientation(*pts[0], *pts[1], *pts[2], *pts[3], off[0], off[1], off[2], off[3]); if ( o[i] != NEGATIVE ) { pts[i] = backup; off[i] = backup_off; continue; } } // Test whether we need to adapt the offset of the query point. // This means, if we get out of the current cover. off_query = combine_offsets(off_query, get_neighbor_offset(c,i,next)); previous = c; c = next; goto try_next_cell; } // Ok, now we have found the cell. It remains to find the dimension of the // intersected simplex. // now p is in c or on its boundary int sum = ( o[0] == COPLANAR ) + ( o[1] == COPLANAR ) + ( o[2] == COPLANAR ) + ( o[3] == COPLANAR ); switch (sum) { case 0: lt = CELL; break; case 1: lt = FACET; li = ( o[0] == COPLANAR ) ? 0 : ( o[1] == COPLANAR ) ? 1 : ( o[2] == COPLANAR ) ? 2 : 3; break; case 2: lt = EDGE; li = ( o[0] != COPLANAR ) ? 0 : ( o[1] != COPLANAR ) ? 1 : 2; lj = ( o[li+1] != COPLANAR ) ? li+1 : ( o[li+2] != COPLANAR ) ? li+2 : li+3; break; case 3: lt = VERTEX; li = ( o[0] != COPLANAR ) ? 0 : ( o[1] != COPLANAR ) ? 1 : ( o[2] != COPLANAR ) ? 2 : 3; break; default: // Vertex can not lie on four facets CGAL_triangulation_assertion(false); } return c; } #ifndef CGAL_NO_STRUCTURAL_FILTERING template < class GT, class TDS > typename Periodic_3_triangulation_3::Cell_handle Periodic_3_triangulation_3:: inexact_periodic_locate(const Point& p, const Offset& o_p, Cell_handle start, int n_of_turns) const { int cumm_off = 0; Offset off_query = o_p; if (number_of_vertices() == 0) { return Cell_handle(); } CGAL_triangulation_assertion(number_of_vertices() != 0); if (start == Cell_handle()) { start = cells_begin(); } cumm_off = start->offset(0) | start->offset(1) | start->offset(2) | start->offset(3); if (is_1_cover() && cumm_off != 0) { if (((cumm_off & 4) == 4) && (FT(2)*p.x()<(_domain.xmax()+_domain.xmin()))) off_query += Offset(1,0,0); if (((cumm_off & 2) == 2) && (FT(2)*p.y()<(_domain.ymax()+_domain.ymin()))) off_query += Offset(0,1,0); if (((cumm_off & 1) == 1) && (FT(2)*p.z()<(_domain.zmax()+_domain.zmin()))) off_query += Offset(0,0,1); } CGAL_triangulation_postcondition(start!=Cell_handle()); CGAL_triangulation_assertion(start->neighbor(0)->neighbor( start->neighbor(0)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(1)->neighbor( start->neighbor(1)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(2)->neighbor( start->neighbor(2)->index(start))==start); CGAL_triangulation_assertion(start->neighbor(3)->neighbor( start->neighbor(3)->index(start))==start); // We implement the remembering visibility/stochastic walk. // Remembers the previous cell to avoid useless orientation tests. Cell_handle previous = Cell_handle(); Cell_handle c = start; // Now treat the cell c. try_next_cell: --n_of_turns; cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); bool simplicity_criterion = (cumm_off == 0) && (off_query.is_null()); // We know that the 4 vertices of c are positively oriented. // So, in order to test if p is seen outside from one of c's facets, // we just replace the corresponding point by p in the orientation // test. We do this using the arrays below. Offset off[4]; const Point* pts[4] = { &(c->vertex(0)->point()), &(c->vertex(1)->point()), &(c->vertex(2)->point()), &(c->vertex(3)->point()) }; if (!simplicity_criterion && is_1_cover() ) { for (int i=0; i<4; i++) { off[i] = int_to_off(c->offset(i)); } } if (!is_1_cover()) { // Just fetch the vertices of c as points with offsets for (int i=0; i<4; i++) { pts[i] = &(c->vertex(i)->point()); off[i] = get_offset(c,i); } } for (int i=0; i != 4; ++i) { Cell_handle next = c->neighbor(i); if (previous == next) { continue; } // We temporarily put p at i's place in pts. const Point* backup = pts[i]; pts[i] = &p; if (simplicity_criterion && is_1_cover() ) { if ( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3]) != NEGATIVE ) { pts[i] = backup; continue; } } else { Offset backup_off; backup_off = off[i]; off[i] = off_query; if ( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3], off[0], off[1], off[2], off[3]) != NEGATIVE ) { pts[i] = backup; off[i] = backup_off; continue; } } // Test whether we need to adapt the offset of the query point. // This means, if we get out of the current cover. off_query = combine_offsets(off_query, get_neighbor_offset(c,i,next)); previous = c; c = next; if (n_of_turns) goto try_next_cell; } return c; } #endif /** * returns * ON_BOUNDED_SIDE if p inside the cell * (for an infinite cell this means that p lies strictly in the half space * limited by its finite facet) * ON_BOUNDARY if p on the boundary of the cell * (for an infinite cell this means that p lies on the *finite* facet) * ON_UNBOUNDED_SIDE if p lies outside the cell * (for an infinite cell this means that p is not in the preceding * two cases) * * lt has a meaning only when ON_BOUNDED_SIDE or ON_BOUNDARY */ // TODO: currently off is not used. It could probably be optimized // using off. template < class GT, class TDS > inline Bounded_side Periodic_3_triangulation_3::side_of_cell( const Point & q, const Offset &off, Cell_handle c, Locate_type & lt, int & i, int & j) const { CGAL_triangulation_precondition( number_of_vertices() != 0 ); Orientation o0,o1,o2,o3; o0 = o1 = o2 = o3 = ZERO; int cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); if ((cumm_off == 0) && (is_1_cover())) { CGAL_triangulation_assertion(off == Offset()); const Point &p0 = c->vertex(0)->point(); const Point &p1 = c->vertex(1)->point(); const Point &p2 = c->vertex(2)->point(); const Point &p3 = c->vertex(3)->point(); if (((o0 = orientation(q ,p1,p2,p3)) == NEGATIVE) || ((o1 = orientation(p0,q ,p2,p3)) == NEGATIVE) || ((o2 = orientation(p0,p1,q ,p3)) == NEGATIVE) || ((o3 = orientation(p0,p1,p2,q )) == NEGATIVE) ) { return ON_UNBOUNDED_SIDE; } } else { // Special case for the periodic space. Offset off_q; Offset offs[4]; const Point *p[4]; for (int i=0; i<4; i++) { p[i] = &(c->vertex(i)->point()); offs[i] = get_offset(c,i); } CGAL_triangulation_assertion(orientation(*p[0], *p[1], *p[2], *p[3], offs[0], offs[1], offs[2], offs[3]) == POSITIVE); bool found=false; for (int i=0; (i<8)&&(!found); i++) { if ((cumm_off | ((~i)&7)) == 7) { o0 = o1 = o2 = o3 = NEGATIVE; off_q = combine_offsets(off, int_to_off(i)); if (((o0 = orientation( q, *p[1], *p[2], *p[3], off_q ,offs[1],offs[2],offs[3])) != NEGATIVE)&& ((o1 = orientation( *p[0], q, *p[2], *p[3], offs[0], off_q,offs[2],offs[3])) != NEGATIVE)&& ((o2 = orientation( *p[0], *p[1], q, *p[3], offs[0],offs[1], off_q,offs[3])) != NEGATIVE)&& ((o3 = orientation( *p[0], *p[1], *p[2], q, offs[0],offs[1],offs[2], off_q)) != NEGATIVE)) { found = true; } } } if (!found) return ON_UNBOUNDED_SIDE; } // now all the oi's are >=0 // sum gives the number of facets p lies on int sum = ( (o0 == ZERO) ? 1 : 0 ) + ( (o1 == ZERO) ? 1 : 0 ) + ( (o2 == ZERO) ? 1 : 0 ) + ( (o3 == ZERO) ? 1 : 0 ); switch (sum) { case 0: { lt = CELL; return ON_BOUNDED_SIDE; } case 1: { lt = FACET; // i = index such that p lies on facet(i) i = ( o0 == ZERO ) ? 0 : ( o1 == ZERO ) ? 1 : ( o2 == ZERO ) ? 2 : 3; return ON_BOUNDARY; } case 2: { lt = EDGE; // i = smallest index such that p does not lie on facet(i) // i must be < 3 since p lies on 2 facets i = ( o0 == POSITIVE ) ? 0 : ( o1 == POSITIVE ) ? 1 : 2; // j = larger index such that p not on facet(j) // j must be > 0 since p lies on 2 facets j = ( o3 == POSITIVE ) ? 3 : ( o2 == POSITIVE ) ? 2 : 1; return ON_BOUNDARY; } case 3: { lt = VERTEX; // i = index such that p does not lie on facet(i) i = ( o0 == POSITIVE ) ? 0 : ( o1 == POSITIVE ) ? 1 : ( o2 == POSITIVE ) ? 2 : 3; return ON_BOUNDARY; } default: { // impossible : cannot be on 4 facets for a real tetrahedron CGAL_triangulation_assertion(false); return ON_BOUNDARY; } } } // side_of_cell template< class GT, class TDS > template< class CellIt > inline void Periodic_3_triangulation_3:: insert_too_long_edges(Vertex_handle v, const CellIt begin, const CellIt end) { CGAL_triangulation_precondition(number_of_vertices() != 0); // add newly added edges to too_long_edges, if necessary. Point p1,p2; Offset omin; std::pair< Vertex_handle, Vertex_handle > edge_to_add; std::pair< Offset, Offset > edge_to_add_off; std::list empty_list; too_long_edges[v] = empty_list; // Iterate over all cells of the new star. for (CellIt it = begin ; it != end ; ++it) { // Consider all possible vertex pairs. for (int k=0; k<4 ; k++) { for (int j=0; j<4 ; j++) { if (j==k) continue; if (&*((*it)->vertex(j)) > &*((*it)->vertex(k))) continue; // make the offsets canonical (wrt. to some notion) // add to too_long_edges, if not yet added and if "too long" CGAL_triangulation_precondition( &*((*it)->vertex(j))< &*((*it)->vertex(k))); edge_to_add = std::make_pair((*it)->vertex(j), (*it)->vertex(k)); p1 = construct_point((*it)->vertex(j)->point(), get_offset(*it, j)); p2 = construct_point((*it)->vertex(k)->point(), get_offset(*it, k)); if ((squared_distance(p1,p2) > edge_length_threshold) && (find(too_long_edges[(*it)->vertex(j)].begin(), too_long_edges[(*it)->vertex(j)].end(), edge_to_add.second) == too_long_edges[(*it)->vertex(j)].end()) ){ too_long_edges[(*it)->vertex(j)].push_back(edge_to_add.second); too_long_edge_counter++; } } } } } template < class GT, class TDS > template < class CellIt > inline void Periodic_3_triangulation_3:: delete_too_long_edges(const CellIt begin, const CellIt end) { std::pair< Vertex_handle, Vertex_handle > edge_to_delete, edge_to_delete2; typename std::list< Vertex_handle >::iterator sit; // Iterate over all cells that are in the star. That means that those cells // are going to be deleted. Therefore, all of them have to be deleted from // too_long_edges, if they are contained in it. for (CellIt it = begin ; it != end ; ++it) { for (int j=0; j<4 ; j++) { for (int k=0; k<4; k++) { if (&*((*it)->vertex(j)) < &*((*it)->vertex(k))) { edge_to_delete = std::make_pair((*it)->vertex(j),(*it)->vertex(k)); } else { edge_to_delete = std::make_pair((*it)->vertex(k),(*it)->vertex(j)); } Vertex_handle v_no = edge_to_delete.first; sit = find(too_long_edges[v_no].begin(), too_long_edges[v_no].end(), edge_to_delete.second); if (sit != too_long_edges[v_no].end()) { too_long_edges[v_no].erase(sit); too_long_edge_counter--; } } } } } /*! \brief Insert point. * * Inserts the point p into the triangulation. It assumes that * the cell c containing p is already known. * * Implementation: * - some precondition checking * - find and mark conflicting cells --> find_conflicts * Conflicting cells are stored in the vector cells. * - backup hidden points * - Delete the edges of the marked cells from too_long_edges * - Insert the new vertex in the hole obtained by removing the * conflicting cells (star-approach) --> _tds._insert_in_hole * - find out about offsets * - Insert the newly added edges that are "too long" * to too_long_edges * - reinsert hidden points */ template < class GT, class TDS > template < class Conflict_tester, class Point_hider > inline typename Periodic_3_triangulation_3::Vertex_handle Periodic_3_triangulation_3::periodic_insert( const Point & p, const Offset& o, Locate_type lt, Cell_handle c, const Conflict_tester &tester, Point_hider &hider, Vertex_handle vh) { Vertex_handle v; CGAL_triangulation_assertion(number_of_vertices() != 0); CGAL_triangulation_precondition_code( Locate_type lt_assert; int i_assert; int j_assert;); CGAL_triangulation_assertion(side_of_cell(tester.point(),o, c, lt_assert, i_assert, j_assert) != ON_UNBOUNDED_SIDE); tester.set_offset(o); // This only holds for Delaunay CGAL_triangulation_assertion(lt != VERTEX); CGAL_USE(lt); // Choose the periodic copy of tester.point() that is inside c. Offset current_off = get_location_offset(tester, c); CGAL_triangulation_assertion(side_of_cell(tester.point(), combine_offsets(o,current_off),c,lt_assert,i_assert,j_assert) != ON_UNBOUNDED_SIDE); // If the new point is not in conflict with its cell, it is hidden. if (!tester.test_initial_cell(c, current_off)) { hider.hide_point(c,p); return Vertex_handle(); } // Ok, we really insert the point now. // First, find the conflict region. std::vector cells; cells.reserve(32); Facet facet; find_conflicts(c, current_off, tester, make_triple(Oneset_iterator(facet), std::back_inserter(cells), Emptyset_iterator())); // Remember the points that are hidden by the conflicting cells, // as they will be deleted during the insertion. hider.set_vertices(cells.begin(), cells.end()); if (!is_1_cover()) delete_too_long_edges(cells.begin(), cells.end()); // Insertion. Attention: facets[0].first MUST be in conflict! // Compute the star and put it into the data structure. // Store the new cells from the star in nbs. v = _tds._insert_in_hole(cells.begin(), cells.end(), facet.first, facet.second); v->set_point(p); //TODO: this could be done within the _insert_in_hole without losing any //time because each cell is visited in any case. //- Do timings to argue to modify _insert_in_conflicts if need be //- Find the modified _insert_in_hole in the branch svn history of TDS std::vector nbs; incident_cells(v, std::back_inserter(nbs)); // For all neighbors of the newly added vertex v: fetch their offsets from // the tester and reset them in the triangulation data structure. for (typename std::vector::iterator cit = nbs.begin(); cit != nbs.end(); cit++) { Offset off[4]; for (int i=0 ; i<4 ; i++) { off[i] = (*cit)->vertex(i)->offset(); } set_offsets(*cit, off[0], off[1], off[2], off[3]); } for (typename std::vector::iterator voit = v_offsets.begin(); voit != v_offsets.end() ; ++voit) { (*voit)->clear_offset(); } v_offsets.clear(); if (vh != Vertex_handle()) { virtual_vertices[v] = Virtual_vertex(vh,o); virtual_vertices_reverse[vh].push_back(v); } if (!is_1_cover()) insert_too_long_edges(v, nbs.begin(), nbs.end()); // Store the hidden points in their new cells. hider.reinsert_vertices(v); return v; } /** Inserts the first point to a triangulation. * * With inserting the first point the 3-sheeted covering is constructed. * So first, the 27 vertices are inserted and are added to virtual_vertices * Then 6*27 cells are created. * Then all links are set. */ template < class GT, class TDS > inline typename Periodic_3_triangulation_3::Vertex_handle Periodic_3_triangulation_3::create_initial_triangulation( const Point &p) { /// Virtual vertices, one per periodic domain Vertex_handle vir_vertices[3][3][3]; /// Virtual cells, 6 per periodic domain Cell_handle cells[3][3][3][6]; // Initialise vertices: vir_vertices[0][0][0] = _tds.create_vertex(); vir_vertices[0][0][0]->set_point(p); virtual_vertices_reverse[vir_vertices[0][0][0]] = std::vector(); for (int i=0; i<_cover[0]; i++) { for (int j=0; j<_cover[1]; j++) { for (int k=0; k<_cover[2]; k++) { if ((i!=0)||(j!=0)||(k!=0)) { // Initialise virtual vertices out of the domain for debugging vir_vertices[i][j][k] = _tds.create_vertex(); vir_vertices[i][j][k]->set_point(p); //+Offset(i,j,k)); virtual_vertices[vir_vertices[i][j][k]] = Virtual_vertex(vir_vertices[0][0][0], Offset(i,j,k)); virtual_vertices_reverse[vir_vertices[0][0][0]].push_back( vir_vertices[i][j][k]); } CGAL_triangulation_assertion(vir_vertices[i][j][k] != Vertex_handle()); CGAL_triangulation_assertion(vir_vertices[0][0][0]->point() == p); } } } // Create cells: for (int i=0; i<_cover[0]; i++) { for (int j=0; j<_cover[1]; j++) { for (int k=0; k<_cover[2]; k++) { for (int l=0; l<6; l++) { // 6 cells per 'cube' cells[i][j][k][l] = _tds.create_cell(); for (int n=0; n<4; n++) CGAL_triangulation_assertion(cells[i][j][k][l] != Cell_handle()); } } } } // set vertex and neighbor information // index to the right vertex: [number of cells][vertex][offset] int vertex_ind[6][4][3] = { { {0, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 0, 0} }, { {1, 1, 0}, {0, 1, 1}, {1, 0, 1}, {1, 1, 1} }, { {1, 0, 0}, {0, 1, 1}, {0, 1, 0}, {0, 0, 1} }, { {1, 0, 0}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1} }, { {1, 0, 0}, {0, 1, 1}, {1, 0, 1}, {1, 1, 0} }, { {1, 0, 0}, {0, 1, 1}, {1, 1, 0}, {0, 1, 0} } }; int neighb_ind[6][4][4] = { { { 0, 0, 0, 2}, { 0,-1, 0, 5}, { 0, 0,-1, 3}, {-1, 0, 0, 4} }, { { 0, 0, 1, 5}, { 1, 0, 0, 2}, { 0, 1, 0, 3}, { 0, 0, 0, 4} }, { {-1, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 3}, { 0, 0, 0, 5} }, { { 0, 0, 1, 0}, { 0,-1, 0, 1}, { 0, 0, 0, 4}, { 0, 0, 0, 2} }, { { 0, 0, 0, 1}, { 1, 0, 0, 0}, { 0, 0, 0, 5}, { 0, 0, 0, 3} }, { { 0, 1, 0, 0}, { 0, 0,-1, 1}, { 0, 0, 0, 2}, { 0, 0, 0, 4} } }; for (int i=0; i<_cover[0]; i++) { for (int j=0; j<_cover[1]; j++) { for (int k=0; k<_cover[2]; k++) { int offset = (i==_cover[0]-1 ? 4 : 0) | (j==_cover[1]-1 ? 2 : 0) | (k==_cover[2]-1 ? 1 : 0); for (int l=0; l<6; l++) { // cell 0: cells[i][j][k][l]->set_vertices( vir_vertices [(i+vertex_ind[l][0][0])%_cover[0]] [(j+vertex_ind[l][0][1])%_cover[1]] [(k+vertex_ind[l][0][2])%_cover[2]], vir_vertices [(i+vertex_ind[l][1][0])%_cover[0]] [(j+vertex_ind[l][1][1])%_cover[1]] [(k+vertex_ind[l][1][2])%_cover[2]], vir_vertices [(i+vertex_ind[l][2][0])%_cover[0]] [(j+vertex_ind[l][2][1])%_cover[1]] [(k+vertex_ind[l][2][2])%_cover[2]], vir_vertices [(i+vertex_ind[l][3][0])%_cover[0]] [(j+vertex_ind[l][3][1])%_cover[1]] [(k+vertex_ind[l][3][2])%_cover[2]]); set_offsets(cells[i][j][k][l], offset & (vertex_ind[l][0][0]*4 + vertex_ind[l][0][1]*2 + vertex_ind[l][0][2]*1), offset & (vertex_ind[l][1][0]*4 + vertex_ind[l][1][1]*2 + vertex_ind[l][1][2]*1), offset & (vertex_ind[l][2][0]*4 + vertex_ind[l][2][1]*2 + vertex_ind[l][2][2]*1), offset & (vertex_ind[l][3][0]*4 + vertex_ind[l][3][1]*2 + vertex_ind[l][3][2]*1)); cells[i][j][k][l]->set_neighbors( cells [(i+_cover[0]+neighb_ind[l][0][0])%_cover[0]] [(j+_cover[1]+neighb_ind[l][0][1])%_cover[1]] [(k+_cover[2]+neighb_ind[l][0][2])%_cover[2]] [ neighb_ind[l][0][3] ], cells [(i+_cover[0]+neighb_ind[l][1][0])%_cover[0]] [(j+_cover[1]+neighb_ind[l][1][1])%_cover[1]] [(k+_cover[2]+neighb_ind[l][1][2])%_cover[2]] [ neighb_ind[l][1][3] ], cells [(i+_cover[0]+neighb_ind[l][2][0])%_cover[0]] [(j+_cover[1]+neighb_ind[l][2][1])%_cover[1]] [(k+_cover[2]+neighb_ind[l][2][2])%_cover[2]] [ neighb_ind[l][2][3] ], cells [(i+_cover[0]+neighb_ind[l][3][0])%_cover[0]] [(j+_cover[1]+neighb_ind[l][3][1])%_cover[1]] [(k+_cover[2]+neighb_ind[l][3][2])%_cover[2]] [ neighb_ind[l][3][3] ] ); } } } } // set pointers from the vertices to incident cells. for (int i=0; i<_cover[0]; i++) { for (int j=0; j<_cover[1]; j++) { for (int k=0; k<_cover[2]; k++) { vir_vertices[i][j][k]->set_cell(cells[i][j][k][0]); } } } _tds.set_dimension(3); // create the base for too_long_edges; CGAL_triangulation_assertion( too_long_edges.empty() ); CGAL_triangulation_assertion(too_long_edge_counter == 0); for (Vertex_iterator vit = vertices_begin() ; vit !=vertices_end() ; ++vit ) too_long_edges[vit] = std::list();; std::vector temp_inc_cells; for (Vertex_iterator vit = vertices_begin() ; vit !=vertices_end() ; ++vit ) { temp_inc_cells.clear(); incident_cells(vit, std::back_inserter(temp_inc_cells)); for (unsigned int i=0 ; iindex(vit); for (int j=0; j<4 ; j++) { if (j==k) continue; if (&*vit > &*(temp_inc_cells[i]->vertex(j))) continue; if ((find(too_long_edges[vit].begin(), too_long_edges[vit].end(), temp_inc_cells[i]->vertex(j)) == too_long_edges[vit].end()) ){ too_long_edges[vit].push_back(temp_inc_cells[i]->vertex(j)); too_long_edge_counter++; } } } } return vir_vertices[0][0][0]; } #include /** finds all cells that are in conflict with the currently added point * (stored in tester). * * The result will be a hole of which the following data is returned: * - boundary facets * - cells * - internal facets. * * c is the current cell, which must be in conflict. * tester is the function object that tests if a cell is in conflict. */ template template Triple Periodic_3_triangulation_3:: find_conflicts(Cell_handle d, const Offset ¤t_off, const Conflict_test &tester, Triple it) const { CGAL_triangulation_precondition( number_of_vertices() != 0 ); CGAL_triangulation_precondition( tester(d, current_off) ); std::stack > cell_stack; cell_stack.push(std::make_pair(d,current_off)); d->tds_data().mark_in_conflict(); *it.second++ = d; do { Cell_handle c = cell_stack.top().first; Offset current_off2 = cell_stack.top().second; cell_stack.pop(); for (int i=0; i< 4; ++i) { Cell_handle test = c->neighbor(i); if (test->tds_data().is_in_conflict()) { if (c < test) { *it.third++ = Facet(c, i); // Internal facet. } continue; // test was already in conflict. } if (test->tds_data().is_clear()) { Offset o_test = current_off2 + get_neighbor_offset(c, i, test); if (tester(test,o_test)) { if (c < test) *it.third++ = Facet(c, i); // Internal facet. cell_stack.push(std::make_pair(test,o_test)); test->tds_data().mark_in_conflict(); *it.second++ = test; continue; } test->tds_data().mark_on_boundary(); // test is on the boundary. } *it.first++ = Facet(c, i); for (int j = 0 ; j<4 ; j++){ if (j==i) continue; if (!c->vertex(j)->get_offset_flag()) { c->vertex(j)->set_offset(int_to_off(c->offset(j))-current_off2); v_offsets.push_back(c->vertex(j)); } } } } while (!cell_stack.empty()); return it; } /*! \brief Insert point into triangulation. * * Inserts the point p into the triangulation. It expects * - a cell to start the point location * - a testing function to determine cells in conflict * - a testing function to determine if a vertex is hidden. * * Implementation: * - If the triangulation is empty call a special function * (create_initial_triangulation) to construct the basic * triangulation. * - Run point location to get the cell c containing point p. * - Call periodic_insert to insert p into the 3-cover. * - Also insert the eight periodic copies of p. */ template < class GT, class TDS > template < class Conflict_tester, class Point_hider > inline typename Periodic_3_triangulation_3::Vertex_handle Periodic_3_triangulation_3::insert_in_conflict(const Point & p, Locate_type lt, Cell_handle c, int li, int lj, const Conflict_tester &tester, Point_hider &hider) { CGAL_triangulation_assertion((_domain.xmin() <= p.x()) && (p.x() < _domain.xmax())); CGAL_triangulation_assertion((_domain.ymin() <= p.y()) && (p.y() < _domain.ymax())); CGAL_triangulation_assertion((_domain.zmin() <= p.z()) && (p.z() < _domain.zmax())); if (number_of_vertices() == 0) { return create_initial_triangulation(p); } if ((lt == VERTEX) && (tester.compare_weight(c->vertex(li)->point(),p)==0) ) { return c->vertex(li); } Vertex_handle vstart; if (!is_1_cover()) { Virtual_vertex_map_it vvmit = virtual_vertices.find(c->vertex(0)); if (vvmit == virtual_vertices.end()) vstart = c->vertex(0); else vstart = vvmit->second.first; CGAL_triangulation_assertion(virtual_vertices.find(vstart) ==virtual_vertices.end()); CGAL_triangulation_assertion(virtual_vertices_reverse.find(vstart) != virtual_vertices_reverse.end()); } CGAL_triangulation_assertion( number_of_vertices() != 0 ); CGAL_triangulation_expensive_assertion(is_valid()); Vertex_handle vh = periodic_insert(p, Offset(), lt, c, tester, hider); if (is_1_cover()) { return vh; } for (Cell_iterator it = all_cells_begin() ; it != all_cells_end() ; it++){ CGAL_triangulation_assertion(it->neighbor(0)->neighbor( it->neighbor(0)->index(it))==it); CGAL_triangulation_assertion(it->neighbor(1)->neighbor( it->neighbor(1)->index(it))==it); CGAL_triangulation_assertion(it->neighbor(2)->neighbor( it->neighbor(2)->index(it))==it); CGAL_triangulation_assertion(it->neighbor(3)->neighbor( it->neighbor(3)->index(it))==it); } std::vector start_vertices = virtual_vertices_reverse.find(vstart)->second; Cell_handle start; virtual_vertices_reverse[vh] = std::vector(); // insert 26 periodic copies for (int i=0; i<_cover[0]; i++) { for (int j=0; j<_cover[1]; j++) { for (int k=0; k<_cover[2]; k++) { if ((i!=0)||(j!=0)||(k!=0)) { start = start_vertices[i*9+j*3+k-1]->cell(); c = periodic_locate(p, Offset(i,j,k), lt, li, lj, start); periodic_insert(p, Offset(i,j,k), lt, c, tester, hider,vh); } } } } CGAL_triangulation_expensive_assertion(is_valid()); // Fall back to 1-cover if the criterion that the longest edge is shorter // than sqrt(0.166) is fulfilled. if ( too_long_edge_counter == 0 ) { CGAL_triangulation_expensive_assertion(is_valid()); convert_to_1_sheeted_covering(); CGAL_triangulation_expensive_assertion( is_valid() ); } return vh; } /// tests if two vertices of one cell are just periodic copies of each other template < class GT, class TDS > inline bool Periodic_3_triangulation_3::has_self_edges(Cell_handle c) const { CGAL_triangulation_assertion((c->vertex(0) != c->vertex(1)) || (c->offset(0) != c->offset(1))); CGAL_triangulation_assertion((c->vertex(0) != c->vertex(2)) || (c->offset(0) != c->offset(2))); CGAL_triangulation_assertion((c->vertex(0) != c->vertex(3)) || (c->offset(0) != c->offset(3))); CGAL_triangulation_assertion((c->vertex(1) != c->vertex(2)) || (c->offset(1) != c->offset(2))); CGAL_triangulation_assertion((c->vertex(1) != c->vertex(3)) || (c->offset(1) != c->offset(3))); CGAL_triangulation_assertion((c->vertex(2) != c->vertex(3)) || (c->offset(2) != c->offset(3))); return ((c->vertex(0) == c->vertex(1)) || (c->vertex(0) == c->vertex(2)) || (c->vertex(0) == c->vertex(3)) || (c->vertex(1) == c->vertex(2)) || (c->vertex(1) == c->vertex(3)) || (c->vertex(2) == c->vertex(3))); } /*! \brief Tests if the triangulation is valid. * * A triangulation is valid if * - A cell is not its own neighbor. * - A cell has no two equal neighbors * - A cell has no two equal vertex-offset pairs * - A cell is positively oriented. * - The point of a neighbor of cell c that does not belong to c is not inside * the circumcircle of c. */ template < class GT, class TDS > bool Periodic_3_triangulation_3:: is_valid(bool verbose, int level) const { bool error = false; for (Cell_iterator cit = cells_begin(); cit != cells_end(); ++cit) { for (int i=0; i<4; i++) { CGAL_triangulation_assertion(cit != cit->neighbor(i)); for (int j=i+1; j<4; j++) { CGAL_triangulation_assertion(cit->neighbor(i) != cit->neighbor(j)); CGAL_triangulation_assertion(cit->vertex(i) != cit->vertex(j)); } } // Check positive orientation: const Point *p[4]; Offset off[4]; for (int i=0; i<4; i++) { p[i] = &cit->vertex(i)->point(); off[i] = get_offset(cit,i); } if (orientation(*p[0], *p[1], *p[2], *p[3], off[0], off[1], off[2], off[3]) != POSITIVE) { if (verbose) { std::cerr<<"Periodic_3_triangulation_3: wrong orientation:"< bool Periodic_3_triangulation_3::is_valid(Cell_handle ch, bool verbose, int level) const { if ( ! _tds.is_valid(ch,verbose,level) ) return false; bool error = false; const Point *p[4]; Offset off[4]; for (int i=0; i<4; i++) { p[i] = &ch->vertex(i)->point(); off[i] = get_offset(ch,i); } if (orientation(*p[0], *p[1], *p[2], *p[3], off[0], off[1], off[2], off[3]) != POSITIVE) { error = true; } return !error; } template < class GT, class TDS > template < class ConflictTester > bool Periodic_3_triangulation_3:: is_valid_conflict(ConflictTester &tester, bool verbose, int level) const { Cell_iterator it; for ( it = cells_begin(); it != cells_end(); ++it ) { is_valid(it, verbose, level); for (int i=0; i<4; i++ ) { Offset o_nb = get_neighbor_offset(it,i,it->neighbor(i)); Offset o_vt = get_offset(it->neighbor(i), it->neighbor(i)->index(it)); if (tester(it, it->neighbor(i)->vertex(it->neighbor(i)->index(it))->point(), o_vt-o_nb)) { if (verbose) std::cerr << "non-empty sphere: " <vertex(0)->point()<<'\t' <vertex(1)->point()<<'\t' <vertex(2)->point()<<'\t' <vertex(3)->point()<<'\n' <neighbor(i)->vertex(it->neighbor(i)->index(it))->point() <<'\t'< inline void Periodic_3_triangulation_3::make_hole(Vertex_handle v, std::map &outer_map, std::vector &hole) { //CGAL_triangulation_precondition( all_vertices_begin()++ // != all_vertices_end() ); incident_cells(v, std::back_inserter(hole)); for (typename std::vector::iterator cit = hole.begin(); cit != hole.end(); ++cit) { int indv = (*cit)->index(v); Cell_handle opp_cit = (*cit)->neighbor( indv ); Facet f(opp_cit, opp_cit->index(*cit)); Vertex_triple vt = make_vertex_triple(f); make_canonical(vt); outer_map[vt] = f; for (int i=0; i<4; i++) if ( i != indv ) (*cit)->vertex(i)->set_cell(opp_cit); } } /*! \brief Remove a vertex from the triangulation. * * Removes vertex v from the triangulation. */ template < class GT, class TDS > template < class PointRemover, class Conflict_tester> inline void Periodic_3_triangulation_3::remove(Vertex_handle v, PointRemover &r, Conflict_tester &t) { CGAL_expensive_precondition(is_vertex(v)); std::vector vhrem; if (!is_1_cover()) { if (number_of_vertices() == 1) { clear(); return; } Virtual_vertex_map_it vvmit = virtual_vertices.find(v); if (vvmit != virtual_vertices.end()) v = vvmit->second.first; CGAL_triangulation_assertion(virtual_vertices_reverse.find(v) != virtual_vertices_reverse.end()); vhrem = virtual_vertices_reverse.find(v)->second; virtual_vertices_reverse.erase(v); CGAL_triangulation_assertion(vhrem.size()==26); for (int i=0 ; i<26 ; i++) { periodic_remove(vhrem[i],r); virtual_vertices.erase(vhrem[i]); CGAL_triangulation_expensive_assertion(is_valid()); } periodic_remove(v,r); } else { periodic_remove(v,r); if (!is_1_cover()) remove(v,r,t); } } /*! \brief Remove a vertex from the triangulation. * * Removes vertex v from the triangulation. * It expects a reference to an instance of a PointRemover. * * Implementation: * - Compute the hole, that is, all cells incident to v. Cells outside of * this hole are not affected by the deletion of v. * - Triangulate the hole. This is done computing the triangulation * in Euclidean space for the points on the border of the hole. * - Sew this triangulation into the hole. * - Test for all newly added edges, whether they are shorter than the * edge_length_threshold. If not, convert to 3-cover. */ template < class GT, class TDS > template < class PointRemover > inline void Periodic_3_triangulation_3::periodic_remove(Vertex_handle v, PointRemover &remover) { // Construct the set of vertex triples on the boundary // with the facet just behind typedef std::map Vertex_triple_Facet_map; typedef PointRemover Point_remover; typedef typename Point_remover::CellE_handle CellE_handle; typedef typename Point_remover::VertexE_handle VertexE_handle; typedef typename Point_remover::FacetE FacetE; typedef typename Point_remover::VertexE_triple VertexE_triple; typedef typename Point_remover::Finite_cellsE_iterator Finite_cellsE_iterator; typedef typename Point_remover::Vertex_triple_FacetE_map Vertex_triple_FacetE_map; // First compute the hole and its boundary vertices. std::vector hole; hole.reserve(64); Vertex_triple_Facet_map outer_map; Vertex_triple_FacetE_map inner_map; make_hole(v, outer_map, hole); CGAL_triangulation_assertion(outer_map.size()==hole.size()); CGAL_triangulation_assertion(remover.hidden_points_begin() == remover.hidden_points_end()); if (!is_1_cover()) { delete_too_long_edges(hole.begin(), hole.end()); } // Output the hidden points. for (typename std::vector::iterator hi = hole.begin(), hend = hole.end(); hi != hend; ++hi) { remover.add_hidden_points(*hi); } // Build up the map between Vertices on the boundary and offsets // collect all vertices on the boundary std::vector vertices; vertices.reserve(64); // The set is needed to ensure that each vertex is inserted only once. std::set tmp_vertices; // The map connects vertices to offsets in the hole std::map vh_off_map; for(typename std::vector::iterator cit = hole.begin(); cit != hole.end(); ++cit) { // Put all incident vertices in tmp_vertices. for (int j=0; j<4; ++j){ if ((*cit)->vertex(j) != v){ tmp_vertices.insert((*cit)->vertex(j)); vh_off_map[(*cit)->vertex(j)] = int_to_off((*cit)->offset(j)) - int_to_off((*cit)->offset((*cit)->index(v))); } } } // Now output the vertices. std::copy(tmp_vertices.begin(), tmp_vertices.end(), std::back_inserter(vertices)); // create a Delaunay/Regular triangulation of the points on the boundary // in Euclidean space and make a map from the vertices in remover.tmp // towards the vertices in *this Unique_hash_map vmap; CellE_handle ch; remover.tmp.clear(); for(unsigned int i=0; i < vertices.size(); i++){ typedef typename Point_remover::Triangulation_R3::Point TRPoint; CGAL_triangulation_assertion(get_offset(vertices[i]) + combine_offsets(Offset(), vh_off_map[vertices[i]]) == combine_offsets(get_offset(vertices[i]),vh_off_map[vertices[i]])); TRPoint trp = std::make_pair(vertices[i]->point(), combine_offsets( get_offset(vertices[i]), vh_off_map[vertices[i]]) ); VertexE_handle vh = remover.tmp.insert(trp, ch); vmap[vh] = vertices[i]; CGAL_triangulation_assertion(vmap.is_defined(vh)); } CGAL_triangulation_assertion(remover.tmp.number_of_vertices() != 0); // Construct the set of vertex triples of tmp // We reorient the vertex triple so that it matches those from outer_map // Also note that we use the vertices of *this, not of tmp for(Finite_cellsE_iterator it = remover.tmp.finite_cells_begin(); it != remover.tmp.finite_cells_end(); ++it){ VertexE_triple vt_aux; for(int i=0; i < 4; i++){ FacetE f = std::pair(it,i); vt_aux = VertexE_triple( f.first->vertex(vertex_triple_index(f.second,0)), f.first->vertex(vertex_triple_index(f.second,1)), f.first->vertex(vertex_triple_index(f.second,2))); if (vmap.is_defined(vt_aux.first) && vmap.is_defined(vt_aux.second) && vmap.is_defined(vt_aux.third) ) { Vertex_triple vt(vmap[vt_aux.first],vmap[vt_aux.third], vmap[vt_aux.second]); make_canonical(vt); inner_map[vt]= f; } } } // A structure for storing the new neighboring relations typedef boost::tuple Neighbor_relation; std::vector nr_vec; std::vector new_cells; // Grow inside the hole, by extending the surface while(! outer_map.empty()){ typename Vertex_triple_Facet_map::iterator oit = outer_map.begin(); typename Vertex_triple_Facet_map::value_type o_vt_f_pair = *oit; Cell_handle o_ch = o_vt_f_pair.second.first; unsigned int o_i = o_vt_f_pair.second.second; typename Vertex_triple_FacetE_map::iterator iit = inner_map.find(o_vt_f_pair.first); CGAL_triangulation_assertion(iit != inner_map.end()); typename Vertex_triple_FacetE_map::value_type i_vt_f_pair = *iit; CellE_handle i_ch = i_vt_f_pair.second.first; unsigned int i_i = i_vt_f_pair.second.second; // create a new cell to glue to the outer surface Cell_handle new_ch = _tds.create_cell(); new_cells.push_back(new_ch); new_ch->set_vertices(vmap[i_ch->vertex(0)], vmap[i_ch->vertex(1)], vmap[i_ch->vertex(2)], vmap[i_ch->vertex(3)]); set_offsets(new_ch, vh_off_map[vmap[i_ch->vertex(0)]], vh_off_map[vmap[i_ch->vertex(1)]], vh_off_map[vmap[i_ch->vertex(2)]], vh_off_map[vmap[i_ch->vertex(3)]]); // Update the edge length management for( int i=0 ; i < 4 ; i++ ) { for (int j=0 ; j < 4 ; j++) { if (j==i) continue; if (&*(new_ch->vertex(i)) > &*(new_ch->vertex(j))) continue; Point p1 = construct_point(new_ch->vertex(i)->point(), get_offset(new_ch, i)); Point p2 = construct_point(new_ch->vertex(j)->point(), get_offset(new_ch, j)); Vertex_handle v_no = new_ch->vertex(i); if (squared_distance(p1,p2) > edge_length_threshold) { // If the cell does not fulfill the edge-length criterion // revert all changes to the triangulation and transform it // to a triangulation in the needed covering space. if (is_1_cover()) { _tds.delete_cells(new_cells.begin(), new_cells.end()); convert_to_27_sheeted_covering(); return; } else if (find(too_long_edges[v_no].begin(), too_long_edges[v_no].end(), new_ch->vertex(j)) == too_long_edges[v_no].end()) { too_long_edges[v_no].push_back(new_ch->vertex(j)); too_long_edge_counter++; } } } } // The neighboring relation needs to be stored temporarily in // nr_vec. It cannot be applied directly because then we could not // easily cancel the removing process if a cell is encountered // that does not obey the edge-length criterion. nr_vec.push_back(boost::make_tuple(o_ch,o_i,new_ch)); nr_vec.push_back(boost::make_tuple(new_ch,i_i,o_ch)); // for the other faces check, if they can also be glued for(unsigned int i = 0; i < 4; i++){ if(i != i_i){ Facet f = std::pair(new_ch,i); Vertex_triple vt = make_vertex_triple(f); make_canonical(vt); std::swap(vt.second,vt.third); typename Vertex_triple_Facet_map::iterator oit2 = outer_map.find(vt); if(oit2 == outer_map.end()){ std::swap(vt.second,vt.third); outer_map[vt]= f; } else { // glue the faces typename Vertex_triple_Facet_map::value_type o_vt_f_pair2 = *oit2; Cell_handle o_ch2 = o_vt_f_pair2.second.first; int o_i2 = o_vt_f_pair2.second.second; nr_vec.push_back(boost::make_tuple(o_ch2,o_i2,new_ch)); nr_vec.push_back(boost::make_tuple(new_ch,i,o_ch2)); outer_map.erase(oit2); } } } outer_map.erase(oit); } // finally set the neighboring relations for (unsigned int i=0 ; i()->set_neighbor(nr_vec[i].template get<1>(),nr_vec[i].template get<2>()); } _tds.delete_vertex(v); _tds.delete_cells(hole.begin(), hole.end()); CGAL_triangulation_expensive_assertion(is_valid()); } // ############################################################################ // ############################################################################ // ############################################################################ template < class GT, class TDS > template < class Cmp > class Periodic_3_triangulation_3::Perturbation_order { typedef GT Geometric_traits; typedef typename Geometric_traits::Point_3 Point; typedef typename Geometric_traits::Compare_xyz_3 Compare_xyz_3; typedef typename Periodic_3_triangulation_3::Periodic_point Periodic_point; Cmp _cmp; public: Perturbation_order(const Cmp & cmp) : _cmp(cmp) {} bool operator()(const Periodic_point *p, const Periodic_point *q) const { return (_cmp(p->first, q->first, p->second, q->second) == SMALLER); } }; // ############################################################################ // ############################################################################ // ############################################################################ /** \brief Delete each redundant cell and the not anymore needed data * structures. * * This function consists of four iterations over all cells and one * iteration over all vertices: * 1. cell iteration: mark all cells that are to delete * 2. cell iteration: redirect neighbors of remaining cells * 3. cell iteration: redirect vertices of remaining cells * 4. cell iteration: delete all cells marked in the 1. iteration * Vertex iteration: delete all vertices outside the original domain. */ template < class GT, class TDS > inline void Periodic_3_triangulation_3::convert_to_1_sheeted_covering() { // ################################################################### // ### First cell iteration ########################################## // ################################################################### { if (is_1_cover()) return; bool to_delete, has_simplifiable_offset; Virtual_vertex_map_it vvmit; // First iteration over all cells: Mark the cells that are to delete. // Cells are to delete if they cannot be translated anymore in the // direction of one of the axes without yielding negative offsets. for( Cell_iterator it = all_cells_begin() ; it != all_cells_end() ; ++it ) { to_delete = false; // for all directions in 3D Space for( int j=0 ; j<3 ; j++ ) { has_simplifiable_offset = true; // for all vertices of cell it for( int i=0 ; i<4 ; i++ ) { vvmit = virtual_vertices.find(it->vertex(i)); // if it->vertex(i) lies inside the original domain: if (vvmit == virtual_vertices.end()) { // the cell cannot be moved any more because if we did, then // it->vertex(i) will get at least one negative offset. has_simplifiable_offset = false; // if it->vertex(i) lies outside the original domain: } else { // The cell can certainly be deleted if the offset contains a 2 to_delete = to_delete || (vvmit->second.second[j] == 2) ; // The cell can be moved into one direction only if the offset of // all for vertices is >=1 for this direction. Since we already // tested for 2 it is sufficient to test here for 1. has_simplifiable_offset = has_simplifiable_offset && (vvmit->second.second[j] == 1) ; } } // if the offset can be simplified, i.e. the cell can be moved, then // it can be deleted. if (has_simplifiable_offset) to_delete = true; } // Mark all cells that are to delete. They cannot be deleted yet, // because neighboring information still needs to be extracted. if (to_delete) { it->set_additional_flag(1); } } } // ################################################################### // ### Second cell iteration ######################################### // ################################################################### { Vertex_handle vert[4], nbv[4]; Offset off[4]; Cell_handle nb, new_neighbor; std::vector > new_neighbor_relations; // Second iteration over all cells: redirect neighbors where necessary for (Cell_iterator it = all_cells_begin() ; it != all_cells_end() ; ++it) { // Skip all cells that are to delete. if (it->get_additional_flag() == 1) continue; // Redirect neighbors: Only neighbors that are marked by the // additional_flag have to be substituted by one of their periodic // copies. The unmarked neighbors stay the same. for ( int i = 0 ; i < 4 ; i++ ) { if ( it->neighbor(i)->get_additional_flag() != 1 ) continue; nb = it->neighbor(i); for ( int j = 0 ; j < 4 ; j++ ) { off[j] = Offset(); get_vertex( nb, j, vert[j], off[j]); } int x,y,z; x = (std::min) ( (std::min) ( off[0][0], off[1][0] ), (std::min) ( off[2][0], off[3][0] ) ); y = (std::min) ( (std::min) ( off[0][1], off[1][1] ), (std::min) ( off[2][1], off[3][1] ) ); z = (std::min) ( (std::min) ( off[0][2], off[1][2] ), (std::min) ( off[2][2], off[3][2] ) ); // The vector from nb to the "original" periodic copy of nb, that is // the copy that will not be deleted. Offset difference_offset(x,y,z); CGAL_triangulation_assertion( !difference_offset.is_null() ); // We now have to find the "original" periodic copy of nb from // its vertices. Therefore, we first have to find the vertices. for ( int j = 0 ; j < 4 ; j++ ) { CGAL_triangulation_assertion( (off[j]-difference_offset)[0] >= 0); CGAL_triangulation_assertion( (off[j]-difference_offset)[1] >= 0); CGAL_triangulation_assertion( (off[j]-difference_offset)[2] >= 0); CGAL_triangulation_assertion( (off[j]-difference_offset)[0] < 3); CGAL_triangulation_assertion( (off[j]-difference_offset)[1] < 3); CGAL_triangulation_assertion( (off[j]-difference_offset)[2] < 3); // find the Vertex_handles of the vertices of the "original" // periodic copy of nb. If the vertex is inside the original // domain, there is nothing to do if ( (off[j]-difference_offset).is_null() ) { nbv[j] = vert[j]; // If the vertex is outside the original domain, we have to search // in virtual_vertices in the "wrong" direction. That means, we // cannot use virtual_vertices.find but have to use // virtual_vertices_reverse. } else { Offset nbo = off[j]-difference_offset; nbv[j] = virtual_vertices_reverse.find(vert[j]) ->second[nbo[0]*9+nbo[1]*3+nbo[2]-1]; } } // Find the new neighbor by its 4 vertices new_neighbor = get_cell( nbv ); // Store the new neighbor relation. This cannot be applied yet because // it would disturb the functioning of get_cell( ... ) new_neighbor_relations.push_back(make_triple(it, i, new_neighbor)); } } // Apply the new neighbor relations now. for (unsigned int i=0 ; iset_neighbor( new_neighbor_relations[i].second, new_neighbor_relations[i].third); } } // ################################################################### // ### Third cell iteration ########################################## // ################################################################### { Vertex_handle vert[4]; Offset off[4]; // Third iteration over all cells: redirect vertices where necessary for (Cell_iterator it = all_cells_begin() ; it != all_cells_end() ; ++it) { // Skip all cells that are marked to delete if (it->get_additional_flag() == 1) continue; // Find the corresponding vertices of it in the original domain // and set them as new vertices of it. for ( int i = 0 ; i < 4 ; i++ ) { off[i] = Offset(); get_vertex( it, i, vert[i], off[i]); it->set_vertex( i, vert[i]); CGAL_triangulation_assertion(vert[i]->point()[0] < _domain.xmax()); CGAL_triangulation_assertion(vert[i]->point()[1] < _domain.ymax()); CGAL_triangulation_assertion(vert[i]->point()[2] < _domain.zmax()); CGAL_triangulation_assertion(vert[i]->point()[0] >= _domain.xmin()); CGAL_triangulation_assertion(vert[i]->point()[1] >= _domain.ymin()); CGAL_triangulation_assertion(vert[i]->point()[2] >= _domain.zmin()); // redirect also the cell pointer of the vertex. it->vertex(i)->set_cell(it); } // Set the offsets. set_offsets(it, off[0], off[1], off[2], off[3] ); CGAL_triangulation_assertion( int_to_off(it->offset(0)) == off[0] ); CGAL_triangulation_assertion( int_to_off(it->offset(1)) == off[1] ); CGAL_triangulation_assertion( int_to_off(it->offset(2)) == off[2] ); CGAL_triangulation_assertion( int_to_off(it->offset(3)) == off[3] ); } } // ################################################################### // ### Fourth cell iteration ######################################### // ################################################################### { // Delete the marked cells. std::vector cells_to_delete; for ( Cell_iterator cit = all_cells_begin() ; cit != all_cells_end() ; ++cit ) { if ( cit->get_additional_flag() == 1 ) cells_to_delete.push_back( cit ); } _tds.delete_cells(cells_to_delete.begin(), cells_to_delete.end()); } // ################################################################### // ### Vertex iteration ############################################## // ################################################################### { // Delete all the vertices in virtual_vertices, that is all vertices // outside the original domain. std::vector vertices_to_delete; for ( Vertex_iterator vit = all_vertices_begin() ; vit != all_vertices_end() ; ++vit ) { if ( virtual_vertices.count( vit ) != 0 ) { CGAL_triangulation_assertion( virtual_vertices.count( vit ) == 1 ); vertices_to_delete.push_back( vit ) ; } } _tds.delete_vertices(vertices_to_delete.begin(), vertices_to_delete.end()); } _cover = make_array(1,1,1); virtual_vertices.clear(); virtual_vertices_reverse.clear(); } template < class GT, class TDS > inline void Periodic_3_triangulation_3::convert_to_27_sheeted_covering() { if (_cover == make_array(3,3,3)) return; CGAL_triangulation_precondition(is_1_cover()); // Create 27 copies of each vertex and write virtual_vertices and // virtual_vertices_reverse std::list original_vertices; // try to use std::copy instead of the following loop. for (Vertex_iterator vit = vertices_begin() ; vit != vertices_end() ; ++vit) original_vertices.push_back(vit); for (typename std::list::iterator vit = original_vertices.begin() ; vit != original_vertices.end() ; ++vit) { Vertex_handle v_cp; std::vector copies; for (int i=0; i<3; i++) for (int j=0; j<3; j++) for (int k=0; k<3; k++) { if (i==0 && j==0 && k==0) continue; v_cp = _tds.create_vertex(*vit); copies.push_back(v_cp); virtual_vertices.insert(std::make_pair(v_cp, std::make_pair(*vit,Offset(i,j,k)))); } virtual_vertices_reverse.insert(std::make_pair(*vit,copies)); } // Create 27 copies of each cell from the respective copies of the // vertices and write virtual_cells and virtual_cells_reverse. typedef std::map > Virtual_cell_map; typedef std::map > Virtual_cell_reverse_map; typedef typename Virtual_cell_reverse_map::const_iterator VCRMIT; Virtual_cell_map virtual_cells; Virtual_cell_reverse_map virtual_cells_reverse; std::list original_cells; for (Cell_iterator cit = cells_begin() ; cit != cells_end() ; ++cit) original_cells.push_back(cit); // Store vertex offsets in a separate data structure std::list< Offset > off_v; for (typename std::list::iterator vit = original_vertices.begin() ; vit != original_vertices.end() ; ++vit) { Cell_handle ccc = (*vit)->cell(); int v_index = ccc->index(*vit); off_v.push_back(int_to_off(ccc->offset(v_index))); } // Store neighboring offsets in a separate data structure std::list< array > off_nb; for (typename std::list::iterator cit = original_cells.begin() ; cit != original_cells.end() ; ++cit) { array off_nb_c; for (int i=0; i<4; i++){ Cell_handle ccc = *cit; Cell_handle nnn = ccc->neighbor(i); off_nb_c[i] = get_neighbor_offset(ccc,i,nnn); } off_nb.push_back(off_nb_c); } // Create copies of cells for (typename std::list::iterator cit = original_cells.begin() ; cit != original_cells.end() ; ++cit) { Cell_handle c_cp; Vertex_handle v0,v1,v2,v3; std::vector copies; Virtual_vertex_reverse_map_it vvrmit[4]; Offset vvoff[4]; for (int i=0; i<4; i++) { vvrmit[i] = virtual_vertices_reverse.find((*cit)->vertex(i)); CGAL_triangulation_assertion( vvrmit[i] != virtual_vertices_reverse.end()); vvoff[i] = int_to_off((*cit)->offset(i)); } Vertex_handle vvh[4]; for (int n=0; n<26; n++) { for (int i=0; i<4; i++) { // Decomposition of n into an offset (nx,ny,nz): // nx = (n+1)/9, ny = ((n+1)/3)%3, nz = (n+1)%3 int o_i = ((n+1)/9+vvoff[i].x()+3)%3; int o_j = ((n+1)/3+vvoff[i].y()+3)%3; int o_k = ((n+1)+vvoff[i].z()+3)%3; int n_c = 9*o_i+3*o_j+o_k-1; CGAL_triangulation_assertion(n_c >= -1); if (n_c == -1) vvh[i] = (*cit)->vertex(i); else vvh[i] = vvrmit[i]->second[n_c]; } c_cp = _tds.create_cell(vvh[0], vvh[1], vvh[2], vvh[3]); copies.push_back(c_cp); } virtual_cells_reverse.insert(std::make_pair(*cit,copies)); } // Set new vertices of boundary cells of the original domain. for (typename std::list::iterator cit = original_cells.begin() ; cit != original_cells.end() ; ++cit) { for (int i=0; i<4; i++) { Virtual_vertex_reverse_map_it vvrmit = virtual_vertices_reverse.find((*cit)->vertex(i)); CGAL_triangulation_assertion(vvrmit != virtual_vertices_reverse.end()); Offset vvoff = int_to_off((*cit)->offset(i)); if (!vvoff.is_null()) { int n_c = 9*vvoff.x()+3*vvoff.y()+vvoff.z()-1; CGAL_triangulation_assertion(n_c >= 0); CGAL_triangulation_assertion(static_cast(n_c) < vvrmit->second.size()); (*cit)->set_vertex(i,vvrmit->second[n_c]); } } } // Set neighboring relations of cell copies typename std::list< array >::iterator oit = off_nb.begin() ; for (typename std::list::iterator cit = original_cells.begin(); cit != original_cells.end() ; ++cit, ++oit) { CGAL_triangulation_assertion( oit != off_nb.end() ); VCRMIT c_cp = virtual_cells_reverse.find(*cit); CGAL_triangulation_assertion(c_cp != virtual_cells_reverse.end()); for (int i=0; i<4; i++) { Cell_handle cit_nb = (*cit)->neighbor(i); VCRMIT c_cp_nb = virtual_cells_reverse.find(cit_nb); CGAL_triangulation_assertion(c_cp_nb != virtual_cells_reverse.end()); Offset nboff = (*oit)[i]; for (int n=0; n<26; n++) { int n_nb; if (nboff.is_null()) n_nb = n; else { int o_i = ((n+1)/9-nboff.x()+3)%3; int o_j = ((n+1)/3-nboff.y()+3)%3; int o_k = (n+1-nboff.z()+3)%3; n_nb = 9*o_i+3*o_j+o_k-1; } if (n_nb == -1) { CGAL_triangulation_assertion(cit_nb->has_vertex( c_cp->second[n]->vertex((i+1)%4)) ); CGAL_triangulation_assertion(cit_nb->has_vertex( c_cp->second[n]->vertex((i+2)%4)) ); CGAL_triangulation_assertion(cit_nb->has_vertex( c_cp->second[n]->vertex((i+3)%4)) ); c_cp->second[n]->set_neighbor(i,cit_nb); } else { CGAL_triangulation_assertion(n_nb >= 0); CGAL_triangulation_assertion(static_cast(n_nb) <= c_cp_nb->second.size()); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex(c_cp->second[n]->vertex((i+1)%4)) ); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex(c_cp->second[n]->vertex((i+2)%4)) ); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex(c_cp->second[n]->vertex((i+3)%4)) ); c_cp->second[n]->set_neighbor(i,c_cp_nb->second[n_nb]); } } } } // Set neighboring relations of original cells oit = off_nb.begin(); for (typename std::list::iterator cit = original_cells.begin(); cit != original_cells.end() ; ++cit, ++oit) { CGAL_triangulation_assertion( oit != off_nb.end() ); for (int i=0; i<4; i++) { Offset nboff = (*oit)[i]; if (!nboff.is_null()) { Cell_handle cit_nb = (*cit)->neighbor(i); VCRMIT c_cp_nb = virtual_cells_reverse.find(cit_nb); CGAL_triangulation_assertion(c_cp_nb != virtual_cells_reverse.end()); int o_i = (3-nboff.x())%3; int o_j = (3-nboff.y())%3; int o_k = (3-nboff.z())%3; int n_nb = 9*o_i+3*o_j+o_k-1; CGAL_triangulation_assertion(n_nb >= 0); CGAL_triangulation_assertion(static_cast(n_nb) <= c_cp_nb->second.size()); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex((*cit)->vertex((i+1)%4)) ); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex((*cit)->vertex((i+2)%4)) ); CGAL_triangulation_assertion(c_cp_nb->second[n_nb] ->has_vertex((*cit)->vertex((i+3)%4)) ); (*cit)->set_neighbor(i,c_cp_nb->second[n_nb]); } } } // Set incident cells for (Cell_iterator cit = cells_begin() ; cit != cells_end() ; ++cit) { for (int i=0 ; i<4 ; i++) { cit->vertex(i)->set_cell(cit); } } // Set offsets where necessary for (typename std::list::iterator cit = original_cells.begin() ; cit != original_cells.end() ; ++cit) { VCRMIT c_cp = virtual_cells_reverse.find(*cit); CGAL_triangulation_assertion( c_cp != virtual_cells_reverse.end()); Offset off[4]; for (int i=0; i<4; i++) off[i] = int_to_off((*cit)->offset(i)); if (off[0].is_null() && off[1].is_null() && off[2].is_null() && off[3].is_null()) continue; for (int n=0; n<26; n++) { Offset off_cp[4]; int o_i = (n+1)/9; int o_j = ((n+1)/3)%3; int o_k = (n+1)%3; if (o_i!=2 && o_j!=2 && o_k !=2) continue; for (int i=0; i<4; i++) { off_cp[i] = Offset((o_i==2)?off[i].x():0, (o_j==2)?off[i].y():0, (o_k==2)?off[i].z():0); CGAL_triangulation_assertion(off_cp[i].x() == 0 || off_cp[i].x() == 1); CGAL_triangulation_assertion(off_cp[i].y() == 0 || off_cp[i].y() == 1); CGAL_triangulation_assertion(off_cp[i].z() == 0 || off_cp[i].z() == 1); } set_offsets(c_cp->second[n],off_cp[0],off_cp[1],off_cp[2],off_cp[3]); } } // Iterate over all original cells and reset offsets. for (typename std::list::iterator cit = original_cells.begin() ; cit != original_cells.end() ; ++cit) { //This statement does not seem to have any effect set_offsets(*cit, 0,0,0,0); CGAL_triangulation_assertion((*cit)->offset(0) == 0); CGAL_triangulation_assertion((*cit)->offset(1) == 0); CGAL_triangulation_assertion((*cit)->offset(2) == 0); CGAL_triangulation_assertion((*cit)->offset(3) == 0); } _cover = make_array(3,3,3); CGAL_triangulation_expensive_assertion(is_valid()); // Set up too long edges data structure int i=0; for (Vertex_iterator vit = vertices_begin(); vit != vertices_end(); ++vit) { too_long_edges[vit] = std::list(); ++i; } too_long_edge_counter = find_too_long_edges(too_long_edges); } // iterate over all edges and store the ones that are longer than // edge_length_threshold in edges. Return the number of too long edges. template < class GT, class TDS > inline int Periodic_3_triangulation_3::find_too_long_edges( std::map >& edges) const { Point p1, p2; int counter = 0; Vertex_handle v_no,vh; for (Edge_iterator eit = edges_begin(); eit != edges_end() ; eit++) { p1 = construct_point(eit->first->vertex(eit->second)->point(), get_offset(eit->first, eit->second)); p2 = construct_point(eit->first->vertex(eit->third)->point(), get_offset(eit->first, eit->third)); if (squared_distance(p1,p2) > edge_length_threshold) { if (&*(eit->first->vertex(eit->second)) < &*(eit->first->vertex(eit->third))) { v_no = eit->first->vertex(eit->second); vh = eit->first->vertex(eit->third); } else { v_no = eit->first->vertex(eit->third); vh = eit->first->vertex(eit->second); } edges[v_no].push_back(vh); counter++; } } return counter; } template < class GT, class TDS > class Periodic_3_triangulation_3::Finder { const Self* _t; const Point & _p; public: Finder(const Self* t, const Point &p) : _t(t), _p(p) {} bool operator()(const Vertex_handle v) { return _t->equal(v->point(), _p); } }; /** Find the cell that consists of the four given vertices * * Iterates over all cells and compare the four vertices of each cell * with the four vertices in vh. */ template < class GT, class TDS > inline typename Periodic_3_triangulation_3::Cell_handle Periodic_3_triangulation_3::get_cell(const Vertex_handle* vh) const { bool contains_v[4]; std::vector cells; incident_cells(vh[3],std::back_inserter(cells)); for ( typename std::vector::iterator it = cells.begin(); it != cells.end(); it++ ) { CGAL_triangulation_assertion( (*it)->vertex(0) == vh[3] || (*it)->vertex(1) == vh[3] ||(*it)->vertex(2) == vh[3] || (*it)->vertex(3) == vh[3]) ; for ( int j=0 ; j<3 ; j++ ) { contains_v[j] = false; contains_v[j] = ( (*it)->vertex(0) == vh[j] ) || ( (*it)->vertex(1) == vh[j] ) || ( (*it)->vertex(2) == vh[j] ) || ( (*it)->vertex(3) == vh[j] ); } if (contains_v[0] && contains_v[1] && contains_v[2]) { return (*it); } } CGAL_triangulation_assertion(false); return Cell_handle(); } /*! \brief Get the offset of tester.point() such that * this point is in conflict with c w.r.t tester.get_offset(). * * Implementation: Just try all eight possibilities. */ template < class GT, class TDS > template < class Conflict_tester > inline typename Periodic_3_triangulation_3::Offset Periodic_3_triangulation_3::get_location_offset( const Conflict_tester& tester, Cell_handle c) const { CGAL_triangulation_precondition( number_of_vertices() != 0 ); // CGAL_triangulation_precondition_code(Locate_type lt; int i; int j;); // CGAL_triangulation_precondition(side_of_cell(q,o,c,lt,i,j) // != ON_UNBOUNDED_SIDE); int cumm_off = c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3); if (cumm_off == 0) { // default case: return Offset(); } else { // Main idea seems to just test all possibilities. for (int i=0; i<8; i++) { if (((cumm_off | (~i))&7) == 7) { if (tester(c,int_to_off(i))) { return int_to_off(i); } } } } CGAL_triangulation_assertion(false); return Offset(); } /** Get the offset between the origins of the internal offset coordinate * systems of two neighboring cells with respect from ch to nb. * * - Find two corresponding vertices from each cell * - Return the difference of their offsets. */ template < class GT, class TDS > inline typename Periodic_3_triangulation_3::Offset Periodic_3_triangulation_3::get_neighbor_offset( Cell_handle ch, int i, Cell_handle nb) const { // Redundance in the signature! CGAL_triangulation_precondition(ch->neighbor(i) == nb); CGAL_triangulation_precondition(nb->neighbor(nb->index(ch)) == ch); Vertex_handle vertex_ch; int index_ch, index_nb; // ensure that vertex_ch \in nb and vertex_nb \in ch index_ch = (i==0? 1 : 0); vertex_ch = ch->vertex(index_ch); index_nb = nb->index(vertex_ch); return int_to_off(nb->offset(index_nb)) - int_to_off(ch->offset(index_ch)); } /** * - ch->offset(i) is an bit triple encapsulated in an integer. Each bit * represents the offset in one direction --> 2-cover! * - it_to_off(int) decodes this again. * - Finally the offset vector is multiplied by cover. * So if we are working in 3-cover we translate it to the neighboring * 3-cover and not only to the neighboring domain. */ template < class GT, class TDS > inline void Periodic_3_triangulation_3::get_vertex( Cell_handle ch, int i, Vertex_handle &vh, Offset &off) const { off = combine_offsets(Offset(),int_to_off(ch->offset(i))); vh = ch->vertex(i); if (is_1_cover()) return; Vertex_handle vh_i = vh; get_vertex(vh_i, vh, off); return; } template < class GT, class TDS > inline void Periodic_3_triangulation_3::get_vertex( Vertex_handle vh_i, Vertex_handle &vh, Offset &off) const { Virtual_vertex_map_it it = virtual_vertices.find(vh_i); if (it == virtual_vertices.end()) { // if ch->vertex(i) is not contained in virtual_vertices, then it is in // the original domain. vh = vh_i; CGAL_triangulation_assertion(vh != Vertex_handle()); } else { // otherwise it has to be looked up as well as its offset. vh = it->second.first; off += it->second.second; CGAL_triangulation_assertion(vh->point().x() < _domain.xmax()); CGAL_triangulation_assertion(vh->point().y() < _domain.ymax()); CGAL_triangulation_assertion(vh->point().z() < _domain.zmax()); CGAL_triangulation_assertion(vh->point().x() >= _domain.xmin()); CGAL_triangulation_assertion(vh->point().y() >= _domain.ymin()); CGAL_triangulation_assertion(vh->point().z() >= _domain.zmin()); } } template < class GT, class TDS > std::istream & operator>> (std::istream& is, Periodic_3_triangulation_3 &tr) // reads // the current covering that guarantees the triangulation to be a // simplicial complex // the number of vertices // the non combinatorial information on vertices (points in case of 1-sheeted // covering, point-offset pairs otherwise) // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY // the number of cells // the cells by the indices of their vertices in the preceding list // of vertices, plus the non combinatorial information on each cell // the neighbors of each cell by their index in the preceding list of cells { CGAL_triangulation_precondition(is.good()); typedef Periodic_3_triangulation_3 Triangulation; typedef typename GT::FT FT; typedef typename Triangulation::size_type size_type; typedef typename Triangulation::Vertex_handle Vertex_handle; typedef typename Triangulation::Cell_handle Cell_handle; typedef typename Triangulation::Offset Offset; typedef typename Triangulation::Iso_cuboid Iso_cuboid; tr.clear(); Iso_cuboid domain(0,0,0,1,1,1); int cx=0, cy=0, cz=0; size_type n=0; if (is_ascii(is)) { is >> domain; is >> cx >> cy >> cz; is >> n; } else { is >> domain; read(is,cx); read(is,cy); read(is,cz); read(is,n); } CGAL_triangulation_assertion((n/(cx*cy*cz))*cx*cy*cz == n); tr.tds().set_dimension((n==0?-2:3)); tr._domain = domain; tr._gt.set_domain(domain); tr._cover = make_array(cx,cy,cz); if ( n==0 ) return is; std::vector< Vertex_handle > V(n); if (cx==1 && cy==1 && cz==1) { for (std::size_t i=0; i < n; i++) { V[i] = tr.tds().create_vertex(); is >> *V[i]; } } else { Vertex_handle v,w; std::vector vv; Offset off; for (std::size_t i=0; i < n; i++) { v = tr.tds().create_vertex(); V[i] = v; is >> *V[i] >> off; vv.clear(); for (int j=1; j> *V[i] >> off; vv.push_back(w); tr.virtual_vertices[w]=std::make_pair(v,off); } tr.virtual_vertices_reverse[v]=vv; } } std::vector< Cell_handle > C; std::size_t m; tr._tds.read_cells(is, V, m, C); // read offsets int off[4] = {0,0,0,0}; for (std::size_t j=0 ; j < m; j++) { if (is_ascii(is)) is >> off[0] >> off[1] >> off[2] >> off[3]; else { read(is,off[0]); read(is,off[1]); read(is,off[2]); read(is,off[3]); } tr.set_offsets(C[j],off[0],off[1],off[2],off[3]); } // read potential other information for (std::size_t j=0 ; j < m; j++) is >> *(C[j]); typedef typename Triangulation::Vertex_iterator VI; int i=0; for (VI vi = tr.vertices_begin(); vi != tr.vertices_end(); ++vi) { tr.too_long_edges[vi]=std::list(); ++i; } tr.edge_length_threshold = FT(0.166) * (tr._domain.xmax()-tr._domain.xmin()) * (tr._domain.xmax()-tr._domain.xmin()); tr.too_long_edge_counter = tr.find_too_long_edges(tr.too_long_edges); CGAL_triangulation_expensive_assertion( tr.is_valid() ); return is; } template < class GT, class TDS > std::ostream & operator<< (std::ostream& os,const Periodic_3_triangulation_3 &tr) // writes : // the number of vertices // the domain as six coordinates: xmin ymin zmin xmax ymax zmax // the current covering that guarantees the triangulation to be a // simplicial complex // the non combinatorial information on vertices (points in case of 1-sheeted // covering, point-offset pairs otherwise) // ALL PERIODIC COPIES OF ONE VERTEX MUST BE STORED CONSECUTIVELY // the number of cells // the cells by the indices of their vertices in the preceding list // of vertices, plus the non combinatorial information on each cell // the neighbors of each cell by their index in the preceding list of cells { typedef Periodic_3_triangulation_3 Triangulation; typedef typename Triangulation::size_type size_type; typedef typename Triangulation::Vertex_handle Vertex_handle; typedef typename Triangulation::Vertex_iterator Vertex_iterator; typedef typename Triangulation::Cell_handle Cell_handle; typedef typename Triangulation::Cell_iterator Cell_iterator; typedef typename Triangulation::Covering_sheets Covering_sheets; typedef typename Triangulation::Offset Offset; typedef typename Triangulation::Virtual_vertex_map_it Virtual_vertex_map_it; typedef typename Triangulation::Iso_cuboid Iso_cuboid; // outputs dimension, domain and number of vertices Iso_cuboid domain = tr.domain(); Covering_sheets cover = tr.number_of_sheets(); size_type n = tr.number_of_vertices(); if (is_ascii(os)) os << domain << std::endl << cover[0] << " " << cover[1] << " " << cover[2] << std::endl << n*cover[0]*cover[1]*cover[2] << std::endl; else { os << domain; write(os,cover[0]); write(os,cover[1]); write(os,cover[2]); write(os,n*cover[0]*cover[1]*cover[2]); } if (n == 0) return os; // write the vertices Unique_hash_map V; std::size_t i=0; if (tr.is_1_cover()) { for (Vertex_iterator it=tr.vertices_begin(); it!=tr.vertices_end(); ++it) { V[it] = i++; os << it->point(); if (is_ascii(os)) os << std::endl; } } else { Virtual_vertex_map_it vit, vvit; std::vector vv; for (Vertex_iterator it=tr.vertices_begin(); it!=tr.vertices_end(); ++it) { vit = tr.virtual_vertices.find(it); if (vit != tr.virtual_vertices.end()) continue; V[it]=i++; if (is_ascii(os)) os << it->point() << std::endl << Offset(0,0,0) << std::endl; else os << it->point() << Offset(0,0,0); CGAL_triangulation_assertion(tr.virtual_vertices_reverse.find(it) != tr.virtual_vertices_reverse.end()); vv = tr.virtual_vertices_reverse.find(it)->second; CGAL_triangulation_assertion(vv.size() == 26); for (std::size_t j=0; jpoint() << std::endl << vvit->second.second << std::endl; else os << vv[j]->point() << vvit->second.second; } } } CGAL_triangulation_postcondition(i==tr._cover[0]*tr._cover[1]*tr._cover[2]*n); // asks the tds for the combinatorial information tr.tds().print_cells(os, V); // write offsets //for (unsigned int i=0 ; ioffset(j); if ( j==3 ) os << std::endl; else os << ' '; } else write(os,ch->offset(j)); } } // write the non combinatorial information on the cells // using the << operator of Cell // works because the iterator of the tds traverses the cells in the // same order as the iterator of the triangulation if(tr.number_of_vertices() != 0) { for(Cell_iterator it=tr.cells_begin(); it != tr.cells_end(); ++it) { os << *it; // other information if(is_ascii(os)) os << std::endl; } } return os ; } namespace internal { /// Internal function used by operator==. // This function tests and registers the 4 neighbors of c1/c2, // and performs a bfs traversal // Returns false if an inequality has been found. //TODO: introduce offsets template bool test_next(const Periodic_3_triangulation_3 &t1, const Periodic_3_triangulation_3 & /* needed_for_deducing_TDS2 */, typename Periodic_3_triangulation_3::Cell_handle c1, typename Periodic_3_triangulation_3::Cell_handle c2, std::map::Cell_handle, typename Periodic_3_triangulation_3::Cell_handle> &Cmap, std::map::Vertex_handle, typename Periodic_3_triangulation_3::Vertex_handle> &Vmap) { typedef Periodic_3_triangulation_3 Tr1; typedef Periodic_3_triangulation_3 Tr2; typedef typename Tr1::Vertex_handle Vertex_handle1; typedef typename Tr1::Cell_handle Cell_handle1; typedef typename Tr2::Vertex_handle Vertex_handle2; typedef typename Tr2::Cell_handle Cell_handle2; typedef typename std::map::const_iterator Cit; typedef typename std::map::const_iterator Vit; std::vector > queue; queue.push_back(std::make_pair(c1,c2)); while(! queue.empty()){ boost::tie(c1,c2) = queue.back(); queue.pop_back(); // Precondition: c1, c2 have been registered as well as their 4 vertices. CGAL_triangulation_precondition(t1.number_of_vertices() != 0); CGAL_triangulation_precondition(Cmap[c1] == c2); CGAL_triangulation_precondition(Vmap.find(c1->vertex(0)) != Vmap.end()); CGAL_triangulation_precondition(Vmap.find(c1->vertex(1)) != Vmap.end()); CGAL_triangulation_precondition(Vmap.find(c1->vertex(2)) != Vmap.end()); CGAL_triangulation_precondition(Vmap.find(c1->vertex(3)) != Vmap.end()); for (int i=0; i <= 3; ++i) { Cell_handle1 n1 = c1->neighbor(i); Cit cit = Cmap.find(n1); Vertex_handle1 v1 = c1->vertex(i); Vertex_handle2 v2 = Vmap[v1]; Cell_handle2 n2 = c2->neighbor(c2->index(v2)); if (cit != Cmap.end()) { // n1 was already registered. if (cit->second != n2) return false; continue; } // n1 has not yet been registered. // We check that the new vertices match geometrically. // And we register them. Vertex_handle1 vn1 = n1->vertex(n1->index(c1)); Vertex_handle2 vn2 = n2->vertex(n2->index(c2)); Vit vit = Vmap.find(vn1); if (vit != Vmap.end()) { // vn1 already registered if (vit->second != vn2) return false; } else { if (t1.geom_traits().compare_xyz_3_object()(vn1->point(), vn2->point()) != 0) return false; // We register vn1/vn2. Vmap.insert(std::make_pair(vn1, vn2)); } // We register n1/n2. Cmap.insert(std::make_pair(n1, n2)); queue.push_back(std::make_pair(n1, n2)); } } return true; } } // namespace internal template < class GT, class TDS1, class TDS2 > bool operator==(const Periodic_3_triangulation_3 &t1, const Periodic_3_triangulation_3 &t2) { typedef typename Periodic_3_triangulation_3::Vertex_handle Vertex_handle1; typedef typename Periodic_3_triangulation_3::Cell_handle Cell_handle1; typedef typename Periodic_3_triangulation_3::Vertex_handle Vertex_handle2; typedef typename Periodic_3_triangulation_3::Vertex_handle Vertex_iterator2; typedef typename Periodic_3_triangulation_3::Cell_handle Cell_handle2; typedef typename Periodic_3_triangulation_3::Point Point; typedef typename Periodic_3_triangulation_3::Offset Offset; // typedef typename Periodic_3_triangulation_3 // ::Geometric_traits::Compare_xyz_3 Compare_xyz_3; // Compare_xyz_3 cmp1 = t1.geom_traits().compare_xyz_3_object(); // Compare_xyz_3 cmp2 = t2.geom_traits().compare_xyz_3_object(); // Some quick checks. if ( t1.domain() != t2.domain() || t1.number_of_sheets() != t2.number_of_sheets()) return false; if ( t1.number_of_vertices() != t2.number_of_vertices() || t1.number_of_cells() != t2.number_of_cells()) return false; // Special case for empty triangulations if (t1.number_of_vertices() == 0) return true; // We will store the mapping between the 2 triangulations vertices and // cells in 2 maps. std::map Vmap; std::map Cmap; // find a common point Vertex_handle1 v1 = static_cast(t1.vertices_begin()); Vertex_handle2 iv2; for (Vertex_iterator2 vit2 = t2.vertices_begin() ; vit2 != t2.vertices_end(); ++vit2) { if (!t1.equal(vit2->point(), v1->point(), t2.get_offset(vit2), t1.get_offset(v1))) continue; iv2 = static_cast(vit2); break; } if (iv2 == Vertex_handle2()) return false; Vmap.insert(std::make_pair(v1, iv2)); // We pick one cell of t1, and try to match it against the // cells of t2. Cell_handle1 c = v1->cell(); Vertex_handle1 v2 = c->vertex((c->index(v1)+1)%4); Vertex_handle1 v3 = c->vertex((c->index(v1)+2)%4); Vertex_handle1 v4 = c->vertex((c->index(v1)+3)%4); Point p2 = v2->point(); Point p3 = v3->point(); Point p4 = v4->point(); Offset o2 = t1.get_offset(v2); Offset o3 = t1.get_offset(v3); Offset o4 = t1.get_offset(v4); std::vector ics; t2.incident_cells(iv2, std::back_inserter(ics)); for (typename std::vector::const_iterator cit = ics.begin(); cit != ics.end(); ++cit) { int inf = (*cit)->index(iv2); if (t1.equal(p2, (*cit)->vertex((inf+1)%4)->point(), o2, t2.get_offset((*cit)->vertex((inf+1)%4)))) Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+1)%4))); else if (t1.equal(p2, (*cit)->vertex((inf+2)%4)->point(), o2, t2.get_offset((*cit)->vertex((inf+2)%4)))) Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+2)%4))); else if (t1.equal(p2, (*cit)->vertex((inf+3)%4)->point(), o2, t2.get_offset((*cit)->vertex((inf+3)%4)))) Vmap.insert(std::make_pair(v2, (*cit)->vertex((inf+3)%4))); else continue; // None matched v2. if (t1.equal(p3, (*cit)->vertex((inf+1)%4)->point(), o3, t2.get_offset((*cit)->vertex((inf+1)%4)))) Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+1)%4))); else if (t1.equal(p3, (*cit)->vertex((inf+2)%4)->point(), o3, t2.get_offset((*cit)->vertex((inf+2)%4)))) Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+2)%4))); else if (t1.equal(p3, (*cit)->vertex((inf+3)%4)->point(), o3, t2.get_offset((*cit)->vertex((inf+3)%4)))) Vmap.insert(std::make_pair(v3, (*cit)->vertex((inf+3)%4))); else continue; // None matched v3. if (t1.equal(p4, (*cit)->vertex((inf+1)%4)->point(), o4, t2.get_offset((*cit)->vertex((inf+1)%4)))) Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+1)%4))); else if (t1.equal(p4, (*cit)->vertex((inf+2)%4)->point(), o4, t2.get_offset((*cit)->vertex((inf+2)%4)))) Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+2)%4))); else if (t1.equal(p4, (*cit)->vertex((inf+3)%4)->point(), o4, t2.get_offset((*cit)->vertex((inf+3)%4)))) Vmap.insert(std::make_pair(v4,(*cit)->vertex((inf+3)%4))); else continue; // None matched v4. // Found it ! Cmap.insert(std::make_pair(c, *cit)); break; } if (Cmap.size() == 0) return false; // We now have one cell, we need to compare in a bfs graph traversal return internal::test_next(t1, t2, Cmap.begin()->first, Cmap.begin()->second, Cmap, Vmap); } template < class GT, class TDS1, class TDS2 > inline bool operator!=(const Periodic_3_triangulation_3 &t1, const Periodic_3_triangulation_3 &t2) { return ! (t1 == t2); } } //namespace CGAL #endif // CGAL_PERIODIC_3_TRIANGULATION_3_H