diff --git a/NewKernel_d/include/CGAL/Epick_d.h b/NewKernel_d/include/CGAL/Epick_d.h index cf780630c56..6443853994f 100644 --- a/NewKernel_d/include/CGAL/Epick_d.h +++ b/NewKernel_d/include/CGAL/Epick_d.h @@ -26,6 +26,7 @@ #include #include #include +#include namespace CGAL { diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_base.h b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_base.h index 28893d64830..80ee4b33580 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_base.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_base.h @@ -87,6 +87,7 @@ struct Cartesian_LA_base_d : public Dimension_base ::add::type ::add::type ::add::type + ::add::type Object_list; typedef typeset< Point_cartesian_const_iterator_tag>::type diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Coaffine.h b/NewKernel_d/include/CGAL/NewKernel_d/Coaffine.h index 1c2efe462d8..62ad90afa9c 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Coaffine.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Coaffine.h @@ -89,6 +89,7 @@ template struct Construct_flat_orientation : private Store_kernel std::vector& rest=o.rest; rest.reserve(dim+1); for(int i=0; i +#include +#include + +#include + +typedef CGAL::Epick_d K; +typedef CGAL::Delaunay_triangulation DT; + +void test(int dim) +{ + std::stringstream input_filename; + input_filename << "data/points_" << dim << ".cin"; + std::ifstream in(input_filename.str()); + + DT::Point p; + std::vector points; + + int dim_from_file; + in >> dim_from_file; + while(in >> p) + points.push_back(p); + + // Build the Regular Triangulation + DT dt(dim_from_file); + dt.insert(points.begin(), points.end()); + CGAL_assertion(dt.is_valid(true)); + + // Export + std::stringstream output_filename; + output_filename << "data/dt_dim" << dim << ".off"; + std::ofstream off_stream(output_filename.str()); + CGAL::export_triangulation_to_off(off_stream, dt); +} + +int main() +{ + //test(2); + //test(3); + test(10); + return 0; +} diff --git a/Triangulation/applications/Triangulation/points_to_RT_to_off.cpp b/Triangulation/applications/Triangulation/points_to_RT_to_off.cpp new file mode 100644 index 00000000000..94a633a505d --- /dev/null +++ b/Triangulation/applications/Triangulation/points_to_RT_to_off.cpp @@ -0,0 +1,43 @@ +#include +#include +#include +#include + +#include + +typedef CGAL::Epick_d K; +typedef CGAL::Regular_triangulation_euclidean_traits Traits; +typedef CGAL::Regular_triangulation RT; + +void test(int dim) +{ + std::stringstream input_filename; + input_filename << "data/points_" << dim << ".cin"; + std::ifstream in(input_filename.str()); + + RT::Weighted_point wp; + std::vector wpoints; + + int dim_from_file; + in >> dim_from_file; + while(in >> wp) + wpoints.push_back(wp); + + // Build the Regular Triangulation + RT rt(dim_from_file); + rt.insert(wpoints.begin(), wpoints.end()); + CGAL_assertion(rt.is_valid(true)); + + // Export + std::stringstream output_filename; + output_filename << "data/rt_dim" << dim << ".off"; + std::ofstream off_stream(output_filename.str()); + CGAL::export_triangulation_to_off(off_stream, rt); +} + +int main() +{ + test(2); + test(3); + return 0; +} diff --git a/Triangulation/benchmark/Triangulation/delaunay.cpp b/Triangulation/benchmark/Triangulation/delaunay.cpp index ea6c160c807..7d1a5b379ee 100644 --- a/Triangulation/benchmark/Triangulation/delaunay.cpp +++ b/Triangulation/benchmark/Triangulation/delaunay.cpp @@ -65,6 +65,5 @@ int main(int argc, char **argv) go<7>(N); go<8>(N); - return 0; } diff --git a/Triangulation/dont_submit b/Triangulation/dont_submit index f5c7a0a6d72..b2de9879971 100644 --- a/Triangulation/dont_submit +++ b/Triangulation/dont_submit @@ -1,3 +1,2 @@ TODO include/CGAL/Convex_hull.h -include/CGAL/Regular_triangulation.h diff --git a/Triangulation/include/CGAL/Delaunay_triangulation.h b/Triangulation/include/CGAL/Delaunay_triangulation.h index 21e0cc8d292..75eabc243e0 100644 --- a/Triangulation/include/CGAL/Delaunay_triangulation.h +++ b/Triangulation/include/CGAL/Delaunay_triangulation.h @@ -55,7 +55,7 @@ class Delaunay_triangulation public: // PUBLIC NESTED TYPES typedef DCTraits Geom_traits; - typedef typename Base::Triangulation_ds Triangulation_ds; + typedef typename Base::Triangulation_ds Triangulation_ds; typedef typename Base::Vertex Vertex; typedef typename Base::Full_cell Full_cell; @@ -71,21 +71,25 @@ public: // PUBLIC NESTED TYPES 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::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::size_type size_type; typedef typename Base::difference_type difference_type; typedef typename Base::Locate_type Locate_type; + //Tag to distinguish triangulations with weighted_points + typedef Tag_false 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; @@ -95,11 +99,11 @@ public: //using Base::incident_full_cells; using Base::geom_traits; using Base::index_of_covertex; + //using Base::index_of_second_covertex; using Base::infinite_vertex; using Base::insert_in_hole; using Base::insert_outside_convex_hull_1; using Base::is_infinite; - using Base::is_valid; using Base::locate; using Base::points_begin; using Base::set_neighbors; @@ -143,36 +147,9 @@ private: }; public: -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - UTILITIES - - // A co-dimension 2 sub-simplex. called a Rotor because we can rotate - // the two "covertices" around the sub-simplex. Useful for traversing the - // boundary of a hole. NOT DOCUMENTED - typedef cpp11::tuple Rotor; - - /*Full_cell_handle full_cell(const Rotor & r) const // NOT DOCUMENTED - { - return cpp11::get<0>(r); - } - int index_of_covertex(const Rotor & r) const // NOT DOCUMENTED - { - return cpp11::get<1>(r); - } - int index_of_second_covertex(const Rotor & r) const // NOT DOCUMENTED - { - return cpp11::get<2>(r); - }*/ - Rotor rotate_rotor(Rotor & r) // NOT DOCUMENTED... - { - int opposite = cpp11::get<0>(r)->mirror_index(cpp11::get<1>(r)); - Full_cell_handle s = cpp11::get<0>(r)->neighbor(cpp11::get<1>(r)); - int new_second = s->index(cpp11::get<0>(r)->vertex(cpp11::get<2>(r))); - return Rotor(s, new_second, opposite); - } - // - - - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS - Delaunay_triangulation(int dim, const Geom_traits k = Geom_traits()) + Delaunay_triangulation(int dim, const Geom_traits &k = Geom_traits()) : Base(dim, k) { } @@ -185,7 +162,7 @@ public: Delaunay_triangulation( int dim, const std::pair &preset_flat_orientation, - const Geom_traits k = Geom_traits()) + const Geom_traits &k = Geom_traits()) : Base(dim, preset_flat_orientation, k) { } @@ -340,6 +317,10 @@ public: return pred_(dc_.full_cell(f)->neighbor(dc_.index_of_covertex(f))); } }; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY + + bool is_valid(bool verbose = false, int level = 0) const; private: // Some internal types to shorten notation @@ -353,27 +334,6 @@ private: Conflict_traversal_pred_in_subspace; typedef Conflict_traversal_predicate Conflict_traversal_pred_in_fullspace; - - // This is used in the |remove(v)| member function to manage sets of Full_cell_handles - template< typename FCH > - struct Full_cell_set : public std::vector - { - typedef std::vector Base_set; - using Base_set::begin; - using Base_set::end; - void make_searchable() - { // sort the full cell handles - std::sort(begin(), end()); - } - bool contains(const FCH & fch) const - { - return std::binary_search(begin(), end(), fch); - } - bool contains_1st_and_not_2nd(const FCH & fst, const FCH & snd) const - { - return ( ! contains(snd) ) && ( contains(fst) ); - } - }; }; // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = @@ -428,7 +388,7 @@ Delaunay_triangulation // THE CASE cur_dim >= 2 // Gather the finite vertices sharing an edge with |v| - typedef Full_cell_set Simplices; + typedef typename Base::template Full_cell_set Simplices; Simplices simps; std::back_insert_iterator out(simps); tds().incident_full_cells(v, out); @@ -564,7 +524,7 @@ Delaunay_triangulation Dark_s_handle dark_ret_s = dark_s; Full_cell_handle ret_s; - typedef Full_cell_set Dark_full_cells; + typedef typename Base::template Full_cell_set Dark_full_cells; Dark_full_cells conflict_zone; std::back_insert_iterator dark_out(conflict_zone); @@ -780,9 +740,8 @@ Delaunay_triangulation ::insert_in_conflicting_cell(const Point & p, const Full_cell_handle s) { typedef std::vector Full_cell_h_vector; - static Full_cell_h_vector cs; // for storing conflicting full_cells. - cs.clear(); - // cs.reserve(64); + Full_cell_h_vector cs; // for storing conflicting full_cells. + cs.reserve(64); std::back_insert_iterator out(cs); Facet ft = compute_conflict_zone(p, s, out); return insert_in_hole(p, cs.begin(), cs.end(), ft); @@ -890,6 +849,48 @@ Delaunay_triangulation } } +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY + +template< typename DCTraits, typename TDS > +bool +Delaunay_triangulation +::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)) + { + Side_of_oriented_sphere_d side = + geom_traits().side_of_oriented_sphere_d_object(); + if (side(Point_const_iterator(ch->vertices_begin()), + Point_const_iterator(ch->vertices_end()), + opposite_vh->point()) == ON_BOUNDED_SIDE) + { + if (verbose) + CGAL_warning_msg(false, "Non-empty sphere"); + return false; + } + } + } + } + } + return true; +} + + } //namespace CGAL #endif // CGAL_DELAUNAY_COMPLEX_H diff --git a/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h b/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h new file mode 100644 index 00000000000..590c2a6b0e6 --- /dev/null +++ b/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h @@ -0,0 +1,255 @@ +// 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 Lesser 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_TRIANGULATION_IO_H +#define CGAL_TRIANGULATION_IO_H + +#include +#include +#include +#include + +namespace CGAL { + +namespace Triangulation_IO +{ +// TODO: test if the stream is binary or text? +template +void +output_point(std::ostream & os, const Traits &traits, const P & p) +{ + typedef typename Traits::Compute_coordinate_d Ccd; + const Ccd ccd = traits.compute_coordinate_d_object(); + const int dim = traits.point_dimension_d_object()(p); + if (dim > 0) + { + os << ccd(p, 0); + for (int i = 1 ; i < dim ; ++i) + os << " " << CGAL::to_double(ccd(p, i)); + } +} + +// TODO: test if the stream is binary or text? +/*template +void +input_point(std::istream & is, const Traits &traits, P & p) +{ + typedef typename Traits::FT FT; + std::vector coords; + + std::string line; + for(;;) + { + if (!std::getline(is, line)) + return is; + if (line != "") + break; + } + std::stringstream line_sstr(line); + FT temp; + while (line_sstr >> temp) + coords.push_back(temp); + + p = traits.construct_point_d_object()(coords.begin(), coords.end()); +}*/ + +} // namespace Triangulation_IO + +/////////////////////////////////////////////////////////////// +// TODO: replace these operator>> by an "input_point" function +/////////////////////////////////////////////////////////////// + +// TODO: test if the stream is binary or text? +template +std::istream & +operator>>(std::istream &is, typename Wrap::Point_d & p) +{ + typedef typename Wrap::Point_d P; + typedef typename K::FT FT; + std::vector coords; + + std::string line; + for(;;) + { + if (!std::getline(is, line)) + return is; + if (line != "") + break; + } + std::stringstream line_sstr(line); + FT temp; + while (line_sstr >> temp) + coords.push_back(temp); + + p = P(coords.begin(), coords.end()); + return is; +} + +// TODO: test if the stream is binary or text? +template +std::istream & +operator>>(std::istream &is, typename Wrap::Weighted_point_d & wp) +{ + typedef typename Wrap::Point_d P; + typedef typename Wrap::Weighted_point_d WP; + typedef typename K::FT FT; + + std::string line; + for(;;) + { + if (!std::getline(is, line)) + return is; + if (line != "") + break; + } + std::stringstream line_sstr(line); + FT temp; + std::vector coords; + while (line_sstr >> temp) + coords.push_back(temp); + + std::vector::iterator last = coords.end() - 1; + P p = P(coords.begin(), last); + wp = WP(p, *last); + + return is; +} + +template < class GT, class TDS > +std::ostream & +export_triangulation_to_off(std::ostream & os, + const Triangulation & tr, + bool in_3D_export_surface_only = false) +{ + typedef Triangulation Tr; + typedef typename Tr::Vertex_const_handle Vertex_handle; + typedef typename Tr::Vertex_const_iterator Vertex_iterator; + typedef typename Tr::Finite_vertex_const_iterator Finite_vertex_iterator; + typedef typename Tr::Full_cell_const_handle Full_cell_handle; + typedef typename Tr::Finite_full_cell_const_iterator Finite_full_cell_iterator; + typedef typename Tr::Full_cell_const_iterator Full_cell_iterator; + typedef typename Tr::Full_cell Full_cell; + typedef typename Full_cell::Vertex_handle_const_iterator Full_cell_vertex_iterator; + + if (tr.maximal_dimension() < 2 || tr.maximal_dimension() > 3) + { + std::cerr << "Warning: export_tds_to_off => dimension should be 2 or 3."; + os << "Warning: export_tds_to_off => dimension should be 2 or 3."; + return os; + } + + size_t n = tr.number_of_vertices(); + + std::stringstream output; + + // write the vertices + std::map index_of_vertex; + int i = 0; + for(Finite_vertex_iterator it = tr.finite_vertices_begin(); + it != tr.finite_vertices_end(); ++it, ++i) + { + Triangulation_IO::output_point(output, tr.geom_traits(), it->point()); + if (tr.maximal_dimension() == 2) + output << " 0"; + output << std::endl; + index_of_vertex[it.base()] = i; + } + CGAL_assertion( i == n ); + + size_t number_of_triangles = 0; + if (tr.maximal_dimension() == 2) + { + for (Finite_full_cell_iterator fch = tr.finite_full_cells_begin() ; + fch != tr.finite_full_cells_end() ; ++fch) + { + output << "3 "; + for (Full_cell_vertex_iterator vit = fch->vertices_begin() ; + vit != fch->vertices_end() ; ++vit) + { + output << index_of_vertex[*vit] << " "; + } + output << std::endl; + ++number_of_triangles; + } + } + else if (tr.maximal_dimension() == 3) + { + if (in_3D_export_surface_only) + { + // Parse boundary facets + for (Full_cell_iterator fch = tr.full_cells_begin() ; + fch != tr.full_cells_end() ; ++fch) + { + if (tr.is_infinite(fch)) + { + output << "3 "; + for (Full_cell_vertex_iterator vit = fch->vertices_begin() ; + vit != fch->vertices_end() ; ++vit) + { + if (!tr.is_infinite(*vit)) + output << index_of_vertex[*vit] << " "; + } + output << std::endl; + ++number_of_triangles; + } + } + } + else + { + // Parse finite cells + for (Finite_full_cell_iterator fch = tr.finite_full_cells_begin() ; + fch != tr.finite_full_cells_end() ; ++fch) + { + output << "3 " + << index_of_vertex[fch->vertex(0)] << " " + << index_of_vertex[fch->vertex(1)] << " " + << index_of_vertex[fch->vertex(2)] + << std::endl; + output << "3 " + << index_of_vertex[fch->vertex(0)] << " " + << index_of_vertex[fch->vertex(2)] << " " + << index_of_vertex[fch->vertex(3)] + << std::endl; + output << "3 " + << index_of_vertex[fch->vertex(1)] << " " + << index_of_vertex[fch->vertex(2)] << " " + << index_of_vertex[fch->vertex(3)] + << std::endl; + output << "3 " + << index_of_vertex[fch->vertex(0)] << " " + << index_of_vertex[fch->vertex(1)] << " " + << index_of_vertex[fch->vertex(3)] + << std::endl; + number_of_triangles += 4; + } + } + } + + os << "OFF \n" + << n << " " + << number_of_triangles << " 0\n" + << output.str(); + + return os; +} + +} //namespace CGAL + +#endif // CGAL_TRIANGULATION_IO_H diff --git a/Triangulation/include/CGAL/Regular_triangulation.h b/Triangulation/include/CGAL/Regular_triangulation.h index afc9455c1f8..39a4898ab5d 100644 --- a/Triangulation/include/CGAL/Regular_triangulation.h +++ b/Triangulation/include/CGAL/Regular_triangulation.h @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2014 INRIA Sophia-Antipolis (France). +// Copyright (c) 2014 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -15,7 +15,7 @@ // $URL$ // $Id$ // -// Author(s) : Samuel Hornus +// Author(s) : Clement Jamin #ifndef CGAL_REGULAR_TRIANGULATION_H #define CGAL_REGULAR_TRIANGULATION_H @@ -23,32 +23,1107 @@ #include #include #include +#include +#include + +#include namespace CGAL { template< typename RTTraits, typename TDS_ = Default > class Regular_triangulation -: public Triangulation::type, - Triangulation_vertex, - Triangulation_full_cell > - >::type > +: public Triangulation< + RTTraits, + typename Default::Get, + Triangulation_full_cell > + >::type > { - typedef typename Maximal_dimension::type - Maximal_dimension_; - typedef typename Default::Get, - Triangulation_full_cell > - >::type TDS; - typedef Triangulation Base; - typedef Regular_triangulation Self; + typedef typename RTTraits::Dimension Maximal_dimension_; + typedef typename Default::Get< + TDS_, + Triangulation_data_structure< + Maximal_dimension_, + Triangulation_vertex, + Triangulation_full_cell + > >::type TDS; + typedef Triangulation Base; + typedef Regular_triangulation Self; + + typedef typename RTTraits::Orientation_d Orientation_d; + typedef typename RTTraits::Construct_weighted_point_d Construct_weighted_point_d; + typedef typename RTTraits::Power_test_d Power_test_d; + typedef typename RTTraits::In_flat_power_test_d In_flat_power_test_d; + typedef typename RTTraits::Flat_orientation_d Flat_orientation_d; + typedef typename RTTraits::Construct_flat_orientation_d Construct_flat_orientation_d; + typedef typename RTTraits::In_flat_orientation_d In_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 Bare_point; + typedef typename RTTraits::Weighted_point 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: - typedef Maximal_dimension_ Maximal_dimension; -}; + + 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::insert_outside_convex_hull_1; + 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: + //*** Power_test_in_flat_d *** CJTODO: better name? + // Wrapper + struct Power_test_in_flat_d + { + boost::optional* fop; + Construct_flat_orientation_d cfo; + In_flat_power_test_d ifpt; + + Power_test_in_flat_d( + boost::optional& x, + Construct_flat_orientation_d const&y, + In_flat_power_test_d const&z) + : fop(&x), cfo(y), ifpt(z) {} + + template + 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 &preset_flat_orientation, + const Geom_traits &k = Geom_traits()) + : Base(dim, preset_flat_orientation, k) + { + } + + ~Regular_triangulation() {} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ACCESS + + // Not Documented + Power_test_in_flat_d power_test_in_flat_predicate() const + { + return Power_test_in_flat_d ( + flat_orientation_, + geom_traits().construct_flat_orientation_d_object(), + geom_traits().in_flat_power_test_d_object() + ); + } + + + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS + + 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 > + size_type insert(ForwardIterator start, ForwardIterator end) + { + size_type n = number_of_vertices(); + typedef std::vector WP_vec; + WP_vec points(start, end); + + typedef boost::function_property_map< + typename Geom_traits::Point_drop_weight_d, + Point, + Bare_point> Drop_weight_pmap; + typedef CGAL::Spatial_sort_traits_adapter_d< + typename Geom_traits::Base, Drop_weight_pmap> Search_traits_d; + Search_traits_d st_d(boost::make_function_property_map( + geom_traits().point_drop_weight_d_object())); + spatial_sort(points.begin(), points.end(), st_d); + + 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 &, + const Locate_type, + const Face &, + const Facet &, + const Full_cell_handle); + + Vertex_handle insert(const Weighted_point & p, + const 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, const 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 &, const Full_cell_handle, + const Vertex_handle only_if_this_vertex_is_in_the_cz = Vertex_handle()); + + Vertex_handle insert_if_in_star(const Weighted_point &, + const Vertex_handle, + const Locate_type, + const Face &, + const Facet &, + const Full_cell_handle); + + Vertex_handle insert_if_in_star( + const Weighted_point & p, const Vertex_handle star_center, + const 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, const Vertex_handle star_center, + const Vertex_handle hint) + { + CGAL_assertion( Vertex_handle() != hint ); + return insert_if_in_star(p, 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_test(const Weighted_point &, + Full_cell_const_handle, const OrientationPredicate &) const; + + template< typename OutputIterator > + Facet compute_conflict_zone(const Weighted_point &, const Full_cell_handle, OutputIterator) const; + + template < typename OrientationPredicate, typename PowerTestPredicate > + class Conflict_predicate + { + const Self & rt_; + const Weighted_point & p_; + OrientationPredicate ori_; + PowerTestPredicate power_test_; + int cur_dim_; + public: + Conflict_predicate( + const Self & rt, + const Weighted_point & p, + const OrientationPredicate & ori, + const PowerTestPredicate & power_test) + : rt_(rt), p_(p), ori_(ori), power_test_(power_test), cur_dim_(rt.current_dimension()) {} + + inline + bool operator()(Full_cell_const_handle s) const + { + bool ok; + if( ! rt_.is_infinite(s) ) + { + Oriented_side power_test = power_test_(rt_.points_begin(s), rt_.points_begin(s) + cur_dim_ + 1, p_); + if( ON_POSITIVE_SIDE == power_test ) + ok = true; + else if( ON_NEGATIVE_SIDE == power_test ) + ok = false; + else + ok = ON_POSITIVE_SIDE == rt_.perturbed_power_test(p_, s, ori_); + } + else + { + typedef typename Full_cell::Vertex_handle_const_iterator VHCI; + typedef Substitute_point_in_vertex_iterator 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; + +private: + + template + 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 + 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 + 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 + Conflict_pred_in_subspace; + typedef Conflict_predicate + Conflict_pred_in_fullspace; + typedef Conflict_traversal_predicate + Conflict_traversal_pred_in_subspace; + typedef Conflict_traversal_predicate + Conflict_traversal_pred_in_fullspace; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MEMBER VARIABLES + std::vector m_hidden_points; + +}; // class Regular_triangulation + + +// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = +// FUNCTIONS THAT ARE MEMBER METHODS: + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS + +template< typename RTTraits, typename TDS > +typename Regular_triangulation::Full_cell_handle +Regular_triangulation +::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( is_infinite(left) && left->neighbor(0)->index(left) == 0 ) // we are on the infinite right. + left = left->neighbor(0); + if( 0 == left->index(v) ) + left = left->neighbor(1); + CGAL_assertion( 1 == left->index(v) ); + Full_cell_handle right = left->neighbor(0); + if( ! is_infinite(right) ) + { + tds().associate_vertex_with_full_cell(left, 1, right->vertex(1)); + set_neighbors(left, 0, right->neighbor(0), right->mirror_index(0)); + } + else + { + tds().associate_vertex_with_full_cell(left, 1, left->vertex(0)); + tds().associate_vertex_with_full_cell(left, 0, infinite_vertex()); + set_neighbors(left, 0, left->neighbor(1), left->mirror_index(1)); + set_neighbors(left, 1, right->neighbor(1), right->mirror_index(1)); + } + 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 Simplices; + Simplices simps; + std::back_insert_iterator out(simps); + tds().incident_full_cells(v, out); + typedef std::set 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 Dark_vertex_base; + typedef Triangulation_full_cell< + Geom_traits, + internal::Triangulation::Dark_full_cell_data > Dark_full_cell_base; + typedef Triangulation_data_structure Dark_tds; + typedef Regular_triangulation 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(current_dimension(), flat_orientation_.get_ptr()) + : std::pair(std::numeric_limits::max(), NULL) ); + + Dark_s_handle dark_s; + Dark_v_handle dark_v; + typedef std::map 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()); + if( v_idx != 0 ) + { + // we must put the infinite vertex at index 0. + // OK, now with the new convention that the infinite vertex + // does not have to be at index 0, this is not necessary, + // but still, I prefer to keep this piece of code here. [-- Samuel Hornus] + (*it)->swap_vertices(0, v_idx); + // Now, we preserve the positive orientation of the full_cell + (*it)->swap_vertices(current_dimension() - 1, current_dimension()); + } + } + // 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 = 1; i <= current_dimension(); ++i ) + { + (*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_full_cells; + Dark_full_cells conflict_zone; + std::back_insert_iterator 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 > 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 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 RTTraits, typename TDS > +void +Regular_triangulation +::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() ) + { + // FIXME: infinite vertex is NOT at index 0 a priori. + Full_cell_handle s = infinite_vertex()->full_cell()->neighbor(0); + Orientation o = orientation(s); + CGAL_assertion( ZERO != o ); + if( NEGATIVE == o ) + reorient_full_cells(); + } +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INSERTIONS + +template< typename RTTraits, typename TDS > +typename Regular_triangulation::Vertex_handle +Regular_triangulation +::insert(const Weighted_point & p, const Locate_type lt, const Face & f, const Facet & ft, const 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)); + v->set_point(p); + return v; + break; + } + default: + if( 1 == current_dimension() ) + { + if( Base::OUTSIDE_CONVEX_HULL == lt ) + { + return insert_outside_convex_hull_1(p, s); + } + Vertex_handle v = tds().insert_in_full_cell(s); + v->set_point(p); + return v; + } + else + return insert_in_conflicting_cell(p, s); + break; + } +} + +template< typename RTTraits, typename TDS > +typename Regular_triangulation::Vertex_handle +Regular_triangulation +::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 ) + { + // FIXME: infinite vertex is NOT at index 0 a priori. + Full_cell_handle s = infinite_vertex()->full_cell()->neighbor(0); + Orientation o = orientation(s); + CGAL_assertion( ZERO != o ); + if( NEGATIVE == o ) + reorient_full_cells(); + } + return v; +} + +template< typename RTTraits, typename TDS > +typename Regular_triangulation::Vertex_handle +Regular_triangulation +::insert_if_in_star(const Weighted_point & p, + const Vertex_handle star_center, + const Locate_type lt, const Face & f, const Facet & ft, + const 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)); + v->set_point(p); + return v; + break; + } + default: + { + if( 1 == current_dimension() ) + { + if (s->has_vertex(star_center)) + { + if( Base::OUTSIDE_CONVEX_HULL == lt ) + { + return insert_outside_convex_hull_1(p, s); + } + Vertex_handle v = tds().insert_in_full_cell(s); + v->set_point(p); + return v; + } + } + else + return insert_in_conflicting_cell(p, s, star_center); + break; + } + } + + return Vertex_handle(); +} + +template< typename RTTraits, typename TDS > +typename Regular_triangulation::Vertex_handle +Regular_triangulation +::insert_in_conflicting_cell(const Weighted_point & p, + const Full_cell_handle s, + const Vertex_handle only_if_this_vertex_is_in_the_cz) +{ + typedef std::vector Full_cell_h_vector; + + bool is_in_conflict = true; + if( current_dimension() < maximal_dimension() ) + { + Conflict_pred_in_subspace c( + *this, p, + coaffine_orientation_predicate(), + power_test_in_flat_predicate()); + is_in_conflict = c(s); + } + else + { + Orientation_d ori = geom_traits().orientation_d_object(); + Power_test_d side = geom_traits().power_test_d_object(); + Conflict_pred_in_fullspace c(*this, p, ori, side); + is_in_conflict = c(s); + } + + // If p is not in conflict with s, then p is hidden + // => we don't insert it + if (!is_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 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 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 RTTraits, typename TDS > +template< typename OrientationPred > +Oriented_side +Regular_triangulation +::perturbed_power_test(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 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(*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, + 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 RTTraits, typename TDS > +bool +Regular_triangulation +::is_in_conflict(const Weighted_point & p, Full_cell_const_handle s) const +{ + CGAL_precondition( 2 <= current_dimension() ); + if( current_dimension() < maximal_dimension() ) + { + Conflict_pred_in_subspace c( + *this, p, + coaffine_orientation_predicate(), + power_test_in_flat_predicate()); + return c(s); + } + else + { + Orientation_d ori = geom_traits().orientation_d_object(); + Power_test_d side = geom_traits().power_test_d_object(); + Conflict_pred_in_fullspace c(*this, p, ori, side); + return c(s); + } +} + +template< typename RTTraits, typename TDS > +template< typename OutputIterator > +typename Regular_triangulation::Facet +Regular_triangulation +::compute_conflict_zone(const Weighted_point & p, const Full_cell_handle s, OutputIterator out) const +{ + CGAL_precondition( 2 <= current_dimension() ); + if( current_dimension() < maximal_dimension() ) + { + Conflict_pred_in_subspace c( + *this, p, + coaffine_orientation_predicate(), + power_test_in_flat_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_test_d side = geom_traits().power_test_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 RTTraits, typename TDS > +bool +Regular_triangulation +::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_test_d side = + geom_traits().power_test_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 +#endif //CGAL_REGULAR_TRIANGULATION_H diff --git a/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h b/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h new file mode 100644 index 00000000000..3f5a9630a85 --- /dev/null +++ b/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h @@ -0,0 +1,254 @@ +// 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_EUCLIDEAN_TRAITS_H +#define CGAL_REGULAR_TRIANGULATION_EUCLIDEAN_TRAITS_H + +#include +#include +#include +#include +#include + +#include + +namespace CGAL { + +template < class K, class Weight = typename K::RT > +class Regular_triangulation_euclidean_traits + : public K +{ +public: + typedef K Base; + typedef Regular_triangulation_euclidean_traits Self; + + // Types from K + + typedef K Kernel; + typedef typename K::Dimension Dimension; + typedef typename K::FT FT; + typedef typename K::Point_d Bare_point; + typedef typename K::Weighted_point_d Weighted_point; + typedef Weighted_point Weighted_point_d; + typedef Weighted_point Point_d; + + typedef typename K::Construct_weighted_point_d Construct_weighted_point_d; + typedef typename K::Power_test_d Power_test_d; + typedef typename K::In_flat_power_test_d In_flat_power_test_d; + typedef typename K::Flat_orientation_d Flat_orientation_d; + typedef typename K::Point_drop_weight_d Point_drop_weight_d; + + //============================================================================= + // Custom types + //============================================================================= + + class Orientation_d + { + const K &m_kernel; + + public: + typedef Orientation result_type; + + Orientation_d(const K &kernel) + : m_kernel(kernel) {} + + template + result_type operator()(ForwardIterator start, ForwardIterator end) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.orientation_d_object() ( + boost::make_transform_iterator(start, pdw), + boost::make_transform_iterator(end, pdw) + ); + } + }; + + //============================================================================= + + class Construct_flat_orientation_d + { + const K &m_kernel; + + public: + typedef Flat_orientation_d result_type; + + Construct_flat_orientation_d(const K &kernel) + : m_kernel(kernel) {} + + template + result_type operator()(ForwardIterator start, ForwardIterator end) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.construct_flat_orientation_d_object() ( + boost::make_transform_iterator(start, pdw), + boost::make_transform_iterator(end, pdw) + ); + } + }; + + + //============================================================================= + + class In_flat_orientation_d + { + const K &m_kernel; + + public: + typedef Orientation result_type; + + In_flat_orientation_d(const K &kernel) + : m_kernel(kernel) {} + + template + result_type operator()(Flat_orientation_d orient, + ForwardIterator start, ForwardIterator end) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.in_flat_orientation_d_object() ( + orient, + boost::make_transform_iterator(start, pdw), + boost::make_transform_iterator(end, pdw) + ); + } + }; + + //============================================================================= + + class Contained_in_affine_hull_d + { + const K &m_kernel; + + public: + typedef bool result_type; + + Contained_in_affine_hull_d(const K &kernel) + : m_kernel(kernel) {} + + template + result_type operator()(ForwardIterator start, ForwardIterator end, + const Weighted_point_d & p) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.contained_in_affine_hull_d_object() ( + boost::make_transform_iterator(start, pdw), + boost::make_transform_iterator(end, pdw), + pdw(p) + ); + } + }; + + //============================================================================= + + class Compare_lexicographically_d + { + const K &m_kernel; + + public: + typedef Comparison_result result_type; + + Compare_lexicographically_d(const K &kernel) + : m_kernel(kernel) {} + + result_type operator()( + const Weighted_point_d & p, const Weighted_point_d & q) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.compare_lexicographically_d_object()(pdw(p), pdw(q)); + } + }; + + //============================================================================= + + class Compute_coordinate_d + { + const K &m_kernel; + + public: + typedef FT result_type; + + Compute_coordinate_d(const K &kernel) + : m_kernel(kernel) {} + + result_type operator()( + const Weighted_point_d & p, const int i) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + m_kernel.compute_coordinate_d_object()(pdw(p), i); + return m_kernel.compute_coordinate_d_object()(pdw(p), i); + } + }; + + //============================================================================= + + class Point_dimension_d + { + const K &m_kernel; + + public: + typedef int result_type; + + Point_dimension_d(const K &kernel) + : m_kernel(kernel) {} + + result_type operator()( + const Weighted_point_d & p) const + { + Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object(); + return m_kernel.point_dimension_d_object()(pdw(p)); + } + }; + + //============================================================================= + // Object creation + //============================================================================= + + Contained_in_affine_hull_d contained_in_affine_hull_d_object() const + { + return Contained_in_affine_hull_d(*this); + } + Orientation_d orientation_d_object() const + { + return Orientation_d(*this); + } + Construct_flat_orientation_d construct_flat_orientation_d_object() const + { + return Construct_flat_orientation_d(*this); + } + In_flat_orientation_d in_flat_orientation_d_object() const + { + return In_flat_orientation_d(*this); + } + Compare_lexicographically_d compare_lexicographically_d_object() const + { + return Compare_lexicographically_d(*this); + } + Compute_coordinate_d compute_coordinate_d_object() const + { + return Compute_coordinate_d(*this); + } + Point_dimension_d point_dimension_d_object() const + { + return Point_dimension_d(*this); + } +}; + + +} //namespace CGAL + +#endif // CGAL_REGULAR_TRIANGULATION_EUCLIDEAN_TRAITS_H diff --git a/Triangulation/include/CGAL/Triangulation.h b/Triangulation/include/CGAL/Triangulation.h index a95a025e426..db57b817ddc 100644 --- a/Triangulation/include/CGAL/Triangulation.h +++ b/Triangulation/include/CGAL/Triangulation.h @@ -226,7 +226,35 @@ public: { return tds().index_of_covertex(f); } + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - UTILITIES + + // A co-dimension 2 sub-simplex. called a Rotor because we can rotate + // the two "covertices" around the sub-simplex. Useful for traversing the + // boundary of a hole. NOT DOCUMENTED + typedef cpp11::tuple Rotor; + // Commented out because it was causing "internal compiler error" in MSVC + /*Full_cell_handle full_cell(const Rotor & r) const // NOT DOCUMENTED + { + return cpp11::get<0>(r); + } + int index_of_covertex(const Rotor & r) const // NOT DOCUMENTED + { + return cpp11::get<1>(r); + } + int index_of_second_covertex(const Rotor & r) const // NOT DOCUMENTED + { + return cpp11::get<2>(r); + }*/ + Rotor rotate_rotor(Rotor & r) // NOT DOCUMENTED... + { + int opposite = cpp11::get<0>(r)->mirror_index(cpp11::get<1>(r)); + Full_cell_handle s = cpp11::get<0>(r)->neighbor(cpp11::get<1>(r)); + int new_second = s->index(cpp11::get<0>(r)->vertex(cpp11::get<2>(r))); + return Rotor(s, new_second, opposite); + } + // - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS Triangulation(int dim, const Geom_traits k = Geom_traits()) @@ -696,7 +724,7 @@ public: Full_cell_handle n = s->neighbor(i); if( ! t_.is_infinite(n) ) return false; - int inf_v_index = n->index(infinite_vertex()); + int inf_v_index = n->index(t_.infinite_vertex()); n->vertex(inf_v_index)->set_point(p_); bool ok = (POSITIVE == ori_(t_.points_begin(n), t_.points_begin(n) + cur_dim_ + 1)); return ok; @@ -706,6 +734,28 @@ public: // make sure all full_cells have positive orientation void reorient_full_cells(); +protected: + // This is used in the |remove(v)| member function to manage sets of Full_cell_handles + template< typename FCH > + struct Full_cell_set : public std::vector + { + typedef std::vector Base_set; + using Base_set::begin; + using Base_set::end; + void make_searchable() + { // sort the full cell handles + std::sort(begin(), end()); + } + bool contains(const FCH & fch) const + { + return std::binary_search(begin(), end(), fch); + } + bool contains_1st_and_not_2nd(const FCH & fst, const FCH & snd) const + { + return ( ! contains(snd) ) && ( contains(fst) ); + } + }; + }; // Triangulation<...> // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = diff --git a/Triangulation/include/CGAL/Triangulation_data_structure.h b/Triangulation/include/CGAL/Triangulation_data_structure.h index 3a2681f9658..bdbc16e7ad0 100644 --- a/Triangulation/include/CGAL/Triangulation_data_structure.h +++ b/Triangulation/include/CGAL/Triangulation_data_structure.h @@ -417,7 +417,6 @@ private: void clear_visited_marks(Full_cell_handle) const; // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - DANGEROUS UPDATE OPERATIONS - // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - DANGEROUS UPDATE OPERATIONS private: @@ -1523,7 +1522,9 @@ operator>>(std::istream & is, Triangulation_data_structure & tr) // - the neighbors of each full_cell by their index in the preceding list { typedef Triangulation_data_structure TDS; - typedef typename TDS::Vertex_handle Vertex_handle; + typedef typename TDS::Full_cell_handle Full_cell_handle; + typedef typename TDS::Full_cell_iterator Full_cell_iterator; + typedef typename TDS::Vertex_handle Vertex_handle; // read current dimension and number of vertices std::size_t n; @@ -1573,8 +1574,10 @@ operator<<(std::ostream & os, const Triangulation_data_structure // - the neighbors of each full_cell by their index in the preceding list { typedef Triangulation_data_structure TDS; - typedef typename TDS::Vertex_const_handle Vertex_handle; - typedef typename TDS::Vertex_const_iterator Vertex_iterator; + typedef typename TDS::Full_cell_const_handle Full_cell_handle; + typedef typename TDS::Full_cell_const_iterator Full_cell_iterator; + typedef typename TDS::Vertex_const_handle Vertex_handle; + typedef typename TDS::Vertex_const_iterator Vertex_iterator; // outputs dimension and number of vertices std::size_t n = tr.number_of_vertices(); @@ -1594,7 +1597,7 @@ operator<<(std::ostream & os, const Triangulation_data_structure int i = 0; for( Vertex_iterator it = tr.vertices_begin(); it != tr.vertices_end(); ++it, ++i ) { - os << *it; // write the vertex + os << *it << std::endl; // write the vertex index_of_vertex[it] = i; } CGAL_assertion( (std::size_t) i == n ); diff --git a/Triangulation/include/CGAL/Triangulation_ds_vertex.h b/Triangulation/include/CGAL/Triangulation_ds_vertex.h index b141936d789..381b97e1aa2 100644 --- a/Triangulation/include/CGAL/Triangulation_ds_vertex.h +++ b/Triangulation/include/CGAL/Triangulation_ds_vertex.h @@ -61,7 +61,6 @@ public: /// Set 's' as an incident full_cell void set_full_cell(Full_cell_handle s) /* Concept */ { - CGAL_precondition( Full_cell_handle() != s ); full_cell_ = s; } diff --git a/Triangulation/include/CGAL/internal/Triangulation/utilities.h b/Triangulation/include/CGAL/internal/Triangulation/utilities.h index 00092cb34c8..a1ffc7750c1 100644 --- a/Triangulation/include/CGAL/internal/Triangulation/utilities.h +++ b/Triangulation/include/CGAL/internal/Triangulation/utilities.h @@ -98,7 +98,7 @@ public: template< class T > struct Compare_points_for_perturbation { - typedef typename T::Point_d Point; + typedef typename T::Geom_traits::Point_d Point; const T & t_; @@ -119,8 +119,8 @@ public: template< class T > struct Point_from_pointer { - typedef const typename T::Point_d * argument_type; - typedef const typename T::Point_d result_type; + typedef const typename T::Geom_traits::Point_d * argument_type; + typedef const typename T::Geom_traits::Point_d result_type; result_type & operator()(argument_type & x) const { return (*x); diff --git a/Triangulation/test/Triangulation/CMakeLists.txt b/Triangulation/test/Triangulation/CMakeLists.txt index 89abb618ce7..beede01c0df 100644 --- a/Triangulation/test/Triangulation/CMakeLists.txt +++ b/Triangulation/test/Triangulation/CMakeLists.txt @@ -1,7 +1,6 @@ # Created by the script cgal_create_cmake_script # This is the CMake script for compiling a CGAL application. - project( Triangulation_test ) cmake_minimum_required(VERSION 2.6.2) @@ -27,18 +26,18 @@ if ( CGAL_FOUND ) include_directories (BEFORE "../../include") include_directories (BEFORE "include") + create_single_source_cgal_program( "test_triangulation.cpp" ) create_single_source_cgal_program( "test_delaunay.cpp" ) + create_single_source_cgal_program( "test_regular.cpp" ) create_single_source_cgal_program( "test_tds.cpp" ) create_single_source_cgal_program( "test_torture.cpp" ) - create_single_source_cgal_program( "test_triangulation.cpp" ) + create_single_source_cgal_program( "test_insert_if_in_star.cpp" ) else() message(STATUS "NOTICE: Some of the executables in this directory need Eigen 3.1 (or greater) and will not be compiled.") endif() else() - message(STATUS "This program requires the CGAL library, and will not be compiled.") - endif() diff --git a/Triangulation/test/Triangulation/test_delaunay.cpp b/Triangulation/test/Triangulation/test_delaunay.cpp index 9e7076f2194..3b98f66660d 100644 --- a/Triangulation/test/Triangulation/test_delaunay.cpp +++ b/Triangulation/test/Triangulation/test_delaunay.cpp @@ -24,79 +24,80 @@ void test(const int d, const string & type, const int N) typedef CGAL::Random_points_in_cube_d Random_points_iterator; - DC pc(d); + DC dt(d); cerr << "\nBuilding Delaunay triangulation of (" << type << d << ") dimension with " << N << " points"; - assert(pc.empty()); + assert(dt.empty()); vector points; - CGAL::Random rng; - Random_points_iterator rand_it(d, 2.0, rng); + //CGAL::Random rng; + //Random_points_iterator rand_it(d, 2.0, rng); //CGAL::cpp11::copy_n(rand_it, N, back_inserter(points)); - vector coords(d); + srand(10); for( int i = 0; i < N; ++i ) { + vector coords(d); for( int j = 0; j < d; ++j ) - coords[j] = rand() % 100000; + coords[j] = static_cast(rand() % 100000)/10000; points.push_back(Point(d, coords.begin(), coords.end())); } - pc.insert(points.begin(), points.end()); + dt.insert(points.begin(), points.end()); cerr << "\nChecking topology and geometry..."; - assert( pc.is_valid() ); + assert( dt.is_valid() ); cerr << "\nTraversing finite full_cells... "; size_t nbfs(0), nbis(0); - Finite_full_cell_const_iterator fsit = pc.finite_full_cells_begin(); - while( fsit != pc.finite_full_cells_end() ) + Finite_full_cell_const_iterator fsit = dt.finite_full_cells_begin(); + while( fsit != dt.finite_full_cells_end() ) ++fsit, ++nbfs; cerr << nbfs << " + "; vector infinite_full_cells; - pc.tds().incident_full_cells(pc.infinite_vertex(), back_inserter(infinite_full_cells)); + dt.tds().incident_full_cells(dt.infinite_vertex(), back_inserter(infinite_full_cells)); nbis = infinite_full_cells.size(); cerr << nbis << " = " << (nbis+nbfs) - << " = " << pc.number_of_full_cells(); - cerr << "\nThe triangulation has current dimension " << pc.current_dimension(); - CGAL_assertion( pc.number_of_full_cells() == nbis+nbfs); + << " = " << dt.number_of_full_cells(); + cerr << "\nThe triangulation has current dimension " << dt.current_dimension(); + CGAL_assertion( dt.number_of_full_cells() == nbis+nbfs); cerr << "\nTraversing finite vertices... "; size_t nbfv(0); - Finite_vertex_iterator fvit = pc.finite_vertices_begin(); - while( fvit != pc.finite_vertices_end() ) + Finite_vertex_iterator fvit = dt.finite_vertices_begin(); + while( fvit != dt.finite_vertices_end() ) ++fvit, ++nbfv; cerr << nbfv < 1 ) + if( dt.maximal_dimension() > 1 ) { typedef vector Faces; Faces edges; back_insert_iterator out(edges); - pc.tds().incident_faces(pc.infinite_vertex(), 1, out); + dt.tds().incident_faces(dt.infinite_vertex(), 1, out); cout << "\nThere are " << edges.size() << " vertices on the convex hull."; edges.clear(); } - else // pc.maximal_dimension() == 1 + else // dt.maximal_dimension() == 1 { typedef vector Cells; Cells cells; back_insert_iterator out(cells); - pc.tds().incident_full_cells(pc.infinite_vertex(), out); + dt.tds().incident_full_cells(dt.infinite_vertex(), out); cout << "\nThere are " << cells.size() << " vertices on the convex hull."; cells.clear(); } // Remove all ! - cerr << "\nBefore removal: " << pc.number_of_vertices() << " vertices. After: "; + cerr << "\nBefore removal: " << dt.number_of_vertices() << " vertices. After: "; random_shuffle(points.begin(), points.end()); - pc.remove(points.begin(), points.end()); - assert( pc.is_valid() ); - cerr << pc.number_of_vertices() << " vertices."; - // assert( pc.empty() ); NOT YET ! + dt.remove(points.begin(), points.end()); + assert( dt.is_valid() ); + cerr << dt.number_of_vertices() << " vertices."; + // assert( dt.empty() ); NOT YET ! // CLEAR - pc.clear(); - assert( -1 == pc.current_dimension() ); - assert( pc.empty() ); - assert( pc.is_valid() ); + dt.clear(); + assert( -1 == dt.current_dimension() ); + assert( dt.empty() ); + assert( dt.is_valid() ); } template< int D > @@ -112,14 +113,14 @@ void go(const int N) int main(int argc, char **argv) { srand(static_cast(time(NULL))); - int N = 100; + int N = 10; if( argc > 1 ) N = atoi(argv[1]); - go<5>(N); - go<4>(N); - go<3>(N); - go<2>(N); - go<1>(N); + //go<5>(N); + //go<4>(N); + //go<3>(N); + go<2>(N); + //go<1>(N); cerr << endl; return 0; diff --git a/Triangulation/test/Triangulation/test_insert_if_in_star.cpp b/Triangulation/test/Triangulation/test_insert_if_in_star.cpp new file mode 100644 index 00000000000..3621c034fbc --- /dev/null +++ b/Triangulation/test/Triangulation/test_insert_if_in_star.cpp @@ -0,0 +1,94 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace std; + +template +void test(const int d, const string & type, const int N) +{ + typedef typename RTri::Vertex_handle Vertex_handle; + typedef typename RTri::Point Point; + typedef typename RTri::Bare_point Bare_point; + + typedef CGAL::Random_points_in_cube_d Random_points_iterator; + + RTri rt(d); + RTri rt_star_only(d); + cerr << "\nBuilding Regular triangulation of (" << type << d + << ") dimension with " << N << " points\n"; + assert(rt.empty()); + assert(rt_star_only.empty()); + + srand(static_cast(time(NULL))); + + // Insert first point (0, 0...) + vector coords(d); + for( int j = 0; j < d; ++j ) + coords[j] = 0; + + Point p = Point( + Bare_point(d, coords.begin(), coords.end()), + static_cast(rand() % 10000)/100000); + + rt.insert(p); + Vertex_handle first_vertex = rt_star_only.insert(p); + + // Insert the other points + for( int i = 1 ; i < N ; ++i ) + { + for( int j = 0; j < d; ++j ) + coords[j] = 10.*(rand() % RAND_MAX)/RAND_MAX - 5.; + + p = Point( + Bare_point(d, coords.begin(), coords.end()), + static_cast(rand() % 10000)/1000000); + + rt.insert(p); + rt_star_only.insert_if_in_star(p, first_vertex); + } + + cerr << "\nChecking topology and geometry..." + << (rt.is_valid(true) ? "OK.\n" : "Error.\n"); + + cerr << "\nThe triangulation using 'insert' has current dimension " << rt.current_dimension() + << " and " << rt.number_of_full_cells() << " full cells\n"; + + cerr << "\nThe triangulation using 'insert_if_in_star' has current dimension " << rt.current_dimension() + << " and " << rt_star_only.number_of_full_cells() << " full cells\n"; + + // Export + if (d <= 3) + { + std::ofstream off_stream_all("data/test_insert_all.off"); + CGAL::export_triangulation_to_off(off_stream_all, rt); + std::ofstream off_stream_star_only("data/test_insert_if_in_star.off"); + CGAL::export_triangulation_to_off(off_stream_star_only, rt_star_only); + } +} + +template< int D > +void go(const int N) +{ + //typedef CGAL::Epick_d FK; + typedef CGAL::Epick_d > FK; + typedef CGAL::Regular_triangulation_euclidean_traits Traits; + typedef CGAL::Regular_triangulation Triangulation; + //test(D, "dynamic", N); + test(D, "static", N); +} + +int main(int argc, char **argv) +{ + go<2>(100); + return 0; +} diff --git a/Triangulation/test/Triangulation/test_regular.cpp b/Triangulation/test/Triangulation/test_regular.cpp new file mode 100644 index 00000000000..d1a6c568777 --- /dev/null +++ b/Triangulation/test/Triangulation/test_regular.cpp @@ -0,0 +1,132 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace std; + +template +void test(const int d, const string & type, const int N) +{ + typedef typename RTri::Full_cell_handle Full_cell_handle; + typedef typename RTri::Face Face; + typedef typename RTri::Point Point; + typedef typename RTri::Bare_point Bare_point; + typedef typename RTri::Finite_full_cell_const_iterator Finite_full_cell_const_iterator; + typedef typename RTri::Finite_vertex_iterator Finite_vertex_iterator; + + typedef CGAL::Random_points_in_cube_d Random_points_iterator; + + RTri rt(d); + cerr << "\nBuilding Regular triangulation of (" << type << d << ") dimension with " << N << " points"; + assert(rt.empty()); + + vector points; + //CGAL::Random rng; + //Random_points_iterator rand_it(d, 2.0, rng); // CJTODO: unused + + srand(10); + for( int i = 0; i < N; ++i ) + { + vector coords(d); + for( int j = 0; j < d; ++j ) + coords[j] = static_cast(rand() % 100000)/10000; + points.push_back(Point( + Bare_point(d, coords.begin(), coords.end()), + /*static_cast(rand() % 100000)/100000*/static_cast(i)/20 + )); + } + rt.insert(points.begin(), points.end()); + cerr << "\nChecking topology and geometry..."; + assert( rt.is_valid(true) ); + + cerr << "\nTraversing finite full_cells... "; + size_t nbfs(0), nbis(0); + Finite_full_cell_const_iterator fsit = rt.finite_full_cells_begin(); + while( fsit != rt.finite_full_cells_end() ) + ++fsit, ++nbfs; + cerr << nbfs << " + "; + vector infinite_full_cells; + rt.tds().incident_full_cells(rt.infinite_vertex(), back_inserter(infinite_full_cells)); + nbis = infinite_full_cells.size(); + cerr << nbis << " = " << (nbis+nbfs) + << " = " << rt.number_of_full_cells(); + cerr << "\nThe triangulation has current dimension " << rt.current_dimension(); + CGAL_assertion( rt.number_of_full_cells() == nbis+nbfs); + + cerr << "\nTraversing finite vertices... "; + size_t nbfv(0); + Finite_vertex_iterator fvit = rt.finite_vertices_begin(); + while( fvit != rt.finite_vertices_end() ) + ++fvit, ++nbfv; + cerr << nbfv < 1 ) + { + typedef vector Faces; + Faces edges; + back_insert_iterator out(edges); + rt.tds().incident_faces(rt.infinite_vertex(), 1, out); + cout << "\nThere are " << edges.size() << " vertices on the convex hull."; + edges.clear(); + } + else // rt.maximal_dimension() == 1 + { + typedef vector Cells; + Cells cells; + back_insert_iterator out(cells); + rt.tds().incident_full_cells(rt.infinite_vertex(), out); + cout << "\nThere are " << cells.size() << " vertices on the convex hull."; + cells.clear(); + } + + // Remove all ! + cerr << "\nBefore removal: " << rt.number_of_vertices() << " vertices. After: "; + random_shuffle(points.begin(), points.end()); + rt.remove(points.begin(), points.end()); + assert( rt.is_valid() ); + //std::cerr << ((rt.is_valid(true)) ? "VALID!" : "NOT VALID :(") << std::endl; + cerr << rt.number_of_vertices() << " vertices."; + // assert( rt.empty() ); NOT YET ! + // CLEAR + rt.clear(); + assert( -1 == rt.current_dimension() ); + assert( rt.empty() ); + assert( rt.is_valid() ); + //std::cerr << ((rt.is_valid(true)) ? "VALID!" : "NOT VALID :(") << std::endl; +} + +template< int D > +void go(const int N) +{ + //typedef CGAL::Epick_d FK; + typedef CGAL::Epick_d > FK; + typedef CGAL::Regular_triangulation_euclidean_traits Traits; + typedef CGAL::Regular_triangulation Triangulation; + //test(D, "dynamic", N); + test(D, "static", N); +} + +int main(int argc, char **argv) +{ + srand(static_cast(time(NULL))); + int N = 10; + if( argc > 1 ) + N = atoi(argv[1]); + //go<5>(N); + //go<4>(N); + go<3>(N); + go<2>(N); + + cerr << endl; + return 0; +} diff --git a/Triangulation/test/Triangulation/test_triangulation.cpp b/Triangulation/test/Triangulation/test_triangulation.cpp index 15b83c281de..f6230bf4250 100644 --- a/Triangulation/test/Triangulation/test_triangulation.cpp +++ b/Triangulation/test/Triangulation/test_triangulation.cpp @@ -114,10 +114,10 @@ int main(int argc, char **argv) int N = 1000; if( argc > 1 ) N = atoi(argv[1]); - go<5>(N); - go<3>(N); + //go<5>(N); + //go<3>(N); go<2>(N); - go<1>(N); + //go<1>(N); cerr << std::endl; return 0; diff --git a/Triangulation_2/applications/Triangulation_2/data/points.cin b/Triangulation_2/applications/Triangulation_2/data/points.cin new file mode 100644 index 00000000000..5e78e2349bc --- /dev/null +++ b/Triangulation_2/applications/Triangulation_2/data/points.cin @@ -0,0 +1,19 @@ +0 0 6.28953 +-2.85086 -0.471442 6.12896 +1.90972 0.101219 0.988689 +0.637771 2.59367 5.80372 +2.22209 0.903198 2.19478 +-0.487202 -2.71506 4.90996 +1.1193 -1.91787 2.99626 +1.54714 0.109831 0 +0.44556 -2.73047 4.48142 +0.427936 1.28495 6.23624 +-2.67212 0.766674 5.29623 +1.5763 -1.59828 2.58905 +-0.476603 2.2546 6.04797 +1.57172 -0.514711 6.11405 +1.84528 2.10139 5.53936 +-2.99827 -0.101677 5.92246 +-0.482122 -2.39584 4.44264 +-2.25558 -1.492 6.23448 +0.128475 -1.75125 3.18916 \ No newline at end of file diff --git a/Triangulation_2/applications/Triangulation_2/points_to_RT2_to_off.cpp b/Triangulation_2/applications/Triangulation_2/points_to_RT2_to_off.cpp new file mode 100644 index 00000000000..c2989d5a902 --- /dev/null +++ b/Triangulation_2/applications/Triangulation_2/points_to_RT2_to_off.cpp @@ -0,0 +1,28 @@ +#include +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Regular_triangulation_filtered_traits_2 Traits; +typedef CGAL::Regular_triangulation_2 Regular_triangulation; + +int main() +{ + std::ifstream in("data/points.cin"); + + Regular_triangulation::Weighted_point wp; + std::vector wpoints; + + while(in >> wp) + wpoints.push_back(wp); + + Regular_triangulation rt(wpoints.begin(), wpoints.end()); + CGAL_assertion(rt.is_valid(true)); + std::ofstream off_stream("data/rt2.off"); + CGAL::export_triangulation_2_to_off(off_stream, rt); + return 0; +} diff --git a/Triangulation_2/include/CGAL/IO/Triangulation_off_ostream_2.h b/Triangulation_2/include/CGAL/IO/Triangulation_off_ostream_2.h new file mode 100644 index 00000000000..a4cd70e8e0d --- /dev/null +++ b/Triangulation_2/include/CGAL/IO/Triangulation_off_ostream_2.h @@ -0,0 +1,79 @@ +// 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 Lesser 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_TRIANGULATION_OFF_OSTREAM_2_H +#define CGAL_TRIANGULATION_OFF_OSTREAM_2_H + +#include +#include +#include + +namespace CGAL { + +template < class GT, class TDS > +std::ostream & +export_triangulation_2_to_off(std::ostream & os, + const Triangulation_2 & tr) +{ + typedef Triangulation_2 Tr; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Vertex_iterator Vertex_iterator; + typedef typename Tr::Finite_vertices_iterator Finite_vertex_iterator; + typedef typename Tr::Finite_faces_iterator Finite_faces_iterator; + + size_t n = tr.number_of_vertices(); + + std::stringstream output; + + // write the vertices + std::map index_of_vertex; + int i = 0; + for(Finite_vertex_iterator it = tr.finite_vertices_begin(); + it != tr.finite_vertices_end(); ++it, ++i) + { + output << it->point().x() << " " << it->point().y() << " 0" << std::endl; + index_of_vertex[it.base()] = i; + } + CGAL_assertion( i == n ); + + size_t number_of_triangles = 0; + + for (Finite_faces_iterator fit = tr.finite_faces_begin() ; + fit != tr.finite_faces_end() ; ++fit) + { + output << "3 " + << index_of_vertex[fit->vertex(0)] << " " + << index_of_vertex[fit->vertex(1)] << " " + << index_of_vertex[fit->vertex(2)] + << std::endl; + ++number_of_triangles; + } + + os << "OFF \n" + << n << " " + << number_of_triangles << " 0\n" + << output.str(); + + return os; +} + +} //namespace CGAL + +#endif // CGAL_TRIANGULATION_OFF_OSTREAM_2_H \ No newline at end of file diff --git a/Triangulation_3/applications/Triangulation_3/data/points - extreme weights.cin b/Triangulation_3/applications/Triangulation_3/data/points - extreme weights.cin new file mode 100644 index 00000000000..10c0791ad4c --- /dev/null +++ b/Triangulation_3/applications/Triangulation_3/data/points - extreme weights.cin @@ -0,0 +1,10 @@ +0.0071 1.6899 2.521 0 +0.3272 1.3694 3.15 0.05 +1.3697 1.8296 2.654 0.1 +-10.6722 0.3012 0.1548 1000.15 +1.1726 0.1899 0.3658 0.2 +0.4374 20.8541 1.45894 2000.25 +2.5923 0.1904 0.6971 0.3 +10.3083 2.5462 1.3658 1000.35 +1.4981 1.3929 2.949 0.4 +2.1304 2.055 0.6597455 1.45 \ No newline at end of file diff --git a/Triangulation_3/applications/Triangulation_3/data/points - no weights.cin b/Triangulation_3/applications/Triangulation_3/data/points - no weights.cin new file mode 100644 index 00000000000..9dcea4973dc --- /dev/null +++ b/Triangulation_3/applications/Triangulation_3/data/points - no weights.cin @@ -0,0 +1,10 @@ +0.0071 1.6899 2.521 0 +0.3272 1.3694 3.15 0 +1.3697 1.8296 2.654 0 +-10.6722 0.3012 0.1548 0 +1.1726 0.1899 0.3658 0 +0.4374 20.8541 1.45894 0 +2.5923 0.1904 0.6971 0 +10.3083 2.5462 1.3658 0 +1.4981 1.3929 2.949 0 +2.1304 2.055 0.6597455 0 \ No newline at end of file diff --git a/Triangulation_3/applications/Triangulation_3/data/points.cin b/Triangulation_3/applications/Triangulation_3/data/points.cin new file mode 100644 index 00000000000..9dcea4973dc --- /dev/null +++ b/Triangulation_3/applications/Triangulation_3/data/points.cin @@ -0,0 +1,10 @@ +0.0071 1.6899 2.521 0 +0.3272 1.3694 3.15 0 +1.3697 1.8296 2.654 0 +-10.6722 0.3012 0.1548 0 +1.1726 0.1899 0.3658 0 +0.4374 20.8541 1.45894 0 +2.5923 0.1904 0.6971 0 +10.3083 2.5462 1.3658 0 +1.4981 1.3929 2.949 0 +2.1304 2.055 0.6597455 0 \ No newline at end of file diff --git a/Triangulation_3/applications/Triangulation_3/points_to_RT3_to_off.cpp b/Triangulation_3/applications/Triangulation_3/points_to_RT3_to_off.cpp new file mode 100644 index 00000000000..4dffb3ea265 --- /dev/null +++ b/Triangulation_3/applications/Triangulation_3/points_to_RT3_to_off.cpp @@ -0,0 +1,26 @@ +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Regular_triangulation_euclidean_traits_3 Traits; +typedef CGAL::Regular_triangulation_3 Regular_triangulation; + +int main() +{ + std::ifstream in("data/points.cin"); + + Regular_triangulation::Weighted_point wp; + std::vector wpoints; + + while(in >> wp) + wpoints.push_back(wp); + + Regular_triangulation rt(wpoints.begin(), wpoints.end()); + std::ofstream off_stream("data/rt3.off"); + CGAL::export_triangulation_3_to_off(off_stream, rt); + return 0; +} diff --git a/Triangulation_3/include/CGAL/IO/Triangulation_off_ostream_3.h b/Triangulation_3/include/CGAL/IO/Triangulation_off_ostream_3.h new file mode 100644 index 00000000000..b6287f71d52 --- /dev/null +++ b/Triangulation_3/include/CGAL/IO/Triangulation_off_ostream_3.h @@ -0,0 +1,119 @@ +// 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 Lesser 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_TRIANGULATION_OFF_OSTREAM_3_H +#define CGAL_TRIANGULATION_OFF_OSTREAM_3_H + +#include +#include +#include + +namespace CGAL { + +template < class GT, class TDS > +std::ostream & +export_triangulation_3_to_off(std::ostream & os, + const Triangulation_3 & tr, + bool export_surface_only = false) +{ + typedef Triangulation_3 Tr; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Vertex_iterator Vertex_iterator; + typedef typename Tr::Finite_vertices_iterator Finite_vertex_iterator; + typedef typename Tr::All_cells_iterator Cells_iterator; + typedef typename Tr::Finite_cells_iterator Finite_cells_iterator; + + size_t n = tr.number_of_vertices(); + + std::stringstream output; + + // write the vertices + std::map index_of_vertex; + int i = 0; + for(Finite_vertex_iterator it = tr.finite_vertices_begin(); + it != tr.finite_vertices_end(); ++it, ++i) + { + output << it->point().x() << " " + << it->point().y() << " " + << it->point().z() << std::endl; + index_of_vertex[it.base()] = i; + } + CGAL_assertion( i == n ); + + size_t number_of_triangles = 0; + + if (export_surface_only) + { + for (Cells_iterator cit = tr.cells_begin() ; + cit != tr.cells_end() ; ++cit) + { + if (tr.is_infinite(cit)) + { + output << "3 "; + for (int i = 0 ; i < 4 ; ++i) + { + if (!tr.is_infinite(cit->vertex(i))) + output << index_of_vertex[cit->vertex(i)] << " "; + } + output << std::endl; + ++number_of_triangles; + } + } + } + else + { + for (Finite_cells_iterator cit = tr.finite_cells_begin() ; + cit != tr.finite_cells_end() ; ++cit) + { + output << "3 " + << index_of_vertex[cit->vertex(0)] << " " + << index_of_vertex[cit->vertex(1)] << " " + << index_of_vertex[cit->vertex(2)] + << std::endl; + output << "3 " + << index_of_vertex[cit->vertex(0)] << " " + << index_of_vertex[cit->vertex(2)] << " " + << index_of_vertex[cit->vertex(3)] + << std::endl; + output << "3 " + << index_of_vertex[cit->vertex(1)] << " " + << index_of_vertex[cit->vertex(2)] << " " + << index_of_vertex[cit->vertex(3)] + << std::endl; + output << "3 " + << index_of_vertex[cit->vertex(0)] << " " + << index_of_vertex[cit->vertex(1)] << " " + << index_of_vertex[cit->vertex(3)] + << std::endl; + number_of_triangles += 4; + } + } + + os << "OFF \n" + << n << " " + << number_of_triangles << " 0\n" + << output.str(); + + return os; +} + +} //namespace CGAL + +#endif // CGAL_TRIANGULATION_OFF_OSTREAM_3_H \ No newline at end of file