// Copyright (c) 1997-2013 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Periodic_2_triangulation_2/include/CGAL/Periodic_2_triangulation_hierarchy_2.h $ // $Id: Periodic_2_triangulation_hierarchy_2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Olivier Devillers // Mariette Yvinec // Nico Kruithof #ifndef CGAL_PERIODIC_2_TRIANGULATION_HIERARCHY_2_H #define CGAL_PERIODIC_2_TRIANGULATION_HIERARCHY_2_H #include #include #include #include #include #include #include #include #include namespace CGAL { template < class PTr> class Periodic_2_triangulation_hierarchy_2 : public PTr { // parameterization of the hierarchy static const int m_ratio = 30; static const int m_minsize = 20; static const int m_maxlevel = 5; // maximal number of points is 30^5 = 24 millions ! public: typedef PTr PTr_Base; typedef typename PTr::Geom_traits Geom_traits; typedef typename PTr::Point Point; typedef typename PTr::Iso_rectangle Iso_rectangle; typedef typename PTr::size_type size_type; typedef typename PTr::Vertex_handle Vertex_handle; typedef typename PTr::Face_handle Face_handle; typedef typename PTr::Vertex Vertex; typedef typename PTr::Locate_type Locate_type; typedef typename PTr::Finite_vertices_iterator Finite_vertices_iterator; //typedef typename PTr::Finite_faces_iterator Finite_faces_iterator; typedef typename PTr::Weighted_tag Weighted_tag; typedef typename PTr::Periodic_tag Periodic_tag; #ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 using PTr_Base::geom_traits; #endif private: // here is the stack of triangulations which form the hierarchy PTr_Base* hierarchy[m_maxlevel]; boost::rand48 random; int level_mult_cover; public: Periodic_2_triangulation_hierarchy_2( const Iso_rectangle& domain = Iso_rectangle(0, 0, 1, 1), const Geom_traits& traits = Geom_traits()); Periodic_2_triangulation_hierarchy_2( const Periodic_2_triangulation_hierarchy_2& tr); template < typename InputIterator > Periodic_2_triangulation_hierarchy_2(InputIterator first, InputIterator last, const Iso_rectangle& domain = Iso_rectangle(0, 0, 1, 1), const Geom_traits& traits = Geom_traits()) : PTr_Base(domain, traits), level_mult_cover(0) { hierarchy[0] = this; for(int i = 1; i < m_maxlevel; ++i) hierarchy[i] = new PTr_Base(domain, traits); insert(first, last); } Periodic_2_triangulation_hierarchy_2 &operator=(const Periodic_2_triangulation_hierarchy_2& tr); ~Periodic_2_triangulation_hierarchy_2(); //Helping void copy_triangulation(const Periodic_2_triangulation_hierarchy_2 &tr); void swap(Periodic_2_triangulation_hierarchy_2 &tr); void clear(); // CHECKING bool is_valid(bool verbose = false, int level = 0) const; // INSERT REMOVE MOVE Vertex_handle insert(const Point &p, Face_handle start = Face_handle() ); Vertex_handle insert(const Point& p, Locate_type lt, Face_handle loc, int li ); Vertex_handle push_back(const Point &p); template < class InputIterator > std::ptrdiff_t insert(InputIterator first, InputIterator last, bool /* is_large_point_set */ = false) { std::ptrdiff_t n = this->number_of_vertices(); std::vector points (first, last); CGAL::spatial_sort (points.begin(), points.end(), geom_traits()); // hints[i] is the face of the previously inserted point in level i. // Thanks to spatial sort, they are better hints than what the hierarchy // would give us. Face_handle hints[m_maxlevel]; for (typename std::vector::const_iterator p = points.begin(), end = points.end(); p != end; ++p) { int vertex_level = random_level(); Vertex_handle v = hierarchy[0]->insert (*p, hints[0]); hints[0] = v->face(); Vertex_handle prev = v; for (int level = 1; level <= vertex_level; ++level) { v = hierarchy[level]->insert (*p, hints[level]); hints[level] = v->face(); v->set_down (prev); if (hierarchy[level]->number_of_sheets()[0] != 1) { std::vector vtc = hierarchy[level]->periodic_copies(v); for (unsigned int i = 0 ; i < vtc.size() ; i++) vtc[i]->set_down(prev); } prev->set_up (v); prev = v; } } std::ptrdiff_t m = this->number_of_vertices(); return m - n; } void remove(Vertex_handle v); void remove_first(Vertex_handle v); void remove_degree_3(Vertex_handle v); Vertex_handle move_point(Vertex_handle v, const Point &p); Vertex_handle move_if_no_collision(Vertex_handle v, const Point &p); //LOCATE Face_handle locate(const Point& p, Locate_type& lt, int& li, Face_handle start = Face_handle()) const; Face_handle locate(const Point &p, Face_handle start = Face_handle()) const; Vertex_handle nearest_vertex(const Point& p, Face_handle start = Face_handle()) const { return PTr_Base::nearest_vertex(p, start != Face_handle() ? start : locate(p)); } private: void locate_in_all(const Point& p, Locate_type& lt, int& li, Face_handle loc, Face_handle pos[m_maxlevel]) const; int random_level(); }; template Periodic_2_triangulation_hierarchy_2:: Periodic_2_triangulation_hierarchy_2(const Iso_rectangle& domain, const Geom_traits& traits) : PTr_Base(domain, traits) { level_mult_cover = 0; hierarchy[0] = this; for(int i = 1; i < m_maxlevel; ++i) hierarchy[i] = new PTr_Base(domain, traits); } // copy constructor duplicates vertices and faces template Periodic_2_triangulation_hierarchy_2:: Periodic_2_triangulation_hierarchy_2(const Periodic_2_triangulation_hierarchy_2 &tr) : PTr_Base() { // create an empty triangulation to be able to delete it ! hierarchy[0] = this; for(int i = 1; i < m_maxlevel; ++i) hierarchy[i] = new PTr_Base(tr.domain(), tr.geom_traits()); copy_triangulation(tr); } //Assignement template Periodic_2_triangulation_hierarchy_2 & Periodic_2_triangulation_hierarchy_2:: operator=(const Periodic_2_triangulation_hierarchy_2 &tr) { Periodic_2_triangulation_hierarchy_2 tmp(tr); swap(tmp); return *this; } template void Periodic_2_triangulation_hierarchy_2:: copy_triangulation(const Periodic_2_triangulation_hierarchy_2 &tr) { { for(int i = 0; i < m_maxlevel; ++i) hierarchy[i]->copy_triangulation(*tr.hierarchy[i]); } //up and down have been copied in straightforward way // compute a map at lower level std::map V; { for(Finite_vertices_iterator it = hierarchy[0]->finite_vertices_begin(); it != hierarchy[0]->finite_vertices_end(); ++it) { if (it->up() != Vertex_handle()) V[ it->up()->down() ] = it; } } { for(int i = 1; i < m_maxlevel; ++i) { for( Finite_vertices_iterator it = hierarchy[i]->finite_vertices_begin(); it != hierarchy[i]->finite_vertices_end(); ++it) { if (hierarchy[i]->is_virtual(it)) { // down pointer goes in original instead in copied triangulation it->set_down(V[it->down()]); // make reverse link it->down()->set_up(it); // I think the next line is unnecessary (my) // make map for next level if (it->up() != Vertex_handle() ) V[ it->up()->down() ] = it; } } } } } template void Periodic_2_triangulation_hierarchy_2:: swap(Periodic_2_triangulation_hierarchy_2 &tr) { PTr_Base::swap(tr); for(int i = 1; i < m_maxlevel; ++i) std::swap(hierarchy[i], tr.hierarchy[i]); } template Periodic_2_triangulation_hierarchy_2:: ~Periodic_2_triangulation_hierarchy_2() { clear(); for(int i = 1; i < m_maxlevel; ++i) delete hierarchy[i]; } template void Periodic_2_triangulation_hierarchy_2:: clear() { for(int i = 0; i < m_maxlevel; ++i) { CGAL_assertion(hierarchy[i] != nullptr); hierarchy[i]->clear(); } } template bool Periodic_2_triangulation_hierarchy_2:: is_valid(bool verbose, int level) const { bool result = true; int i; Finite_vertices_iterator it; //verify correctness of triangulation at all levels for(i = 0; i < m_maxlevel; ++i) { if(verbose) // print number of vertices at each level std::cout << "number_of_vertices " << hierarchy[i]->number_of_vertices() << std::endl; result = result && hierarchy[i]->is_valid(verbose, level); } //verify that lower level has no down pointers for( it = hierarchy[0]->finite_vertices_begin(); it != hierarchy[0]->finite_vertices_end(); ++it) if (!hierarchy[0]->is_virtual(it)) result = result && (it->down() == Vertex_handle()); //verify that other levels have down pointer and reciprocal link is fine for(i = 1; i < m_maxlevel; ++i) for( it = hierarchy[i]->finite_vertices_begin(); it != hierarchy[i]->finite_vertices_end(); ++it) if (!hierarchy[i]->is_virtual(it)) result = result && (&*(it->down()->up()) == &*(it)); //verify that levels have up pointer and reciprocal link is fine for(i = 0; i < m_maxlevel - 1; ++i) for( it = hierarchy[i]->finite_vertices_begin(); it != hierarchy[i]->finite_vertices_end(); ++it) if (!hierarchy[i]->is_virtual(it)) result = result && ( it->up() == Vertex_handle() || &*it == &*(it->up())->down() ); return result; } template typename Periodic_2_triangulation_hierarchy_2::Vertex_handle Periodic_2_triangulation_hierarchy_2:: insert(const Point &p, Face_handle loc) { int vertex_level = random_level(); Locate_type lt; int i; // locate using hierarchy Face_handle positions[m_maxlevel]; locate_in_all(p, lt, i, loc, positions); Vertex_handle vertex = hierarchy[0]->PTr_Base::insert(p, lt, positions[0], i); Vertex_handle previous = vertex; Vertex_handle first = vertex; int level = 1; while (level <= vertex_level ) { vertex = hierarchy[level]->PTr_Base::insert(p, positions[level]); vertex->set_down(previous);// link with level above if (hierarchy[level]->number_of_sheets()[0] != 1) { std::vector vtc = hierarchy[level]->periodic_copies(vertex); for (unsigned int i = 0 ; i < vtc.size() ; i++) vtc[i]->set_down(previous); } previous->set_up(vertex); previous = vertex; level++; } return first; } template typename Periodic_2_triangulation_hierarchy_2::Vertex_handle Periodic_2_triangulation_hierarchy_2:: insert(const Point& p, Locate_type lt, Face_handle loc, int li ) { int vertex_level = random_level(); //insert at level 0 Vertex_handle vertex = hierarchy[0]->PTr_Base::insert(p, lt, loc, li); Vertex_handle previous = vertex; Vertex_handle first = vertex; if (vertex_level > 0) { // locate using hierarchy Locate_type ltt; int lii; Face_handle positions[m_maxlevel]; locate_in_all(p, ltt, lii, loc, positions); //insert in higher levels int level = 1; while (level <= vertex_level ) { vertex = hierarchy[level]->PTr_Base::insert(p, positions[level]); vertex->set_down(previous);// link with level above if (hierarchy[level]->number_of_sheets()[0] != 1) { std::vector vtc = hierarchy[level]->periodic_copies(vertex); for (unsigned int i = 0 ; i < vtc.size() ; i++) vtc[i]->set_down(previous); } previous->set_up(vertex); previous = vertex; level++; } } return first; } template inline typename Periodic_2_triangulation_hierarchy_2::Vertex_handle Periodic_2_triangulation_hierarchy_2:: push_back(const Point &p) { return insert(p); } template void Periodic_2_triangulation_hierarchy_2:: remove(Vertex_handle v ) { Vertex_handle u = v->up(); int l = 0 ; while(1) { hierarchy[l++]->remove(v); if (u == Vertex_handle()) break; if (l >= m_maxlevel) break; v = u; u = v->up(); } } template inline void Periodic_2_triangulation_hierarchy_2:: remove_degree_3(Vertex_handle v ) { remove(v); } template inline void Periodic_2_triangulation_hierarchy_2:: remove_first(Vertex_handle v ) { remove(v); } template typename Periodic_2_triangulation_hierarchy_2::Vertex_handle Periodic_2_triangulation_hierarchy_2:: move_if_no_collision(Vertex_handle v, const Point &p) { CGAL_triangulation_precondition(v != Vertex_handle()); Vertex_handle old, ret; for (int l = 0; l < m_maxlevel; ++l) { Vertex_handle u = v->up(); CGAL_triangulation_assertion(hierarchy[l]->is_valid()); Vertex_handle w = hierarchy[l]->move_if_no_collision(v, p); if (l == 0) { ret = w; } else { old->set_up(w); w->set_down(old); if (hierarchy[l]->number_of_sheets()[0] != 1) { std::vector vtc = hierarchy[l]->periodic_copies(w); for (unsigned int i = 0 ; i < vtc.size() ; i++) vtc[i]->set_down(old); } } if (u == Vertex_handle()) break; old = w; v = u; } return ret; } template typename Periodic_2_triangulation_hierarchy_2::Vertex_handle Periodic_2_triangulation_hierarchy_2:: move_point(Vertex_handle v, const Point &p) { CGAL_triangulation_precondition(v != Vertex_handle()); Vertex_handle old, ret; for (int l = 0; l < m_maxlevel; ++l) { Vertex_handle u = v->up(); CGAL_triangulation_assertion(hierarchy[l]->is_valid()); Vertex_handle w = hierarchy[l]->move_point(v, p); if (l == 0) { ret = w; } else { old->set_up(w); w->set_down(old); if (hierarchy[l]->number_of_sheets()[0] != 1) { std::vector vtc = hierarchy[l]->periodic_copies(w); for (unsigned int i = 0 ; i < vtc.size() ; i++) vtc[i]->set_down(old); } } if (u == Vertex_handle()) break; old = w; v = u; } return ret; } template typename Periodic_2_triangulation_hierarchy_2::Face_handle Periodic_2_triangulation_hierarchy_2:: locate(const Point& p, Locate_type& lt, int& li, Face_handle loc) const { Face_handle positions[m_maxlevel]; locate_in_all(p, lt, li, loc, positions); return positions[0]; } template typename Periodic_2_triangulation_hierarchy_2::Face_handle Periodic_2_triangulation_hierarchy_2:: locate(const Point& p, Face_handle loc ) const { if (loc != Face_handle ()) return PTr_Base::locate (p, loc); Locate_type lt; int li; return locate(p, lt, li, loc); } template void Periodic_2_triangulation_hierarchy_2:: locate_in_all(const Point& p, Locate_type& lt, int& li, Face_handle loc, Face_handle pos[m_maxlevel]) const { Face_handle position; Vertex_handle nearest; int level = m_maxlevel; typename Geom_traits::Compare_distance_2 closer = this->geom_traits().compare_distance_2_object(); // find the highest level with enough vertices that is at the same time 2D while ( (hierarchy[--level]->number_of_vertices() < static_cast (m_minsize )) || (hierarchy[level]->dimension() < 2) ) { if ( ! level) break; // do not go below 0 } if((level > 0) && (hierarchy[level]->dimension() < 2)) { level--; } for (int i = level + 1; i < m_maxlevel; ++i) pos[i] = 0; while(level > 0) { pos[level] = position = hierarchy[level]->locate(p, position); // locate at that level from "position" // result is stored in "position" for the next level // find the nearest between vertices 0 and 1 if (hierarchy[level]->is_infinite(position->vertex(0))) { nearest = position->vertex(1); } else if (hierarchy[level]->is_infinite(position->vertex(1))) { nearest = position->vertex(0); } else if ( closer(p, position->vertex(0)->point(), position->vertex(1)->point()) == SMALLER) { nearest = position->vertex(0); } else { nearest = position->vertex(1); } // compare to vertex 2, but only if the triangulation is 2D, because otherwise vertex(2) is nullptr if ( (hierarchy[level]->dimension() == 2) && (! hierarchy[level]->is_infinite(position->vertex(2)))) { if ( closer( p, position->vertex(2)->point(), nearest->point()) == SMALLER ) { nearest = position->vertex(2); } } // go at the same vertex on level below nearest = nearest->down(); position = nearest->face(); // incident face --level; } pos[0] = hierarchy[0]->locate(p, lt, li, loc == Face_handle() ? position : loc); // at level 0 } template int Periodic_2_triangulation_hierarchy_2:: random_level() { if ( level_mult_cover < m_maxlevel && hierarchy[level_mult_cover]->number_of_sheets() == make_array(1, 1) ) ++level_mult_cover; boost::geometric_distribution<> proba(1.0 / m_ratio); boost::variate_generator > die(random, proba); return (std::min)(die() - 1, level_mult_cover); } } //namespace CGAL #endif // CGAL_PERIODIC_2_TRIANGULATION_HIERARCHY_2_H