diff options
Diffstat (limited to 'include/gudhi_patches/CGAL/Regular_triangulation.h')
| -rw-r--r-- | include/gudhi_patches/CGAL/Regular_triangulation.h | 1169 |
1 files changed, 1169 insertions, 0 deletions
diff --git a/include/gudhi_patches/CGAL/Regular_triangulation.h b/include/gudhi_patches/CGAL/Regular_triangulation.h new file mode 100644 index 00000000..111c6ac9 --- /dev/null +++ b/include/gudhi_patches/CGAL/Regular_triangulation.h @@ -0,0 +1,1169 @@ +// Copyright (c) 2014 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) : Clement Jamin + +#ifndef CGAL_REGULAR_TRIANGULATION_H +#define CGAL_REGULAR_TRIANGULATION_H + +#include <CGAL/Triangulation.h> +#include <CGAL/Dimension.h> +#include <CGAL/Default.h> +#include <CGAL/spatial_sort.h> +#include <CGAL/Regular_triangulation_traits_adapter.h> + +#include <boost/property_map/function_property_map.hpp> + +namespace CGAL { + +template< typename Traits_, typename TDS_ = Default > +class Regular_triangulation +: public Triangulation< + Regular_triangulation_traits_adapter<Traits_>, + typename Default::Get< + TDS_, + Triangulation_data_structure< + typename Regular_triangulation_traits_adapter<Traits_>::Dimension, + Triangulation_vertex<Regular_triangulation_traits_adapter<Traits_> >, + Triangulation_full_cell<Regular_triangulation_traits_adapter<Traits_> > + > + >::type> +{ + typedef Regular_triangulation_traits_adapter<Traits_> RTTraits; + typedef typename RTTraits::Dimension Maximal_dimension_; + typedef typename Default::Get< + TDS_, + Triangulation_data_structure< + Maximal_dimension_, + Triangulation_vertex<RTTraits>, + Triangulation_full_cell<RTTraits> + > >::type TDS; + typedef Triangulation<RTTraits, TDS> Base; + typedef Regular_triangulation<Traits_, TDS_> Self; + + typedef typename RTTraits::Orientation_d Orientation_d; + typedef typename RTTraits::Power_side_of_power_sphere_d Power_side_of_power_sphere_d; + typedef typename RTTraits::In_flat_power_side_of_power_sphere_d + In_flat_power_side_of_power_sphere_d; + typedef typename RTTraits::Flat_orientation_d Flat_orientation_d; + typedef typename RTTraits::Construct_flat_orientation_d Construct_flat_orientation_d; + +public: // PUBLIC NESTED TYPES + + typedef RTTraits Geom_traits; + typedef typename Base::Triangulation_ds Triangulation_ds; + + typedef typename Base::Vertex Vertex; + typedef typename Base::Full_cell Full_cell; + typedef typename Base::Facet Facet; + typedef typename Base::Face Face; + + typedef Maximal_dimension_ Maximal_dimension; + typedef typename RTTraits::Bare_point_d Bare_point; + typedef typename RTTraits::Weighted_point_d Weighted_point; + + typedef typename Base::Point_const_iterator Point_const_iterator; + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Vertex_iterator Vertex_iterator; + typedef typename Base::Vertex_const_handle Vertex_const_handle; + typedef typename Base::Vertex_const_iterator Vertex_const_iterator; + + typedef typename Base::Full_cell_handle Full_cell_handle; + typedef typename Base::Full_cell_iterator Full_cell_iterator; + typedef typename Base::Full_cell_const_handle Full_cell_const_handle; + typedef typename Base::Full_cell_const_iterator Full_cell_const_iterator; + typedef typename Base::Finite_full_cell_const_iterator + Finite_full_cell_const_iterator; + + typedef typename Base::size_type size_type; + typedef typename Base::difference_type difference_type; + + typedef typename Base::Locate_type Locate_type; + + //Tag to distinguish Delaunay from Regular triangulations + typedef Tag_true Weighted_tag; + +protected: // DATA MEMBERS + + +public: + + using typename Base::Rotor; + using Base::maximal_dimension; + using Base::are_incident_full_cells_valid; + using Base::coaffine_orientation_predicate; + using Base::reset_flat_orientation; + using Base::current_dimension; + using Base::geom_traits; + using Base::index_of_covertex; + //using Base::index_of_second_covertex; + using Base::rotate_rotor; + using Base::infinite_vertex; + using Base::insert_in_hole; + using Base::is_infinite; + using Base::locate; + using Base::points_begin; + using Base::set_neighbors; + using Base::new_full_cell; + using Base::number_of_vertices; + using Base::orientation; + using Base::tds; + using Base::reorient_full_cells; + using Base::full_cell; + using Base::full_cells_begin; + using Base::full_cells_end; + using Base::finite_full_cells_begin; + using Base::finite_full_cells_end; + using Base::vertices_begin; + using Base::vertices_end; + +private: + + // Wrapper + struct Power_side_of_power_sphere_for_non_maximal_dim_d + { + boost::optional<Flat_orientation_d>* fop; + Construct_flat_orientation_d cfo; + In_flat_power_side_of_power_sphere_d ifpt; + + Power_side_of_power_sphere_for_non_maximal_dim_d( + boost::optional<Flat_orientation_d>& x, + Construct_flat_orientation_d const&y, + In_flat_power_side_of_power_sphere_d const&z) + : fop(&x), cfo(y), ifpt(z) {} + + template<class Iter> + CGAL::Orientation operator()(Iter a, Iter b, const Weighted_point & p)const + { + if(!*fop) + *fop=cfo(a,b); + return ifpt(fop->get(),a,b,p); + } + }; + +public: + +// - - - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS + + Regular_triangulation(int dim, const Geom_traits &k = Geom_traits()) + : Base(dim, k) + { + } + + // With this constructor, + // the user can specify a Flat_orientation_d object to be used for + // orienting simplices of a specific dimension + // (= preset_flat_orientation_.first) + // It it used by the dark triangulations created by DT::remove + Regular_triangulation( + int dim, + const std::pair<int, const Flat_orientation_d *> &preset_flat_orientation, + const Geom_traits &k = Geom_traits()) + : Base(dim, preset_flat_orientation, k) + { + } + + ~Regular_triangulation() {} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ACCESS + + // Not Documented + Power_side_of_power_sphere_for_non_maximal_dim_d power_side_of_power_sphere_for_non_maximal_dim_predicate() const + { + return Power_side_of_power_sphere_for_non_maximal_dim_d ( + flat_orientation_, + geom_traits().construct_flat_orientation_d_object(), + geom_traits().in_flat_power_side_of_power_sphere_d_object() + ); + } + + + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS + + // Warning: these functions are not correct since they do not restore hidden + // vertices + + Full_cell_handle remove(Vertex_handle); + Full_cell_handle remove(const Weighted_point & p, Full_cell_handle hint = Full_cell_handle()) + { + Locate_type lt; + Face f(maximal_dimension()); + Facet ft; + Full_cell_handle s = locate(p, lt, f, ft, hint); + if( Base::ON_VERTEX == lt ) + { + return remove(s->vertex(f.index(0))); + } + return Full_cell_handle(); + } + + template< typename ForwardIterator > + void remove(ForwardIterator start, ForwardIterator end) + { + while( start != end ) + remove(*start++); + } + + // Not documented + void remove_decrease_dimension(Vertex_handle); + + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INSERTIONS + + template< typename ForwardIterator > + std::ptrdiff_t insert(ForwardIterator start, ForwardIterator end) + { + size_type n = number_of_vertices(); + typedef std::vector<Weighted_point> WP_vec; + WP_vec points(start, end); + + spatial_sort(points.begin(), points.end(), geom_traits()); + + Full_cell_handle hint; + for(typename WP_vec::const_iterator p = points.begin(); p != points.end(); ++p ) + { + Locate_type lt; + Face f(maximal_dimension()); + Facet ft; + Full_cell_handle c = locate (*p, lt, f, ft, hint); + Vertex_handle v = insert (*p, lt, f, ft, c); + + hint = v == Vertex_handle() ? c : v->full_cell(); + } + return number_of_vertices() - n; + } + + Vertex_handle insert(const Weighted_point &, + Locate_type, + const Face &, + const Facet &, + Full_cell_handle); + + Vertex_handle insert(const Weighted_point & p, + Full_cell_handle start = Full_cell_handle()) + { + Locate_type lt; + Face f(maximal_dimension()); + Facet ft; + Full_cell_handle s = locate(p, lt, f, ft, start); + return insert(p, lt, f, ft, s); + } + + Vertex_handle insert(const Weighted_point & p, Vertex_handle hint) + { + CGAL_assertion( Vertex_handle() != hint ); + return insert(p, hint->full_cell()); + } + + Vertex_handle insert_outside_affine_hull(const Weighted_point &); + Vertex_handle insert_in_conflicting_cell( + const Weighted_point &, Full_cell_handle, + Vertex_handle only_if_this_vertex_is_in_the_cz = Vertex_handle()); + + Vertex_handle insert_if_in_star(const Weighted_point &, + Vertex_handle, + Locate_type, + const Face &, + const Facet &, + Full_cell_handle); + + Vertex_handle insert_if_in_star( + const Weighted_point & p, Vertex_handle star_center, + Full_cell_handle start = Full_cell_handle()) + { + Locate_type lt; + Face f(maximal_dimension()); + Facet ft; + Full_cell_handle s = locate(p, lt, f, ft, start); + return insert_if_in_star(p, star_center, lt, f, ft, s); + } + + Vertex_handle insert_if_in_star( + const Weighted_point & p, Vertex_handle star_center, + Vertex_handle hint) + { + CGAL_assertion( Vertex_handle() != hint ); + return insert_if_in_star(p, star_center, hint->full_cell()); + } + +// - - - - - - - - - - - - - - - - - - - - - - - - - GATHERING CONFLICTING SIMPLICES + + bool is_in_conflict(const Weighted_point &, Full_cell_const_handle) const; + + template< class OrientationPredicate > + Oriented_side perturbed_power_side_of_power_sphere(const Weighted_point &, + Full_cell_const_handle, const OrientationPredicate &) const; + + template< typename OutputIterator > + Facet compute_conflict_zone(const Weighted_point &, Full_cell_handle, OutputIterator) const; + + template < typename OrientationPredicate, typename PowerTestPredicate > + class Conflict_predicate + { + const Self & rt_; + const Weighted_point & p_; + OrientationPredicate ori_; + PowerTestPredicate power_side_of_power_sphere_; + int cur_dim_; + public: + Conflict_predicate( + const Self & rt, + const Weighted_point & p, + const OrientationPredicate & ori, + const PowerTestPredicate & power_side_of_power_sphere) + : rt_(rt), p_(p), ori_(ori), power_side_of_power_sphere_(power_side_of_power_sphere), cur_dim_(rt.current_dimension()) {} + + inline + bool operator()(Full_cell_const_handle s) const + { + bool ok; + if( ! rt_.is_infinite(s) ) + { + Oriented_side power_side_of_power_sphere = power_side_of_power_sphere_(rt_.points_begin(s), rt_.points_begin(s) + cur_dim_ + 1, p_); + if( ON_POSITIVE_SIDE == power_side_of_power_sphere ) + ok = true; + else if( ON_NEGATIVE_SIDE == power_side_of_power_sphere ) + ok = false; + else + ok = ON_POSITIVE_SIDE == rt_.perturbed_power_side_of_power_sphere<OrientationPredicate>(p_, s, ori_); + } + else + { + typedef typename Full_cell::Vertex_handle_const_iterator VHCI; + typedef Substitute_point_in_vertex_iterator<VHCI> F; + F spivi(rt_.infinite_vertex(), &p_); + + Orientation o = ori_( + boost::make_transform_iterator(s->vertices_begin(), spivi), + boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1, + spivi)); + + if( POSITIVE == o ) + ok = true; + else if( o == NEGATIVE ) + ok = false; + else + ok = (*this)(s->neighbor( s->index( rt_.infinite_vertex() ) )); + } + return ok; + } + }; + + template < typename ConflictPredicate > + class Conflict_traversal_predicate + { + const Self & rt_; + const ConflictPredicate & pred_; + public: + Conflict_traversal_predicate(const Self & rt, const ConflictPredicate & pred) + : rt_(rt), pred_(pred) + {} + inline + bool operator()(const Facet & f) const + { + return pred_(rt_.full_cell(f)->neighbor(rt_.index_of_covertex(f))); + } + }; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY + + bool is_valid(bool verbose = false, int level = 0) const; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MISC + + std::size_t number_of_hidden_vertices() const + { + return m_hidden_points.size(); + } + +private: + + template<typename InputIterator> + bool + does_cell_range_contain_vertex(InputIterator cz_begin, InputIterator cz_end, + Vertex_handle vh) const + { + // Check all vertices + while(cz_begin != cz_end) + { + Full_cell_handle fch = *cz_begin; + for (int i = 0 ; i <= current_dimension() ; ++i) + { + if (fch->vertex(i) == vh) + return true; + } + ++cz_begin; + } + return false; + } + + template<typename InputIterator, typename OutputIterator> + void + process_conflict_zone(InputIterator cz_begin, InputIterator cz_end, + OutputIterator vertices_out) const + { + // Get all vertices + while(cz_begin != cz_end) + { + Full_cell_handle fch = *cz_begin; + for (int i = 0 ; i <= current_dimension() ; ++i) + { + Vertex_handle vh = fch->vertex(i); + if (vh->full_cell() != Full_cell_handle()) + { + (*vertices_out++) = vh; + vh->set_full_cell(Full_cell_handle()); + } + } + ++cz_begin; + } + } + + + template<typename InputIterator> + void + process_cz_vertices_after_insertion(InputIterator vertices_begin, + InputIterator vertices_end) + { + // Get all vertices + while(vertices_begin != vertices_end) + { + Vertex_handle vh = *vertices_begin; + if (vh->full_cell() == Full_cell_handle()) + { + m_hidden_points.push_back(vh->point()); + tds().delete_vertex(vh); + } + ++vertices_begin; + } + } + +private: + // Some internal types to shorten notation + using typename Base::Coaffine_orientation_d; + using Base::flat_orientation_; + typedef Conflict_predicate<Coaffine_orientation_d, Power_side_of_power_sphere_for_non_maximal_dim_d> + Conflict_pred_in_subspace; + typedef Conflict_predicate<Orientation_d, Power_side_of_power_sphere_d> + Conflict_pred_in_fullspace; + typedef Conflict_traversal_predicate<Conflict_pred_in_subspace> + Conflict_traversal_pred_in_subspace; + typedef Conflict_traversal_predicate<Conflict_pred_in_fullspace> + Conflict_traversal_pred_in_fullspace; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MEMBER VARIABLES + std::vector<Weighted_point> m_hidden_points; + +}; // class Regular_triangulation + + +// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = +// FUNCTIONS THAT ARE MEMBER METHODS: + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS + + +// Warning: this function is not correct since it does not restore hidden +// vertices +template< typename Traits, typename TDS > +typename Regular_triangulation<Traits, TDS>::Full_cell_handle +Regular_triangulation<Traits, TDS> +::remove( Vertex_handle v ) +{ + CGAL_precondition( ! is_infinite(v) ); + CGAL_expensive_precondition( is_vertex(v) ); + + // THE CASE cur_dim == 0 + if( 0 == current_dimension() ) + { + remove_decrease_dimension(v); + return Full_cell_handle(); + } + else if( 1 == current_dimension() ) + { // THE CASE cur_dim == 1 + if( 2 == number_of_vertices() ) + { + remove_decrease_dimension(v); + return Full_cell_handle(); + } + Full_cell_handle left = v->full_cell(); + if( 0 == left->index(v) ) + left = left->neighbor(1); + CGAL_assertion( 1 == left->index(v) ); + Full_cell_handle right = left->neighbor(0); + tds().associate_vertex_with_full_cell(left, 1, right->vertex(1)); + set_neighbors(left, 0, right->neighbor(0), right->mirror_index(0)); + tds().delete_vertex(v); + tds().delete_full_cell(right); + return left; + } + + // THE CASE cur_dim >= 2 + // Gather the finite vertices sharing an edge with |v| + typedef typename Base::template Full_cell_set<Full_cell_handle> Simplices; + Simplices simps; + std::back_insert_iterator<Simplices> out(simps); + tds().incident_full_cells(v, out); + typedef std::set<Vertex_handle> Vertex_set; + Vertex_set verts; + Vertex_handle vh; + for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it ) + for( int i = 0; i <= current_dimension(); ++i ) + { + vh = (*it)->vertex(i); + if( is_infinite(vh) ) + continue; + if( vh == v ) + continue; + verts.insert(vh); + } + + // After gathering finite neighboring vertices, create their Dark Delaunay triangulation + typedef Triangulation_vertex<Geom_traits, Vertex_handle> Dark_vertex_base; + typedef Triangulation_full_cell< + Geom_traits, + internal::Triangulation::Dark_full_cell_data<TDS> > Dark_full_cell_base; + typedef Triangulation_data_structure<Maximal_dimension, + Dark_vertex_base, + Dark_full_cell_base + > Dark_tds; + typedef Regular_triangulation<Traits, Dark_tds> Dark_triangulation; + typedef typename Dark_triangulation::Face Dark_face; + typedef typename Dark_triangulation::Facet Dark_facet; + typedef typename Dark_triangulation::Vertex_handle Dark_v_handle; + typedef typename Dark_triangulation::Full_cell_handle Dark_s_handle; + + // If flat_orientation_ is defined, we give it the Dark triangulation + // so that the orientation it uses for "current_dimension()"-simplices is + // coherent with the global triangulation + Dark_triangulation dark_side( + maximal_dimension(), + flat_orientation_ ? + std::pair<int, const Flat_orientation_d *>(current_dimension(), flat_orientation_.get_ptr()) + : std::pair<int, const Flat_orientation_d *>(std::numeric_limits<int>::max(), NULL) ); + + Dark_s_handle dark_s; + Dark_v_handle dark_v; + typedef std::map<Vertex_handle, Dark_v_handle> Vertex_map; + Vertex_map light_to_dark; + typename Vertex_set::iterator vit = verts.begin(); + while( vit != verts.end() ) + { + dark_v = dark_side.insert((*vit)->point(), dark_s); + dark_s = dark_v->full_cell(); + dark_v->data() = *vit; + light_to_dark[*vit] = dark_v; + ++vit; + } + + if( dark_side.current_dimension() != current_dimension() ) + { + CGAL_assertion( dark_side.current_dimension() + 1 == current_dimension() ); + // Here, the finite neighbors of |v| span a affine subspace of + // dimension one less than the current dimension. Two cases are possible: + if( (size_type)(verts.size() + 1) == number_of_vertices() ) + { + remove_decrease_dimension(v); + return Full_cell_handle(); + } + else + { // |v| is strictly outside the convex hull of the rest of the points. This is an + // easy case: first, modify the finite full_cells, then, delete the infinite ones. + // We don't even need the Dark triangulation. + Simplices infinite_simps; + { + Simplices finite_simps; + for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it ) + if( is_infinite(*it) ) + infinite_simps.push_back(*it); + else + finite_simps.push_back(*it); + simps.swap(finite_simps); + } // now, simps only contains finite simplices + // First, modify the finite full_cells: + for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it ) + { + int v_idx = (*it)->index(v); + tds().associate_vertex_with_full_cell(*it, v_idx, infinite_vertex()); + } + // Make the handles to infinite full cells searchable + infinite_simps.make_searchable(); + // Then, modify the neighboring relation + for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it ) + { + for( int i = 0 ; i <= current_dimension(); ++i ) + { + if (is_infinite((*it)->vertex(i))) + continue; + (*it)->vertex(i)->set_full_cell(*it); + Full_cell_handle n = (*it)->neighbor(i); + // Was |n| a finite full cell prior to removing |v| ? + if( ! infinite_simps.contains(n) ) + continue; + int n_idx = n->index(v); + set_neighbors(*it, i, n->neighbor(n_idx), n->neighbor(n_idx)->index(n)); + } + } + Full_cell_handle ret_s; + // Then, we delete the infinite full_cells + for( typename Simplices::iterator it = infinite_simps.begin(); it != infinite_simps.end(); ++it ) + tds().delete_full_cell(*it); + tds().delete_vertex(v); + return simps.front(); + } + } + else // From here on, dark_side.current_dimension() == current_dimension() + { + dark_side.infinite_vertex()->data() = infinite_vertex(); + light_to_dark[infinite_vertex()] = dark_side.infinite_vertex(); + } + + // Now, compute the conflict zone of v->point() in + // the dark side. This is precisely the set of full_cells + // that we have to glue back into the light side. + Dark_face dark_f(dark_side.maximal_dimension()); + Dark_facet dark_ft; + typename Dark_triangulation::Locate_type lt; + dark_s = dark_side.locate(v->point(), lt, dark_f, dark_ft); + CGAL_assertion( lt != Dark_triangulation::ON_VERTEX + && lt != Dark_triangulation::OUTSIDE_AFFINE_HULL ); + + // |ret_s| is the full_cell that we return + Dark_s_handle dark_ret_s = dark_s; + Full_cell_handle ret_s; + + typedef typename Base::template Full_cell_set<Dark_s_handle> Dark_full_cells; + Dark_full_cells conflict_zone; + std::back_insert_iterator<Dark_full_cells> dark_out(conflict_zone); + + dark_ft = dark_side.compute_conflict_zone(v->point(), dark_s, dark_out); + // Make the dark simplices in the conflict zone searchable + conflict_zone.make_searchable(); + + // THE FOLLOWING SHOULD MAYBE GO IN TDS. + // Here is the plan: + // 1. Pick any Facet from boundary of the light zone + // 2. Find corresponding Facet on boundary of dark zone + // 3. stitch. + + // 1. Build a facet on the boudary of the light zone: + Full_cell_handle light_s = *simps.begin(); + Facet light_ft(light_s, light_s->index(v)); + + // 2. Find corresponding Dark_facet on boundary of the dark zone + Dark_full_cells dark_incident_s; + for( int i = 0; i <= current_dimension(); ++i ) + { + if( index_of_covertex(light_ft) == i ) + continue; + Dark_v_handle dark_v = light_to_dark[full_cell(light_ft)->vertex(i)]; + dark_incident_s.clear(); + dark_out = std::back_inserter(dark_incident_s); + dark_side.tds().incident_full_cells(dark_v, dark_out); + for(typename Dark_full_cells::iterator it = dark_incident_s.begin(); + it != dark_incident_s.end(); + ++it) + { + (*it)->data().count_ += 1; + } + } + + for( typename Dark_full_cells::iterator it = dark_incident_s.begin(); it != dark_incident_s.end(); ++it ) + { + if( current_dimension() != (*it)->data().count_ ) + continue; + if( ! conflict_zone.contains(*it) ) + continue; + // We found a full_cell incident to the dark facet corresponding to the light facet |light_ft| + int ft_idx = 0; + while( light_s->has_vertex( (*it)->vertex(ft_idx)->data() ) ) + ++ft_idx; + dark_ft = Dark_facet(*it, ft_idx); + break; + } + // Pre-3. Now, we are ready to traverse both boundary and do the stiching. + + // But first, we create the new full_cells in the light triangulation, + // with as much adjacency information as possible. + + // Create new full_cells with vertices + for( typename Dark_full_cells::iterator it = conflict_zone.begin(); it != conflict_zone.end(); ++it ) + { + Full_cell_handle new_s = new_full_cell(); + (*it)->data().light_copy_ = new_s; + for( int i = 0; i <= current_dimension(); ++i ) + tds().associate_vertex_with_full_cell(new_s, i, (*it)->vertex(i)->data()); + if( dark_ret_s == *it ) + ret_s = new_s; + } + + // Setup adjacencies inside the hole + for( typename Dark_full_cells::iterator it = conflict_zone.begin(); it != conflict_zone.end(); ++it ) + { + Full_cell_handle new_s = (*it)->data().light_copy_; + for( int i = 0; i <= current_dimension(); ++i ) + if( conflict_zone.contains((*it)->neighbor(i)) ) + tds().set_neighbors(new_s, i, (*it)->neighbor(i)->data().light_copy_, (*it)->mirror_index(i)); + } + + // 3. Stitch + simps.make_searchable(); + typedef std::queue<std::pair<Facet, Dark_facet> > Queue; + Queue q; + q.push(std::make_pair(light_ft, dark_ft)); + dark_s = dark_side.full_cell(dark_ft); + int dark_i = dark_side.index_of_covertex(dark_ft); + // mark dark_ft as visited: + // TODO try by marking with Dark_v_handle (vertex) + dark_s->neighbor(dark_i)->set_neighbor(dark_s->mirror_index(dark_i), Dark_s_handle()); + while( ! q.empty() ) + { + std::pair<Facet, Dark_facet> p = q.front(); + q.pop(); + light_ft = p.first; + dark_ft = p.second; + light_s = full_cell(light_ft); + int light_i = index_of_covertex(light_ft); + dark_s = dark_side.full_cell(dark_ft); + int dark_i = dark_side.index_of_covertex(dark_ft); + Full_cell_handle light_n = light_s->neighbor(light_i); + set_neighbors(dark_s->data().light_copy_, dark_i, light_n, light_s->mirror_index(light_i)); + for( int di = 0; di <= current_dimension(); ++di ) + { + if( di == dark_i ) + continue; + int li = light_s->index(dark_s->vertex(di)->data()); + Rotor light_r(light_s, li, light_i); + typename Dark_triangulation::Rotor dark_r(dark_s, di, dark_i); + + while( simps.contains(cpp11::get<0>(light_r)->neighbor(cpp11::get<1>(light_r))) ) + light_r = rotate_rotor(light_r); + + while( conflict_zone.contains(cpp11::get<0>(dark_r)->neighbor(cpp11::get<1>(dark_r))) ) + dark_r = dark_side.rotate_rotor(dark_r); + + Dark_s_handle dark_ns = cpp11::get<0>(dark_r); + int dark_ni = cpp11::get<1>(dark_r); + Full_cell_handle light_ns = cpp11::get<0>(light_r); + int light_ni = cpp11::get<1>(light_r); + // mark dark_r as visited: + // TODO try by marking with Dark_v_handle (vertex) + Dark_s_handle outside = dark_ns->neighbor(dark_ni); + Dark_v_handle mirror = dark_ns->mirror_vertex(dark_ni, current_dimension()); + int dn = outside->index(mirror); + if( Dark_s_handle() == outside->neighbor(dn) ) + continue; + outside->set_neighbor(dn, Dark_s_handle()); + q.push(std::make_pair(Facet(light_ns, light_ni), Dark_facet(dark_ns, dark_ni))); + } + } + tds().delete_full_cells(simps.begin(), simps.end()); + tds().delete_vertex(v); + return ret_s; +} + +template< typename Traits, typename TDS > +void +Regular_triangulation<Traits, TDS> +::remove_decrease_dimension(Vertex_handle v) +{ + CGAL_precondition( current_dimension() >= 0 ); + tds().remove_decrease_dimension(v, infinite_vertex()); + // reset the predicates: + reset_flat_orientation(); + if( 1 <= current_dimension() ) + { + Full_cell_handle inf_v_cell = infinite_vertex()->full_cell(); + int inf_v_index = inf_v_cell->index(infinite_vertex()); + Full_cell_handle s = inf_v_cell->neighbor(inf_v_index); + Orientation o = orientation(s); + CGAL_assertion( ZERO != o ); + if( NEGATIVE == o ) + reorient_full_cells(); + } +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INSERTIONS + +template< typename Traits, typename TDS > +typename Regular_triangulation<Traits, TDS>::Vertex_handle +Regular_triangulation<Traits, TDS> +::insert(const Weighted_point & p, Locate_type lt, const Face & f, const Facet & ft, Full_cell_handle s) +{ + switch( lt ) + { + case Base::OUTSIDE_AFFINE_HULL: + return insert_outside_affine_hull(p); + break; + case Base::ON_VERTEX: + { + Vertex_handle v = s->vertex(f.index(0)); + typename RTTraits::Compute_weight_d pw = + geom_traits().compute_weight_d_object(); + + if (pw(p) == pw(v->point())) + return v; + // If dim == 0 and the new point has a bigger weight, + // we just replace the point, and the former point gets hidden + else if (current_dimension() == 0) + { + if (pw(p) > pw(v->point())) + { + m_hidden_points.push_back(v->point()); + v->set_point(p); + return v; + } + // Otherwise, the new point is hidden + else + { + m_hidden_points.push_back(p); + return Vertex_handle(); + } + } + // Otherwise, we apply the "normal" algorithm + + // !NO break here! + } + default: + return insert_in_conflicting_cell(p, s); + } +} + +/* +Inserts the point `p` in the regular triangulation. Returns a handle to the +newly created vertex at that position. +\pre The point `p` +must lie outside the affine hull of the regular triangulation. This implies that +`rt`.`current_dimension()` must be smaller than `rt`.`maximal_dimension()`. +*/ +template< typename Traits, typename TDS > +typename Regular_triangulation<Traits, TDS>::Vertex_handle +Regular_triangulation<Traits, TDS> +::insert_outside_affine_hull(const Weighted_point & p) +{ + // we don't use Base::insert_outside_affine_hull(...) because here, we + // also need to reset the side_of_oriented_subsphere functor. + CGAL_precondition( current_dimension() < maximal_dimension() ); + Vertex_handle v = tds().insert_increase_dimension(infinite_vertex()); + // reset the predicates: + reset_flat_orientation(); + v->set_point(p); + if( current_dimension() >= 1 ) + { + Full_cell_handle inf_v_cell = infinite_vertex()->full_cell(); + int inf_v_index = inf_v_cell->index(infinite_vertex()); + Full_cell_handle s = inf_v_cell->neighbor(inf_v_index); + Orientation o = orientation(s); + CGAL_assertion( ZERO != o ); + if( NEGATIVE == o ) + reorient_full_cells(); + + // We just inserted the second finite point and the right infinite + // cell is like : (inf_v, v), but we want it to be (v, inf_v) to be + // consistent with the rest of the cells + if (current_dimension() == 1) + { + // Is "inf_v_cell" the right infinite cell? Then inf_v_index should be 1 + if (inf_v_cell->neighbor(inf_v_index)->index(inf_v_cell) == 0 + && inf_v_index == 0) + { + inf_v_cell->swap_vertices(current_dimension() - 1, current_dimension()); + } + else + { + inf_v_cell = inf_v_cell->neighbor((inf_v_index + 1) % 2); + inf_v_index = inf_v_cell->index(infinite_vertex()); + // Is "inf_v_cell" the right infinite cell? Then inf_v_index should be 1 + if (inf_v_cell->neighbor(inf_v_index)->index(inf_v_cell) == 0 + && inf_v_index == 0) + { + inf_v_cell->swap_vertices(current_dimension() - 1, current_dimension()); + } + } + } + } + return v; +} + +template< typename Traits, typename TDS > +typename Regular_triangulation<Traits, TDS>::Vertex_handle +Regular_triangulation<Traits, TDS> +::insert_if_in_star(const Weighted_point & p, + Vertex_handle star_center, + Locate_type lt, + const Face & f, + const Facet & ft, + Full_cell_handle s) +{ + switch( lt ) + { + case Base::OUTSIDE_AFFINE_HULL: + return insert_outside_affine_hull(p); + break; + case Base::ON_VERTEX: + { + Vertex_handle v = s->vertex(f.index(0)); + typename RTTraits::Compute_weight_d pw = + geom_traits().compute_weight_d_object(); + if (pw(p) == pw(v->point())) + return v; + // If dim == 0 and the new point has a bigger weight, + // we replace the point + else if (current_dimension() == 0) + { + if (pw(p) > pw(v->point())) + v->set_point(p); + else + return v; + } + // Otherwise, we apply the "normal" algorithm + + // !NO break here! + } + default: + return insert_in_conflicting_cell(p, s, star_center); + } + + return Vertex_handle(); +} + +/* +[Undocumented function] + +Inserts the point `p` in the regular triangulation. `p` must be +in conflict with the second parameter `c`, which is used as a +starting point for `compute_conflict_zone`. +The function is faster than the standard `insert` function since +it does not need to call `locate`. + +If this insertion creates a vertex, this vertex is returned. + +If `p` coincides with an existing vertex and has a greater weight, +then the existing weighted point becomes hidden and `p` replaces it as vertex +of the triangulation. + +If `p` coincides with an already existing vertex (both point and +weights being equal), then this vertex is returned and the triangulation +remains unchanged. + +Otherwise if `p` does not appear as a vertex of the triangulation, +then it is stored as a hidden point and this method returns the default +constructed handle. + +\pre The point `p` must be in conflict with the full cell `c`. +*/ + +template< typename Traits, typename TDS > +typename Regular_triangulation<Traits, TDS>::Vertex_handle +Regular_triangulation<Traits, TDS> +::insert_in_conflicting_cell(const Weighted_point & p, + Full_cell_handle s, + Vertex_handle only_if_this_vertex_is_in_the_cz) +{ + typedef std::vector<Full_cell_handle> Full_cell_h_vector; + + bool in_conflict = is_in_conflict(p, s); + + // If p is not in conflict with s, then p is hidden + // => we don't insert it + if (!in_conflict) + { + m_hidden_points.push_back(p); + return Vertex_handle(); + } + else + { + Full_cell_h_vector cs; // for storing conflicting full_cells. + cs.reserve(64); + std::back_insert_iterator<Full_cell_h_vector> out(cs); + Facet ft = compute_conflict_zone(p, s, out); + + // Check if the CZ contains "only_if_this_vertex_is_in_the_cz" + if (only_if_this_vertex_is_in_the_cz != Vertex_handle() + && !does_cell_range_contain_vertex(cs.begin(), cs.end(), + only_if_this_vertex_is_in_the_cz)) + { + return Vertex_handle(); + } + + // Otherwise, proceed with the insertion + std::vector<Vertex_handle> cz_vertices; + cz_vertices.reserve(64); + process_conflict_zone(cs.begin(), cs.end(), + std::back_inserter(cz_vertices)); + + Vertex_handle ret = insert_in_hole(p, cs.begin(), cs.end(), ft); + + process_cz_vertices_after_insertion(cz_vertices.begin(), cz_vertices.end()); + + return ret; + } +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - GATHERING CONFLICTING SIMPLICES + +// NOT DOCUMENTED +template< typename Traits, typename TDS > +template< typename OrientationPred > +Oriented_side +Regular_triangulation<Traits, TDS> +::perturbed_power_side_of_power_sphere(const Weighted_point & p, Full_cell_const_handle s, + const OrientationPred & ori) const +{ + CGAL_precondition_msg( ! is_infinite(s), "full cell must be finite"); + CGAL_expensive_precondition( POSITIVE == orientation(s) ); + typedef std::vector<const Weighted_point *> Points; + Points points(current_dimension() + 2); + int i(0); + for( ; i <= current_dimension(); ++i ) + points[i] = &(s->vertex(i)->point()); + points[i] = &p; + std::sort(points.begin(), points.end(), + internal::Triangulation::Compare_points_for_perturbation<Self>(*this)); + typename Points::const_reverse_iterator cut_pt = points.rbegin(); + Points test_points; + while( cut_pt != points.rend() ) + { + if( &p == *cut_pt ) + // because the full_cell "s" is assumed to be positively oriented + return ON_NEGATIVE_SIDE; // we consider |p| to lie outside the sphere + test_points.clear(); + Point_const_iterator spit = points_begin(s); + int adjust_sign = -1; + for( i = 0; i < current_dimension(); ++i ) + { + if( &(*spit) == *cut_pt ) + { + ++spit; + adjust_sign = (((current_dimension() + i) % 2) == 0) ? -1 : +1; + } + test_points.push_back(&(*spit)); + ++spit; + } + test_points.push_back(&p); + + typedef typename CGAL::Iterator_project< + typename Points::iterator, + internal::Triangulation::Point_from_pointer<Self>, + const Weighted_point &, const Weighted_point * + > Point_pointer_iterator; + + Orientation ori_value = ori( + Point_pointer_iterator(test_points.begin()), + Point_pointer_iterator(test_points.end())); + + if( ZERO != ori_value ) + return Oriented_side( - adjust_sign * ori_value ); + + ++cut_pt; + } + CGAL_assertion(false); // we should never reach here + return ON_NEGATIVE_SIDE; +} + +template< typename Traits, typename TDS > +bool +Regular_triangulation<Traits, TDS> +::is_in_conflict(const Weighted_point & p, Full_cell_const_handle s) const +{ + CGAL_precondition( 1 <= current_dimension() ); + if( current_dimension() < maximal_dimension() ) + { + Conflict_pred_in_subspace c( + *this, p, + coaffine_orientation_predicate(), + power_side_of_power_sphere_for_non_maximal_dim_predicate()); + return c(s); + } + else + { + Orientation_d ori = geom_traits().orientation_d_object(); + Power_side_of_power_sphere_d side = geom_traits().power_side_of_power_sphere_d_object(); + Conflict_pred_in_fullspace c(*this, p, ori, side); + return c(s); + } +} + +template< typename Traits, typename TDS > +template< typename OutputIterator > +typename Regular_triangulation<Traits, TDS>::Facet +Regular_triangulation<Traits, TDS> +::compute_conflict_zone(const Weighted_point & p, Full_cell_handle s, OutputIterator out) const +{ + CGAL_precondition( 1 <= current_dimension() ); + if( current_dimension() < maximal_dimension() ) + { + Conflict_pred_in_subspace c( + *this, p, + coaffine_orientation_predicate(), + power_side_of_power_sphere_for_non_maximal_dim_predicate()); + Conflict_traversal_pred_in_subspace tp(*this, c); + return tds().gather_full_cells(s, tp, out); + } + else + { + Orientation_d ori = geom_traits().orientation_d_object(); + Power_side_of_power_sphere_d side = geom_traits().power_side_of_power_sphere_d_object(); + Conflict_pred_in_fullspace c(*this, p, ori, side); + Conflict_traversal_pred_in_fullspace tp(*this, c); + return tds().gather_full_cells(s, tp, out); + } +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY + +template< typename Traits, typename TDS > +bool +Regular_triangulation<Traits, TDS> +::is_valid(bool verbose, int level) const +{ + if (!Base::is_valid(verbose, level)) + return false; + + int dim = current_dimension(); + if (dim == maximal_dimension()) + { + for (Finite_full_cell_const_iterator cit = finite_full_cells_begin() ; + cit != finite_full_cells_end() ; ++cit ) + { + Full_cell_const_handle ch = cit.base(); + for(int i = 0; i < dim+1 ; ++i ) + { + // If the i-th neighbor is not an infinite cell + Vertex_handle opposite_vh = + ch->neighbor(i)->vertex(ch->neighbor(i)->index(ch)); + if (!is_infinite(opposite_vh)) + { + Power_side_of_power_sphere_d side = + geom_traits().power_side_of_power_sphere_d_object(); + if (side(Point_const_iterator(ch->vertices_begin()), + Point_const_iterator(ch->vertices_end()), + opposite_vh->point()) == ON_POSITIVE_SIDE) + { + if (verbose) + CGAL_warning_msg(false, "Non-empty sphere"); + return false; + } + } + } + } + } + return true; +} + +} //namespace CGAL + +#endif //CGAL_REGULAR_TRIANGULATION_H |