diff --git a/Mesh_3/demo/Mesh_3/Mesh_function.h b/Mesh_3/demo/Mesh_3/Mesh_function.h index b9e8e8e7147..d6d6a5a7f01 100644 --- a/Mesh_3/demo/Mesh_3/Mesh_function.h +++ b/Mesh_3/demo/Mesh_3/Mesh_function.h @@ -204,7 +204,7 @@ launch() /*mesher_->initialize(); #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif while ( ! mesher_->is_algorithm_done() && continue_ ) diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 3b3f57c15c1..1a0f6e0e697 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -2607,7 +2607,7 @@ rebuild_restricted_delaunay(OutdatedCells& outdated_cells, # ifdef CGAL_MESH_3_PROFILING std::cerr << std::endl << " Updating cells..."; - WallClockTimer t; + Wall_clock_timer t; size_t num_cells = c3t3_.number_of_cells_in_complex(); # endif diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h index 80f52a349ab..f5fcea26540 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h @@ -311,7 +311,7 @@ private: #ifdef CGAL_MESH_3_PROFILING std::cerr << "Computing moves..."; - WallClockTimer t; + Wall_clock_timer t; #endif @@ -690,7 +690,7 @@ operator()(int nb_iterations, Visitor visitor) running_time_.start(); #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif // Fill set containing moving vertices @@ -898,7 +898,7 @@ update_mesh(const Moves_vector& moves, #ifdef CGAL_MESH_3_PROFILING std::cerr << "Moving vertices..."; - WallClockTimer t; + Wall_clock_timer t; #endif #ifdef CGAL_LINKED_WITH_TBB diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h index cea06d3fbf6..56e6ad56eec 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h @@ -320,7 +320,7 @@ Mesher_3::refine_mesh(std::string dump_after_refine_surface_prefix) #ifdef CGAL_MESH_3_PROFILING std::cerr << "Refining facets..." << std::endl; - WallClockTimer t; + Wall_clock_timer t; #endif facets_mesher_.refine(facets_visitor_); #ifdef CGAL_MESH_3_PROFILING @@ -476,7 +476,7 @@ initialize() { #ifdef CGAL_MESH_3_PROFILING std::cerr << "Initializing... "; - WallClockTimer t; + Wall_clock_timer t; #endif //===================================== // Bounding box estimation diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h index 509be7abe19..79943874df4 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h @@ -194,7 +194,7 @@ protected: #ifdef CGAL_MESH_3_PROFILING protected: - WallClockTimer m_timer; + Wall_clock_timer m_timer; #endif public: diff --git a/Mesh_3/include/CGAL/Mesh_3/Profiling_tools.h b/Mesh_3/include/CGAL/Mesh_3/Profiling_tools.h index eb436c93af3..32c707f9d9e 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Profiling_tools.h +++ b/Mesh_3/include/CGAL/Mesh_3/Profiling_tools.h @@ -28,10 +28,10 @@ // TBB timers #ifdef CGAL_LINKED_WITH_TBB #include - struct WallClockTimer + struct Wall_clock_timer { tbb::tick_count t; - WallClockTimer() + Wall_clock_timer() { t = tbb::tick_count::now(); } @@ -48,10 +48,10 @@ #else #include - struct WallClockTimer + struct Wall_clock_timer { CGAL::Real_timer t; - WallClockTimer() + Wall_clock_timer() { t.start(); } diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h index 776cd4f4796..8b6e17cb96b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h @@ -632,7 +632,7 @@ scan_triangulation_impl() typedef typename Tr::Finite_cells_iterator Finite_cell_iterator; #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h index 1bdcc3b66d9..a58627b667a 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h @@ -882,7 +882,7 @@ scan_triangulation_impl() typedef typename Tr::Finite_facets_iterator Finite_facet_iterator; #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif #ifdef CGAL_MESH_3_VERY_VERBOSE diff --git a/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h b/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h index fd3847f41ec..74fb202daf1 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h +++ b/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h @@ -799,7 +799,7 @@ operator()(Visitor visitor) running_time_.start(); #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif // Build priority queue (we use one queue for all steps) diff --git a/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h b/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h index bf5cf89ec60..5d8496fbf5c 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h +++ b/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h @@ -453,7 +453,7 @@ public: // methods operator()(Visitor visitor = Visitor()) { #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif Mesh_optimization_return_code ret = @@ -964,7 +964,7 @@ pump_vertices(double sliver_criterion_limit, Visitor& visitor) { #ifdef CGAL_MESH_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif init(sliver_criterion_limit); diff --git a/NewKernel_d/include/CGAL/Epeck_d.h b/NewKernel_d/include/CGAL/Epeck_d.h new file mode 100644 index 00000000000..52bce84ca3f --- /dev/null +++ b/NewKernel_d/include/CGAL/Epeck_d.h @@ -0,0 +1,53 @@ +// Copyright (c) 2014 +// INRIA Saclay-Ile de France (France) +// +// 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) : Marc Glisse + +#ifndef CGAL_EPECK_D_H +#define CGAL_EPECK_D_H +#include +#include +#include +#include + + +namespace CGAL { +#define CGAL_BASE \ + Cartesian_base_d::Type, Dim> +template +struct Epeck_d_help1 +: CGAL_BASE +{ + CGAL_CONSTEXPR Epeck_d_help1(){} + CGAL_CONSTEXPR Epeck_d_help1(int d):CGAL_BASE(d){} +}; +#undef CGAL_BASE +#define CGAL_BASE \ + Kernel_d_interface< \ + Cartesian_wrap< \ + Epeck_d_help1, \ + Epeck_d > > +template +struct Epeck_d +: CGAL_BASE +{ + CGAL_CONSTEXPR Epeck_d(){} + CGAL_CONSTEXPR Epeck_d(int d):CGAL_BASE(d){} +}; +#undef CGAL_BASE +} +#endif diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h index 403a0e4d68c..179e97bf016 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_filter_K.h @@ -30,7 +30,7 @@ namespace CGAL { template < typename Base_, typename AK_, typename EK_ > struct Cartesian_filter_K : public Base_, - private Store_kernel, private Store_kernel2 + private Store_kernel, private Store_kernel2 { CGAL_CONSTEXPR Cartesian_filter_K(){} CGAL_CONSTEXPR Cartesian_filter_K(int d):Base_(d){} @@ -61,11 +61,11 @@ struct Cartesian_filter_K : public Base_, template struct Type : Get_type {}; template::type> struct Functor : - Inherit_functor {}; + Inherit_functor {}; template struct Functor { - typedef typename Get_functor::type AP; - typedef typename Get_functor::type EP; - typedef Filtered_predicate2 type; + typedef typename Get_functor::type AP; + typedef typename Get_functor::type EP; + typedef Filtered_predicate2 type; }; // TODO: // template struct Functor : diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h b/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h index dcbea88278c..cd8269f07da 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h @@ -73,6 +73,7 @@ template struct Kernel_d_interface : public Base_ { typedef typename Get_functor::type In_flat_power_test_d; typedef typename Get_functor::type Point_to_vector_d; typedef typename Get_functor::type Vector_to_point_d; + typedef typename Get_functor::type Translated_point_d; typedef typename Get_functor::type Scaled_vector_d; typedef typename Get_functor::type Difference_of_vectors_d; typedef typename Get_functor::type Difference_of_points_d; @@ -184,6 +185,7 @@ template struct Kernel_d_interface : public Base_ { In_flat_power_test_d in_flat_power_test_d_object()const{ return In_flat_power_test_d(*this); } Point_to_vector_d point_to_vector_d_object()const{ return Point_to_vector_d(*this); } Vector_to_point_d vector_to_point_d_object()const{ return Vector_to_point_d(*this); } + Translated_point_d translated_point_d_object()const{ return Translated_point_d(*this); } Scaled_vector_d scaled_vector_d_object()const{ return Scaled_vector_d(*this); } Difference_of_vectors_d difference_of_vectors_d_object()const{ return Difference_of_vectors_d(*this); } Difference_of_points_d difference_of_points_d_object()const{ return Difference_of_points_d(*this); } diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h index 7ab10f16b99..aac6259cd09 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h @@ -57,7 +57,6 @@ template struct Construct_sphere : Store_kernel { typedef typename LA::Square_matrix Matrix; typedef typename LA::Vector Vec; typedef typename LA::Construct_vector CVec; - typedef typename Get_type::type Point; typename Get_functor::type c(this->kernel()); typename Get_functor >::type cp(this->kernel()); typename Get_functor::type pd(this->kernel()); diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h index 32e8b586b58..4c76c32d579 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h @@ -75,6 +75,33 @@ template struct Point_weight { } }; +template struct Power_distance : private Store_kernel { + CGAL_FUNCTOR_INIT_STORE(Power_distance) + typedef typename Get_type::type first_argument_type; + typedef first_argument_type second_argument_type; + typedef typename Get_type::type result_type; + + result_type operator()(first_argument_type const&a, second_argument_type const&b)const{ + typename Get_functor::type pdw(this->kernel()); + typename Get_functor::type pw(this->kernel()); + typename Get_functor::type sd(this->kernel()); + return sd(pdw(a),pdw(b))-pw(a)-pw(b); + } +}; +template struct Power_distance_to_point : private Store_kernel { + CGAL_FUNCTOR_INIT_STORE(Power_distance_to_point) + typedef typename Get_type::type first_argument_type; + typedef typename Get_type::type second_argument_type; + typedef typename Get_type::type result_type; + + result_type operator()(first_argument_type const&a, second_argument_type const&b)const{ + typename Get_functor::type pdw(this->kernel()); + typename Get_functor::type pw(this->kernel()); + typename Get_functor::type sd(this->kernel()); + return sd(pdw(a),b)-pw(a); + } +}; + template struct Power_test : private Store_kernel { CGAL_FUNCTOR_INIT_STORE(Power_test) typedef R_ R; @@ -114,6 +141,56 @@ template struct In_flat_power_test : private Store_kernel { } }; +// Construct a point at (weighted) distance 0 from all the input +template struct Power_center : Store_kernel { + CGAL_FUNCTOR_INIT_STORE(Power_center) + typedef typename Get_type::type WPoint; + typedef WPoint result_type; + typedef typename Get_type::type Point; + typedef typename Get_type::type FT; + template + result_type operator()(Iter f, Iter e)const{ + // 2*(x-y).c == (x^2-wx^2)-(y^2-wy^2) + typedef typename R_::LA LA; + typedef typename LA::Square_matrix Matrix; + typedef typename LA::Vector Vec; + typedef typename LA::Construct_vector CVec; + typename Get_functor::type c(this->kernel()); + typename Get_functor >::type cp(this->kernel()); + typename Get_functor::type pd(this->kernel()); + typename Get_functor::type sdo(this->kernel()); + typename Get_functor::type pdp(this->kernel()); + typename Get_functor::type pdw(this->kernel()); + typename Get_functor::type pw(this->kernel()); + typename Get_functor >::type cwp(this->kernel()); + + WPoint const& wp0 = *f; + Point const& p0 = pdw(wp0); + int d = pd(p0); + FT const& n0 = sdo(p0) - pw(wp0); + Matrix m(d,d); + Vec b = typename CVec::Dimension()(d); + // Write the point coordinates in lines. + int i; + for(i=0; ++f!=e; ++i) { + WPoint const& wp=*f; + Point const& p=pdw(wp); + FT const& np = sdo(p) - pw(wp); + for(int j=0;j),(Point_tag),()); CGAL_KD_DEFAULT_FUNCTOR(Construct_ttag,(CartesianDKernelFunctors::Construct_weighted_point),(Weighted_point_tag,Point_tag),()); @@ -121,5 +198,8 @@ CGAL_KD_DEFAULT_FUNCTOR(Point_drop_weight_tag,(CartesianDKernelFunctors::Point_d CGAL_KD_DEFAULT_FUNCTOR(Point_weight_tag,(CartesianDKernelFunctors::Point_weight),(Weighted_point_tag,Point_tag),()); CGAL_KD_DEFAULT_FUNCTOR(Power_test_tag,(CartesianDKernelFunctors::Power_test),(Weighted_point_tag),(Power_test_raw_tag,Point_drop_weight_tag,Point_weight_tag)); CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_test_tag,(CartesianDKernelFunctors::In_flat_power_test),(Weighted_point_tag),(In_flat_power_test_raw_tag,Point_drop_weight_tag,Point_weight_tag)); +CGAL_KD_DEFAULT_FUNCTOR(Power_distance_tag,(CartesianDKernelFunctors::Power_distance),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag)); +CGAL_KD_DEFAULT_FUNCTOR(Power_distance_to_point_tag,(CartesianDKernelFunctors::Power_distance_to_point),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag)); +CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center),(Weighted_point_tag,Point_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag,Construct_ttag,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag)); } // namespace CGAL #endif diff --git a/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h b/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h index afe0720955a..56036b09d0e 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h @@ -218,6 +218,8 @@ namespace CGAL { CGAL_DECL_COMPUTE(Hyperplane_translation); CGAL_DECL_COMPUTE(Value_at); CGAL_DECL_COMPUTE(Point_weight); + CGAL_DECL_COMPUTE(Power_distance); + CGAL_DECL_COMPUTE(Power_distance_to_point); #undef CGAL_DECL_COMPUTE #define CGAL_DECL_ITER_OBJ(X,Y,Z,C) struct X##_tag {}; \ @@ -268,6 +270,7 @@ namespace CGAL { CGAL_DECL_CONSTRUCT(Point_to_vector,Vector); CGAL_DECL_CONSTRUCT(Vector_to_point,Point); CGAL_DECL_CONSTRUCT(Point_drop_weight,Point); + CGAL_DECL_CONSTRUCT(Power_center,Weighted_point); #undef CGAL_DECL_CONSTRUCT #if 0 #define CGAL_DECL_ITER_CONSTRUCT(X,Y) struct X##_tag {}; \ diff --git a/NewKernel_d/test/NewKernel_d/Epick_d.cpp b/NewKernel_d/test/NewKernel_d/Epick_d.cpp index 0635bc4cc3f..96ca0b4c39e 100644 --- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp +++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp @@ -13,6 +13,7 @@ #include #include #include +#include //typedef CGAL::Cartesian_base_d > K0; //typedef CGAL::Cartesian_base_d > KA; @@ -73,6 +74,7 @@ void test2(){ typedef typename K1::Ray_d R; typedef typename K1::Iso_box_d IB; typedef typename K1::Flat_orientation_d FO; + typedef typename K1::Weighted_point_d WP; //typedef K1::Construct_point CP; typedef typename K1::Construct_point_d CP; @@ -120,6 +122,13 @@ void test2(){ typedef typename K1::Scalar_product_d SP; typedef typename K1::Difference_of_vectors_d DV; typedef typename K1::Difference_of_points_d DP; + typedef typename K1::Translated_point_d TP; + typedef typename CGAL::Get_functor::type PC; + typedef typename CGAL::Get_functor::type PoD; + typedef typename K1::Weighted_point_d WP; + typedef typename K1::Construct_weighted_point_d CWP; + typedef typename K1::Point_drop_weight_d PDW; + typedef typename K1::Point_weight_d PW; CGAL_USE_TYPE(AT); CGAL_USE_TYPE(D); @@ -177,6 +186,12 @@ void test2(){ SP spr Kinit(scalar_product_d_object); DV dv Kinit(difference_of_vectors_d_object); DP dp Kinit(difference_of_points_d_object); + TP tp Kinit(translated_point_d_object); + PC pc (k); + CWP cwp Kinit(construct_weighted_point_d_object); + PDW pdw Kinit(point_drop_weight_d_object); + PW pw Kinit(point_weight_d_object); + PoD pod (k); CGAL_USE(bc); CGAL_USE(pol); @@ -185,11 +200,12 @@ void test2(){ CGAL_USE(cli); CGAL_USE(cr); CGAL_USE(cib); + using std::abs; P a=cp(3,4); assert(pd(a)==2); assert(pv(a)[1]==4); P b=vp(cv(5,6,7)); - assert(fabs(b[0]-5./7)<.0001); + assert(abs(b[0]-5./7)<.0001); assert(lc(b,a,1)); assert(!lc(a,b,0)); int rr[]={3,5,2}; @@ -203,8 +219,8 @@ void test2(){ assert(cl(a,c)==CGAL::SMALLER); assert(ll(b,c)); assert(cl(c,b)==CGAL::LARGER); - assert(fabs(m(a,c)[0]-3)<.0001); - assert(fabs(m(a,c)[1]-4.5)<.0001); + assert(abs(m(a,c)[0]-3)<.0001); + assert(abs(m(a,c)[1]-4.5)<.0001); P d=cp(r,r+3,CGAL::Homogeneous_tag()); S s=cs(c,d); std::cout << cc(a,1) << std::endl; @@ -253,9 +269,9 @@ void test2(){ assert(v.size()==1); assert(lr(tab3+0,tab3+2)==1); H h=ch(tab2+1,tab2+3); - assert(fabs(va(h,x2)-1)<.0001); - assert(fabs(va(h,x3)-1)<.0001); - assert(fabs(va(h,x1)+1)<.0001); + assert(abs(va(h,x2)-1)<.0001); + assert(abs(va(h,x3)-1)<.0001); + assert(abs(va(h,x1)+1)<.0001); H h2=ch(tab2+1,tab2+3,x1,CGAL::ON_POSITIVE_SIDE); assert(hops(h2,x1)); assert(os(h2,x1)==CGAL::ON_POSITIVE_SIDE); @@ -312,20 +328,35 @@ void test2(){ Sp sp = csp(tabz+0,tabz+3); P cent0=cos(sp); P cent1=cos(tabz+0,tabz+3); - assert(fabs(cent0[0]-2)<.0001); - assert(fabs(cent0[1]+3)<.0001); - assert(fabs(cent1[0]-2)<.0001); - assert(fabs(cent1[1]+3)<.0001); - assert(fabs(sp.squared_radius()-25)<.0001); + assert(abs(cent0[0]-2)<.0001); + assert(abs(cent0[1]+3)<.0001); + assert(abs(cent1[0]-2)<.0001); + assert(abs(cent1[1]+3)<.0001); + assert(abs(sp.squared_radius()-25)<.0001); +#if 1 + // Fails for an exact kernel P psp0=ps(sp,0); P psp1=ps(sp,1); P psp2=ps(sp,2); assert(!ed(psp0,psp1)); assert(!ed(psp0,psp2)); assert(!ed(psp2,psp1)); - assert(fabs(sd(cent0,psp0)-25)<.0001); - assert(fabs(sd(cent0,psp1)-25)<.0001); - assert(fabs(sd(cent0,psp2)-25)<.0001); + assert(abs(sd(cent0,psp0)-25)<.0001); + assert(abs(sd(cent0,psp1)-25)<.0001); + assert(abs(sd(cent0,psp2)-25)<.0001); +#endif + P x2py1 = tp(x2,y1); + assert(x2py1[1]==-2); + WP tw[]={cwp(cp(5,0),1.5),cwp(cp(2,std::sqrt(3)),1),cwp(cp(2,-std::sqrt(3)),1)}; + WP xw=pc(tw+0,tw+3); + assert(abs(pod(xw,tw[0]))<.0001); + assert(abs(pod(xw,tw[1]))<.0001); + assert(abs(pod(xw,tw[2]))<.0001); + assert(pdw(xw)[0]<2.95); + assert(pdw(xw)[0]>2.5); + assert(pw(xw)<2.95); + assert(pw(xw)>2.5); + Sp un1; CGAL_USE(un1); H un2; CGAL_USE(un2); S un3; CGAL_USE(un3); @@ -418,6 +449,7 @@ void test3(){ SD sd Kinit(squared_distance_d_object); PD pd Kinit(point_dimension_d_object); AI ai Kinit(affinely_independent_d_object); + using std::abs; P a; // Triangulation needs this :-( a=cp(2,3,4); assert(pd(a)==3); @@ -441,7 +473,7 @@ void test3(){ std::cout << *i << ' '; std::cout << '\n'; P e=cp(-2,3,0); - assert(fabs(sd(e,a)-32)<.0001); + assert(abs(sd(e,a)-32)<.0001); P tab[]={a,b,c,d,e}; std::cout << po (&tab[0],tab+4) << ' '; std::cout << sos(&tab[0],tab+5) << ' '; diff --git a/Spatial_searching/include/CGAL/Search_traits_adapter.h b/Spatial_searching/include/CGAL/Search_traits_adapter.h index 1526c084885..c5d629a0358 100644 --- a/Spatial_searching/include/CGAL/Search_traits_adapter.h +++ b/Spatial_searching/include/CGAL/Search_traits_adapter.h @@ -89,16 +89,26 @@ public: struct Construct_cartesian_const_iterator_d: public Base_traits::Construct_cartesian_const_iterator_d{ PointPropertyMap ppmap; - using Base_traits::Construct_cartesian_const_iterator_d::operator(); - + Construct_cartesian_const_iterator_d(const typename Base_traits::Construct_cartesian_const_iterator_d& base, const PointPropertyMap& ppmap_) :Base_traits::Construct_cartesian_const_iterator_d(base), ppmap(ppmap_){} typename Base_traits::Cartesian_const_iterator_d operator()(const Point_with_info& p) const - { return this->operator() (get(ppmap,p)); } + { return Base_traits::Construct_cartesian_const_iterator_d::operator() (get(ppmap,p)); } typename Base_traits::Cartesian_const_iterator_d operator()(const Point_with_info& p, int) const - { return this->operator() (get(ppmap,p),0); } + { return Base_traits::Construct_cartesian_const_iterator_d::operator() (get(ppmap,p),0); } + + // These 2 additional operators forward the call to Base_traits. + // This is needed because of an undocumented requirement of + // Orthogonal_k_neighbor_search and Orthogonal_incremental_neighbor_search: + // Traits::Construct_cartesian_const_iterator should be callable + // on the query point type + typename Base_traits::Cartesian_const_iterator_d operator()(const typename Base_traits::Point_d& p) const + { return Base_traits::Construct_cartesian_const_iterator_d::operator() (p); } + + typename Base_traits::Cartesian_const_iterator_d operator()(const typename Base_traits::Point_d& p, int) const + { return Base_traits::Construct_cartesian_const_iterator_d::operator() (p,0); } }; struct Construct_iso_box_d: public Base::Construct_iso_box_d{ diff --git a/Tangential_complex/benchmark/Tangential_complex/CMakeLists.txt b/Tangential_complex/benchmark/Tangential_complex/CMakeLists.txt new file mode 100644 index 00000000000..7155f246c2b --- /dev/null +++ b/Tangential_complex/benchmark/Tangential_complex/CMakeLists.txt @@ -0,0 +1,77 @@ +# Created by the script cgal_create_cmake_script +# This is the CMake script for compiling a CGAL application. + + +project( Tangential_complex_benchmark ) + +cmake_minimum_required(VERSION 2.6.2) +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6) + if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3) + cmake_policy(VERSION 2.8.4) + else() + cmake_policy(VERSION 2.6) + endif() +endif() + +# Creates a new CMake option, turned ON by default +option(ACTIVATE_MSVC_PRECOMPILED_HEADERS + "Activate precompiled headers in MSVC" + OFF) + +# Macro to add precompiled headers for MSVC +# This function does two things: +# 1. Enable precompiled headers on each file which is listed in "SourcesVar". +# 2. Add the content of "PrecompiledSource" (e.g. "StdAfx.cpp") to "SourcesVar". +MACRO(ADD_MSVC_PRECOMPILED_HEADER PrecompiledHeader PrecompiledSource SourcesVar) + IF(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS) + GET_FILENAME_COMPONENT(PrecompiledBasename ${PrecompiledHeader} NAME_WE) + SET(Sources ${${SourcesVar}}) + + SET_SOURCE_FILES_PROPERTIES(${PrecompiledSource} + PROPERTIES COMPILE_FLAGS "/Yc\"${PrecompiledHeader}\"") + SET_SOURCE_FILES_PROPERTIES(${Sources} + PROPERTIES COMPILE_FLAGS "/Yu\"${PrecompiledHeaders}\" /FI\"${PrecompiledHeader}\"") + # Add precompiled header to SourcesVar + LIST(APPEND ${SourcesVar} ${PrecompiledSource}) + ENDIF(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS) +ENDMACRO(ADD_MSVC_PRECOMPILED_HEADER) + +# The compiler might need more memory because of precompiled headers +if(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS AND NOT(MSVC_VERSION LESS 1310)) + set(CGAL_C_FLAGS "${CGAL_C_FLAGS} /Zm1000") + set(CGAL_CXX_FLAGS "${CGAL_CXX_FLAGS} /Zm1000") +endif() + +find_package(CGAL QUIET COMPONENTS Core ) + +if ( CGAL_FOUND ) + include( ${CGAL_USE_FILE} ) + + find_package( TBB QUIET ) + + if( TBB_FOUND ) + include(${TBB_USE_FILE}) + list(APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES}) + endif() + + + include( CGAL_CreateSingleSourceCGALProgram ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + include_directories (BEFORE "include") + + set (SOURCE_FILES "benchmark_tc.cpp") + ADD_MSVC_PRECOMPILED_HEADER("StdAfx.h" "StdAfx.cpp" SOURCE_FILES) + create_single_source_cgal_program( ${SOURCE_FILES} ) + + 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/Tangential_complex/benchmark/Tangential_complex/StdAfx.cpp b/Tangential_complex/benchmark/Tangential_complex/StdAfx.cpp new file mode 100644 index 00000000000..15668dcadef --- /dev/null +++ b/Tangential_complex/benchmark/Tangential_complex/StdAfx.cpp @@ -0,0 +1,2 @@ +// Build the precompiled headers. +#include "StdAfx.h" \ No newline at end of file diff --git a/Tangential_complex/benchmark/Tangential_complex/StdAfx.h b/Tangential_complex/benchmark/Tangential_complex/StdAfx.h new file mode 100644 index 00000000000..573345d0022 --- /dev/null +++ b/Tangential_complex/benchmark/Tangential_complex/StdAfx.h @@ -0,0 +1,14 @@ +#ifndef STDAFX_H +#define STDAFX_H + +// CGAL +#include +#include +#include +#include +#include +#include +#include +#include + +#endif //STDAFX_H \ No newline at end of file diff --git a/Tangential_complex/benchmark/Tangential_complex/XML_exporter.h b/Tangential_complex/benchmark/Tangential_complex/XML_exporter.h new file mode 100644 index 00000000000..158d3d648d4 --- /dev/null +++ b/Tangential_complex/benchmark/Tangential_complex/XML_exporter.h @@ -0,0 +1,236 @@ +// Copyright (c) 2012 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 +// + +#include +#include +#include +#include +#include + +template +class Simple_XML_exporter +{ +public: + typedef value_type Value_type; + typedef std::vector Element; + typedef std::map Element_with_map; + typedef std::vector List_of_elements; + + Simple_XML_exporter( + const std::string &list_name, + const std::string &element_name, + const std::vector &subelement_names, + bool add_timestamp = true) + : m_list_name(list_name), + m_element_name(element_name), + m_subelement_names(subelement_names), + m_add_timestamp(add_timestamp) + {} + + bool add_element(const Element &element) + { + if (element.size() == m_subelement_names.size()) + { + m_list_of_elements.push_back(element); + return true; + } + else + { + std::cerr << "ERROR: element.size() == m_subelement_names.size()" << std::endl; + return false; + } + } + + bool add_element(Element_with_map &element) + { + Element elt; + + std::vector::const_iterator + it_subelement_name = m_subelement_names.begin(); + std::vector::const_iterator + it_subelement_name_end = m_subelement_names.end(); + for ( ; it_subelement_name != it_subelement_name_end ; ++it_subelement_name) + { + elt.push_back(element[*it_subelement_name]); + } + + return add_element(elt); + } + + bool export_to_xml(const std::string &filename) const + { + std::ofstream xmlfile; + xmlfile.open(filename.c_str()); + xmlfile << "" << std::endl; + xmlfile << "<" << m_list_name << ">" << std::endl; + + typename List_of_elements::const_iterator it_element = m_list_of_elements.begin(); + typename List_of_elements::const_iterator it_element_end = m_list_of_elements.end(); + for (int id = 1 ; it_element != it_element_end ; ++it_element, ++id) + { + xmlfile << " <" << m_element_name << ">" << std::endl; + std::vector::const_iterator + it_subelement_name = m_subelement_names.begin(); + std::vector::const_iterator + it_subelement_name_end = m_subelement_names.end(); + + if (m_add_timestamp) + xmlfile << " " << time(NULL) << " " << std::endl; + + for (int i = 0 ; + it_subelement_name != it_subelement_name_end ; + ++it_subelement_name, ++i) + { + xmlfile + << " <" << *it_subelement_name << "> " + << (*it_element)[i] + << " " << std::endl; + } + xmlfile << " " << std::endl; + } + + xmlfile << "" << std::endl; + xmlfile.close(); + return 0; + + } + +protected: + std::string m_list_name; + std::string m_element_name; + std::vector m_subelement_names; + List_of_elements m_list_of_elements; + bool m_add_timestamp; +}; + + + + + +template +class Streaming_XML_exporter +{ +public: + typedef value_type Value_type; + typedef std::vector Element; + typedef std::map Element_with_map; + typedef std::vector List_of_elements; + + Streaming_XML_exporter( + const std::string &filename, + const std::string &list_name, + const std::string &element_name, + const std::vector &subelement_names, + bool add_timestamp = true) + : m_list_name(list_name), + m_element_name(element_name), + m_subelement_names(subelement_names), + m_add_timestamp(add_timestamp) + { + m_xml_fstream.open(filename.c_str()); + if (m_xml_fstream.good()) + { + m_xml_fstream << "" << std::endl; + m_xml_fstream << "<" << m_list_name << ">" << std::endl; + } + else + { + std::cerr << "Could not open file '" << filename << "'." << std::endl; + } + } + + virtual ~Streaming_XML_exporter() + { + close_file(); + } + + void close_file() + { + m_xml_fstream.close(); + } + + bool add_element(const Element &element) + { + if (element.size() == m_subelement_names.size()) + { + m_xml_fstream << " <" << m_element_name << ">" << std::endl; + std::vector::const_iterator + it_subelement_name = m_subelement_names.begin(); + std::vector::const_iterator + it_subelement_name_end = m_subelement_names.end(); + + if (m_add_timestamp) + { + m_xml_fstream << " " << time(NULL) << " " << std::endl; + } + + for (int i = 0 ; + it_subelement_name != it_subelement_name_end ; + ++it_subelement_name, ++i) + { + m_xml_fstream + << " <" << *it_subelement_name << "> " + << element[i] + << " " << std::endl; + } + m_xml_fstream << " " << std::endl; + + // Save current pointer position + std::ofstream::streampos pos = m_xml_fstream.tellp(); + // Close the XML file (temporarily) so that the XML file is always correct + m_xml_fstream << "" << std::endl; + // Restore the pointer position so that the next "add_element" will overwrite + // the end of the file + m_xml_fstream.seekp(pos); + + m_xml_fstream.flush(); + return true; + } + else + { + std::cerr << "ERROR: element.size() == m_subelement_names.size()" << std::endl; + return false; + } + } + + bool add_element(Element_with_map &element) + { + Element elt; + + std::vector::const_iterator + it_subelement_name = m_subelement_names.begin(); + std::vector::const_iterator + it_subelement_name_end = m_subelement_names.end(); + for ( ; it_subelement_name != it_subelement_name_end ; ++it_subelement_name) + { + elt.push_back(element[*it_subelement_name]); + } + + return add_element(elt); + } + +protected: + std::ofstream m_xml_fstream; + std::string m_list_name; + std::string m_element_name; + std::vector m_subelement_names; + bool m_add_timestamp; +}; diff --git a/Tangential_complex/benchmark/Tangential_complex/benchmark_tc.cpp b/Tangential_complex/benchmark/Tangential_complex/benchmark_tc.cpp new file mode 100644 index 00000000000..26606958ab2 --- /dev/null +++ b/Tangential_complex/benchmark/Tangential_complex/benchmark_tc.cpp @@ -0,0 +1,231 @@ +#include +#include +#include +#include + +#include "../../test/Tangential_complex/test_utilities.h" + + +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +# include +#endif +#include "XML_exporter.h" +#define CGAL_TC_EXPORT_PERFORMANCE_DATA +#define CGAL_TC_SET_PERFORMANCE_DATA(value_name, value) \ + XML_perf_data::set(value_name, value); + +class XML_perf_data +{ +public: + typedef Streaming_XML_exporter XML_exporter; + + XML_perf_data(const std::string &filename) + : m_xml(filename, "ContainerPerformance", "Perf", + construct_subelements_names()) + {} + + virtual ~XML_perf_data() + { + } + + static XML_perf_data &get() + { + static XML_perf_data singleton(build_filename()); + return singleton; + } + + template + static void set(const std::string &name, Value_type value) + { + get().set_data(name, value); + } + + static void commit() + { + get().commit_current_element(); + } + +protected: + static std::string build_filename() + { + std::stringstream sstr; + sstr << "perf_logs/Performance_log_" << time(0) << ".xml"; + return sstr.str(); + } + + static std::vector construct_subelements_names() + { + std::vector subelements; + subelements.push_back("Name"); + subelements.push_back("Intrinsic_dim"); + subelements.push_back("Ambient_dim"); + subelements.push_back("Num_points"); + subelements.push_back("Noise"); + subelements.push_back("Info"); + subelements.push_back("Init_time"); + subelements.push_back("Comput_time"); + subelements.push_back("Fix_time"); + subelements.push_back("Fix_steps"); + + return subelements; + } + + void set_data(const std::string &name, const std::string &value) + { + m_current_element[name] = value; + } + + template + void set_data(const std::string &name, Value_type value) + { + std::stringstream sstr; + sstr << value; + set_data(name, sstr.str()); + } + + void commit_current_element() + { + m_xml.add_element(m_current_element); + m_current_element.clear(); + } + + XML_exporter m_xml; + XML_exporter::Element_with_map m_current_element; +}; + +const double NOISE = 0.01; +#ifdef _DEBUG + const int NUM_POINTS = 50; +#else + const int NUM_POINTS = 30000; +#endif + +int main() +{ +#ifdef CGAL_LINKED_WITH_TBB +# ifdef _DEBUG + tbb::task_scheduler_init init(1); +# else + tbb::task_scheduler_init init(10); +# endif +#endif + + const int INTRINSIC_DIMENSION = 2; + const int AMBIENT_DIMENSION = 4; + + typedef CGAL::Epick_d > Kernel; + typedef Kernel::FT FT; + typedef Kernel::Point_d Point; + + bool stop = false; + //for (int i = 0, NUM_POINTS = 10000 ; + // NUM_POINTS <= 50000 && !stop ; + // ++i, NUM_POINTS += (i%3 == 0 ? 5000 : 0)) + for (int i = 0 ; i < 5 && !stop ; ++i) + { + Kernel k; + Wall_clock_timer t; + CGAL::default_random = CGAL::Random(i); + std::cerr << "Random seed = " << i << std::endl; + +#ifdef CGAL_TC_PROFILING + Wall_clock_timer t_gen; +#endif + //std::vector points; + //load_points_from_file( + // "data/points_10_10k.cin", std::back_inserter(points)/*, 600*/); + + std::vector points = + //generate_points_on_circle_2(NUM_POINTS, 3.); + //generate_points_on_moment_curve(NUM_POINTS, AMBIENT_DIMENSION, 0., 1.); + //generate_points_on_plane(NUM_POINTS); + //generate_points_on_sphere_3(NUM_POINTS, 3.0); + //generate_points_on_sphere_d(NUM_POINTS, AMBIENT_DIMENSION, 3.0); + //generate_points_on_klein_bottle_3D(NUM_POINTS, 4., 3.); + generate_points_on_klein_bottle_4D(NUM_POINTS, 4., 3., NOISE); + //generate_points_on_klein_bottle_variant_5D(NUM_POINTS, 4., 3.); + + + CGAL_TC_SET_PERFORMANCE_DATA("Name", "Klein bottle 4D"); + CGAL_TC_SET_PERFORMANCE_DATA("Intrinsic_dim", INTRINSIC_DIMENSION); + CGAL_TC_SET_PERFORMANCE_DATA("Ambient_dim", AMBIENT_DIMENSION); + CGAL_TC_SET_PERFORMANCE_DATA("Noise", NOISE); + XML_perf_data::set("Info", "" +#ifdef CGAL_TC_ONLY_CHANGE_SIMPLEX_WEIGHTS + "ONLY_CHANGE_SIMPLEX_WEIGHTS " +#endif +#ifdef CGAL_TC_GLOBAL_REFRESH + "GLOBAL_REFRESH " +#endif + ); + +#ifdef CGAL_TC_PROFILING + std::cerr << "Point set generated in " << t_gen.elapsed() + << " seconds." << std::endl; +#endif + + std::cerr << "Number of points before sparsification: " + << points.size() << std::endl; + points = sparsify_point_set( + k, points, FT(INPUT_SPARSITY)*FT(INPUT_SPARSITY)); + std::cerr << "Number of points after sparsification: " + << points.size() << std::endl; + + CGAL::Tangential_complex< + Kernel, + INTRINSIC_DIMENSION, + CGAL::Parallel_tag> tc(points.begin(), points.end(), k); + double init_time = t.elapsed(); t.reset(); + + //tc.estimate_intrinsic_dimension(); + tc.compute_tangential_complex(); + double computation_time = t.elapsed(); t.reset(); + + std::set > incorrect_simplices; + //stop = !tc.check_if_all_simplices_are_in_the_ambient_delaunay(&incorrect_simplices); + + t.reset(); + unsigned int num_fix_steps = tc.fix_inconsistencies(); + double fix_time = t.elapsed(); t.reset(); + + double export_time = -1.; + if (INTRINSIC_DIMENSION <= 3) + { + t.reset(); + std::stringstream output_filename; + output_filename << "data/test_tc_" << INTRINSIC_DIMENSION + << "_in_R" << AMBIENT_DIMENSION << ".off"; + std::ofstream off_stream(output_filename.str().c_str()); + tc.export_to_off(off_stream, true, &incorrect_simplices, true); + double export_time = t.elapsed(); t.reset(); + } + /*else + tc.number_of_inconsistent_simplices();*/ + + + std::cerr << std::endl + << "================================================" << std::endl + << "Number of vertices: " << tc.number_of_vertices() << std::endl + << "Computation times (seconds): " << std::endl + << " * Tangential complex: " << init_time + computation_time < +const double SQ_HALF_SPARSITY = 0.5*0.5*INPUT_SPARSITY*INPUT_SPARSITY; + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef CGAL_TC_PROFILING +# include +#endif + +#include // CJTODO TEMP + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +# include +//# include +#endif + +//#define CGAL_TC_EXPORT_NORMALS // Only for 3D surfaces (k=2, d=3) + +namespace CGAL { + +using namespace Tangential_complex_; + +class Vertex_data +{ +public: + Vertex_data(std::size_t data = std::numeric_limits::max()) + : m_data(data) + {} + operator std::size_t() { return m_data; } + operator std::size_t() const { return m_data; } + +private: + std::size_t m_data; +}; + +/// The class Tangential_complex represents a tangential complex +template < + typename Kernel, + int Intrinsic_dimension, + typename Concurrency_tag = CGAL::Parallel_tag, + typename Tr = Regular_triangulation + < + Regular_triangulation_euclidean_traits< + Epick_d > >, + + Triangulation_data_structure + < + typename Regular_triangulation_euclidean_traits< + Epick_d > >::Dimension, + Triangulation_vertex > >, Vertex_data >, + Triangulation_full_cell > > > + > + > +> +class Tangential_complex +{ + typedef typename Kernel::FT FT; + typedef typename Kernel::Point_d Point; + typedef typename Kernel::Vector_d Vector; + + typedef Tr Triangulation; + typedef typename Triangulation::Geom_traits Tr_traits; + typedef typename Triangulation::Weighted_point Tr_point; + typedef typename Triangulation::Bare_point Tr_bare_point; + typedef typename Triangulation::Vertex_handle Tr_vertex_handle; + typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle; + + typedef std::vector Tangent_space_basis; + + typedef std::vector Points; + typedef std::vector Weights; + typedef Point_cloud_data_structure Points_ds; + typedef typename Points_ds::KNS_range KNS_range; + typedef typename Points_ds::KNS_iterator KNS_iterator; + typedef typename Points_ds::INS_range INS_range; + typedef typename Points_ds::INS_iterator INS_iterator; + + // Store a local triangulation and a handle to its center vertex + struct Tr_and_VH + { + public: + Tr_and_VH() + : m_tr(NULL) {} + Tr_and_VH(int dim) + : m_tr(new Triangulation(dim)) {} + + ~Tr_and_VH() { destroy_triangulation(); } + + Triangulation & construct_triangulation(int dim) + { + delete m_tr; + m_tr = new Triangulation(dim); + return tr(); + } + + void destroy_triangulation() + { + delete m_tr; + m_tr = NULL; + } + + Triangulation & tr() { return *m_tr; } + Triangulation const& tr() const { return *m_tr; } + + + Tr_vertex_handle const& center_vertex() const { return m_center_vertex; } + Tr_vertex_handle & center_vertex() { return m_center_vertex; } + + private: + Triangulation* m_tr; + Tr_vertex_handle m_center_vertex; + }; + + typedef typename std::vector TS_container; + typedef typename std::vector Tr_container; +#ifdef CGAL_LINKED_WITH_TBB + // CJTODO: test other mutexes + // http://www.threadingbuildingblocks.org/docs/help/reference/synchronization/mutexes/mutex_concept.htm + //typedef tbb::queuing_mutex Tr_mutex; +#endif +#ifdef CGAL_TC_EXPORT_NORMALS + typedef typename std::vector Normals; +#endif + +public: + + /// Constructor for a range of points + template + Tangential_complex(InputIterator first, InputIterator last, + const Kernel &k = Kernel()) + : m_k(k), + m_points(first, last), + m_points_ds(m_points) + {} + + /// Destructor + ~Tangential_complex() {} + + std::size_t number_of_vertices() + { + return m_points.size(); + } + + void compute_tangential_complex() + { +#ifdef CGAL_TC_PROFILING + Wall_clock_timer t; +#endif + + // We need to do that because we don't want the container to copy the + // already-computed triangulations (while resizing) since it would + // invalidate the vertex handles stored beside the triangulations + m_triangulations.resize(m_points.size()); +#ifdef CGAL_LINKED_WITH_TBB + //m_tr_mutexes.resize(m_points.size()); +#endif + m_tangent_spaces.resize(m_points.size()); + m_weights.resize(m_points.size(), FT(0)); +#ifdef CGAL_TC_EXPORT_NORMALS + m_normals.resize(m_points.size()); +#endif + +#ifdef CGAL_LINKED_WITH_TBB + // Parallel + if (boost::is_convertible::value) + { + tbb::parallel_for(tbb::blocked_range(0, m_points.size()), + Compute_tangent_triangulation(*this) + ); + } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB + { + for (std::size_t i = 0 ; i < m_points.size() ; ++i) + compute_tangent_triangulation(i); + } + +#ifdef CGAL_TC_PROFILING + std::cerr << "Tangential complex computed in " << t.elapsed() + << " seconds." << std::endl; +#endif + } + + void estimate_intrinsic_dimension() + { + const int amb_dim = m_k.point_dimension_d_object()(*m_points.begin()); + + // Kernel functors + typename Kernel::Compute_coordinate_d coord = + m_k.compute_coordinate_d_object(); + + std::vector sum_eigen_values(amb_dim, FT(0)); + + Points::const_iterator it_p = m_points.begin(); + Points::const_iterator it_p_end = m_points.end(); + // For each point p + for ( ; it_p != it_p_end ; ++it_p) + { + const Point &p = *it_p; + + KNS_range kns_range = m_points_ds.query_ANN( + p, NUM_POINTS_FOR_PCA, false); + + //******************************* PCA ************************************* + + // One row = one point + Eigen::MatrixXd mat_points(NUM_POINTS_FOR_PCA, amb_dim); + KNS_iterator nn_it = kns_range.begin(); + for (int j = 0 ; + j < NUM_POINTS_FOR_PCA && nn_it != kns_range.end() ; + ++j, ++nn_it) + { + for (int i = 0 ; i < amb_dim ; ++i) + mat_points(j, i) = CGAL::to_double(coord(m_points[nn_it->first], i)); + } + Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean(); + Eigen::MatrixXd cov = centered.adjoint() * centered; + Eigen::SelfAdjointEigenSolver eig(cov); + + // The eigenvectors are sorted in increasing order of their corresponding + // eigenvalues + Tangent_space_basis ts; + for (int i = 0 ; i < amb_dim ; ++i) + sum_eigen_values[i] += eig.eigenvalues()[i]; + + //************************************************************************* + } + + // CJTODO: replace this by an actual estimation + for (FT v : sum_eigen_values) // CJTODO C++11 + { + std::cout << v << " "; + } + std::cout << "\n"; + } + + + void refresh_tangential_complex() + { +#ifdef CGAL_TC_PROFILING + Wall_clock_timer t; +#endif + +#ifdef CGAL_LINKED_WITH_TBB + // Parallel + if (boost::is_convertible::value) + { + tbb::parallel_for(tbb::blocked_range(0, m_points.size()), + Compute_tangent_triangulation(*this, + true) //tangent_spaces_are_already_computed + ); + } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB + { + for (std::size_t i = 0 ; i < m_points.size() ; ++i) + { + compute_tangent_triangulation(i, + true); // tangent_spaces_are_already_computed + } + } + + +#ifdef CGAL_TC_PROFILING + std::cerr << "Tangential complex refreshed in " << t.elapsed() + << " seconds." << std::endl; +#endif + } + + unsigned int fix_inconsistencies() + { + typename Kernel::Point_drop_weight_d drop_w = m_k.point_drop_weight_d_object(); + typename Kernel::Construct_weighted_point_d cwp = + m_k.construct_weighted_point_d_object(); + +#ifdef CGAL_TC_VERBOSE + std::cerr << "Fixing inconsistencies..." << std::endl; +#endif + + std::pair stats_before = + number_of_inconsistent_simplices(false); + +#ifdef CGAL_TC_VERBOSE + std::cerr << "Initial number of inconsistencies: " + << stats_before.second << std::endl; +#endif + + if (stats_before.second == 0) + { +#ifdef CGAL_TC_VERBOSE + std::cerr << "Nothing to fix." << std::endl; +#endif + return 0; + } + + bool done = false; + unsigned int num_steps = 0; + while (!done) + { +// CJTODO: the parallel version is not working for now +/*#ifdef CGAL_LINKED_WITH_TBB + // Parallel + if (boost::is_convertible::value) + { + tbb::parallel_for( + tbb::blocked_range(0, m_triangulations.size()), + Try_to_solve_inconsistencies_in_a_local_triangulation(*this) + ); + } + // Sequential + else +#endif // CGAL_LINKED_WITH_TBB*/ + { + for (std::size_t i = 0 ; i < m_triangulations.size() ; ++i) + try_to_solve_inconsistencies_in_a_local_triangulation(i); + } + +#ifdef CGAL_TC_GLOBAL_REFRESH + refresh_tangential_complex(); +#endif + + std::pair stats_after = + number_of_inconsistent_simplices(false); + +#ifdef CGAL_TC_VERBOSE + std::cerr << std::endl + << "================================================" << std::endl + << "Inconsistencies:\n" + << " * Number of vertices: " << m_points.size() << std::endl + << std::endl + << " * BEFORE fix_inconsistencies:" << std::endl + << " - Total number of simplices in stars (incl. duplicates): " + << stats_before.first << std::endl + << " - Number of inconsistent simplices in stars (incl. duplicates): " + << stats_before.second << std::endl + << " - Percentage of inconsistencies: " + << 100. * stats_before.second / stats_before.first << "%" + << std::endl + << std::endl + << " * AFTER fix_inconsistencies:" << std::endl + << " - Total number of simplices in stars (incl. duplicates): " + << stats_after.first << std::endl + << " - Number of inconsistent simplices in stars (incl. duplicates): " + << stats_after.second << std::endl + << " - Percentage of inconsistencies: " + << 100. * stats_after.second / stats_before.first << "%" + << std::endl + << "================================================" << std::endl; +#endif + done = (stats_after.second == 0); + stats_before = stats_after; + ++num_steps; + } + + return num_steps; + } + + + // Return a pair + std::pair number_of_inconsistent_simplices( +#ifdef CGAL_TC_VERBOSE + bool verbose = true +#else + bool verbose = false +#endif + ) + { + std::size_t num_simplices = 0; + std::size_t num_inconsistent_simplices = 0; + typename Tr_container::const_iterator it_tr = m_triangulations.begin(); + typename Tr_container::const_iterator it_tr_end = m_triangulations.end(); + // For each triangulation + for ( ; it_tr != it_tr_end ; ++it_tr) + { + Triangulation const& tr = it_tr->tr(); + Tr_vertex_handle center_vh = it_tr->center_vertex(); + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + typename std::vector::const_iterator it_c = + incident_cells.begin(); + typename std::vector::const_iterator it_c_end = + incident_cells.end(); + // For each cell + for ( ; it_c != it_c_end ; ++it_c) + { + if (tr.is_infinite(*it_c)) // Don't check infinite cells + continue; + + if (!is_simplex_consistent(*it_c)) + ++num_inconsistent_simplices; + ++num_simplices; + } + } + + if (verbose) + { + std::cerr << std::endl + << "================================================" << std::endl + << "Inconsistencies:\n" + << " * Number of vertices: " << m_points.size() << std::endl + << " * Total number of simplices in stars (incl. duplicates): " + << num_simplices << std::endl + << " * Number of inconsistent simplices in stars (incl. duplicates): " + << num_inconsistent_simplices << std::endl + << " * Percentage of inconsistencies: " + << 100 * num_inconsistent_simplices / num_simplices << "%" << std::endl + << "================================================" << std::endl; + } + + return std::make_pair(num_simplices, num_inconsistent_simplices); + } + + std::ostream &export_to_off( + std::ostream & os, + bool color_inconsistencies = false, + std::set > const* excluded_simplices = NULL, + bool show_excluded_vertices_in_color = false) + { + if (m_points.empty()) + return os; + + const int ambient_dim = m_k.point_dimension_d_object()(*m_points.begin()); + if (ambient_dim < 2) + { + std::cerr << "Error: export_to_off => ambient dimension should be >= 2." + << std::endl; + os << "Error: export_to_off => ambient dimension should be >= 2." + << std::endl; + return os; + } + if (ambient_dim > 3) + { + std::cerr << "Warning: export_to_off => ambient dimension should be " + "<= 3. Only the first 3 coordinates will be exported." + << std::endl; + } + + if (Intrinsic_dimension < 1 || Intrinsic_dimension > 3) + { + std::cerr << "Error: export_to_off => intrinsic dimension should be " + "between 1 and 3." + << std::endl; + os << "Error: export_to_off => intrinsic dimension should be " + "between 1 and 3." + << std::endl; + return os; + } + + std::stringstream output; + std::size_t num_simplices, num_vertices; + export_vertices_to_off(output, num_vertices); + export_simplices_to_off( + output, num_simplices, color_inconsistencies, + excluded_simplices, show_excluded_vertices_in_color); + +#ifdef CGAL_TC_EXPORT_NORMALS + os << "N"; +#endif + + os << "OFF \n" + << num_vertices << " " + << num_simplices << " " + << "0 \n" + << output.str(); + + return os; + } + + bool check_if_all_simplices_are_in_the_ambient_delaunay( + std::set > * incorrect_simplices = NULL) + { + if (m_points.empty()) + return true; + + const int ambient_dim = m_k.point_dimension_d_object()(*m_points.begin()); + typedef Delaunay_triangulation + > > DT; + typedef typename DT::Vertex_handle DT_VH; + typedef typename DT::Finite_full_cell_const_iterator FFC_it; + typedef std::set Indexed_simplex; + + //------------------------------------------------------------------------- + // Build the ambient Delaunay triangulation + // Then save its simplices into "amb_dt_simplices" + //------------------------------------------------------------------------- + + DT ambient_dt(ambient_dim); + for (std::size_t i=0; idata() = i; + } + + std::set amb_dt_simplices; + + for (FFC_it cit = ambient_dt.finite_full_cells_begin() ; + cit != ambient_dt.finite_full_cells_end() ; ++cit ) + { + CGAL::Combination_enumerator combi( + Intrinsic_dimension + 1, 0, ambient_dim + 1); + + for ( ; !combi.finished() ; ++combi) + { + Indexed_simplex simplex; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + simplex.insert(cit.base()->vertex(combi[i])->data()); + + amb_dt_simplices.insert(simplex); + } + } + + //------------------------------------------------------------------------- + // Parse the TC and save its simplices into "stars_simplices" + //------------------------------------------------------------------------- + + std::set stars_simplices; + + typename Tr_container::const_iterator it_tr = m_triangulations.begin(); + typename Tr_container::const_iterator it_tr_end = m_triangulations.end(); + // For each triangulation + for ( ; it_tr != it_tr_end ; ++it_tr) + { + Triangulation const& tr = it_tr->tr(); + Tr_vertex_handle center_vh = it_tr->center_vertex(); + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + typename std::vector::const_iterator it_c = + incident_cells.begin(); + typename std::vector::const_iterator it_c_end = + incident_cells.end(); + // For each cell + for ( ; it_c != it_c_end ; ++it_c) + { + if (tr.is_infinite(*it_c)) + { + std::cerr << "Warning: infinite cell in star" << std::endl; + continue; + } + Indexed_simplex simplex; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + simplex.insert((*it_c)->vertex(i)->data()); + + stars_simplices.insert(simplex); + } + } + + //------------------------------------------------------------------------- + // Check if simplices of "stars_simplices" are all in "amb_dt_simplices" + //------------------------------------------------------------------------- + + std::set diff; + if (!incorrect_simplices) + incorrect_simplices = &diff; + set_difference(stars_simplices.begin(), stars_simplices.end(), + amb_dt_simplices.begin(), amb_dt_simplices.end(), + std::inserter(*incorrect_simplices, + incorrect_simplices->begin()) ); + + if (!incorrect_simplices->empty()) + { + std::cerr + << "ERROR check_if_all_simplices_are_in_the_ambient_delaunay:" + << std::endl + << " Number of simplices in ambient DT: " << amb_dt_simplices.size() + << std::endl + << " Number of unique simplices in TC stars: " << stars_simplices.size() + << std::endl + << " Number of wrong simplices: " << incorrect_simplices->size() + << std::endl; + return false; + } + else + return true; + } + +private: + + class Compare_distance_to_ref_point + { + public: + Compare_distance_to_ref_point(Point const& ref, Kernel const& k) + : m_ref(ref), m_k(k) {} + + bool operator()(Point const& p1, Point const& p2) + { + typename Kernel::Squared_distance_d sqdist = + m_k.squared_distance_d_object(); + return sqdist(p1, m_ref) < sqdist(p2, m_ref); + } + + private: + Point const& m_ref; + Kernel const& m_k; + }; + + struct Tr_vertex_to_global_point + { + typedef Tr_vertex_handle argument_type; + typedef Point result_type; + + Tr_vertex_to_global_point(Points const& points) + : m_points(points) {} + + result_type operator()(argument_type const& vh) const + { + return m_points[vh->data()]; + } + + private: + Points const& m_points; + }; + + struct Tr_vertex_to_bare_point + { + typedef Tr_vertex_handle argument_type; + typedef Tr_bare_point result_type; + + Tr_vertex_to_bare_point(Tr_traits const& traits) + : m_traits(traits) {} + + result_type operator()(argument_type const& vh) const + { + typename Tr_traits::Point_drop_weight_d pdw = + m_traits.point_drop_weight_d_object(); + return pdw(vh->point()); + } + + private: + Tr_traits const& m_traits; + }; + +#ifdef CGAL_LINKED_WITH_TBB + // Functor for compute_tangential_complex function + class Compute_tangent_triangulation + { + Tangential_complex & m_tc; + bool m_tangent_spaces_are_already_computed; + + public: + // Constructor + Compute_tangent_triangulation( + Tangential_complex &tc, bool tangent_spaces_are_already_computed = false) + : m_tc(tc), + m_tangent_spaces_are_already_computed(tangent_spaces_are_already_computed) + {} + + // Constructor + Compute_tangent_triangulation(const Compute_tangent_triangulation &ctt) + : m_tc(ctt.m_tc) + {} + + // operator() + void operator()( const tbb::blocked_range& r ) const + { + for( size_t i = r.begin() ; i != r.end() ; ++i) + { + m_tc.compute_tangent_triangulation( + i, m_tangent_spaces_are_already_computed); + } + } + }; +#endif // CGAL_LINKED_WITH_TBB + + void compute_tangent_triangulation( + std::size_t i, bool tangent_spaces_are_already_computed = false) + { + //std::cerr << "***********************************************" << std::endl; + Triangulation &local_tr = + m_triangulations[i].construct_triangulation(Intrinsic_dimension); + const Tr_traits &local_tr_traits = local_tr.geom_traits(); + Tr_vertex_handle ¢er_vertex = m_triangulations[i].center_vertex(); + + // Kernel functor & objects + typename Kernel::Difference_of_points_d k_diff_pts = + m_k.difference_of_points_d_object(); + typename Kernel::Squared_distance_d k_sqdist = + m_k.squared_distance_d_object(); + + // Triangulation's traits functor & objects + typename Tr_traits::Squared_distance_d sqdist = + local_tr_traits.squared_distance_d_object(); + typename Tr_traits::Point_weight_d point_weight = + local_tr_traits.point_weight_d_object(); + typename Tr_traits::Point_drop_weight_d drop_w = + local_tr_traits.point_drop_weight_d_object(); + /*typename Tr_traits::Power_center_d power_center = + local_tr_traits.power_center_d_object();*/ // CJTODO + typename Get_functor::type power_center(local_tr_traits); + + // Estimate the tangent space + const Point ¢er_pt = m_points[i]; + if (!tangent_spaces_are_already_computed) + { +#ifdef CGAL_TC_EXPORT_NORMALS + m_tangent_spaces[i] = compute_tangent_space(center_pt, &m_normals[i]); +#else + m_tangent_spaces[i] = compute_tangent_space(center_pt); +#endif + } + + //*************************************************** + // Build a minimal triangulation in the tangent space + // (we only need the star of p) + //*************************************************** + + // Insert p + Tr_point wp = local_tr_traits.construct_weighted_point_d_object()( + local_tr_traits.construct_point_d_object()(Intrinsic_dimension, ORIGIN), + m_weights[i]); + center_vertex = local_tr.insert(wp); + center_vertex->data() = i; + + //const int NUM_NEIGHBORS = 150; + //KNS_range ins_range = m_points_ds.query_ANN(center_pt, NUM_NEIGHBORS); + INS_range ins_range = m_points_ds.query_incremental_ANN(center_pt); + + // While building the local triangulation, we keep the radius + // of the sphere "star sphere" centered at "center_vertex" + // and which contains all the + // circumspheres of the star of "center_vertex" + boost::optional star_sphere_squared_radius; + + // Insert points until we find a point which is outside "star shere" + for (INS_iterator nn_it = ins_range.begin() ; + nn_it != ins_range.end() ; + ++nn_it) + { + std::size_t neighbor_point_idx = nn_it->first; + + // ith point = p, which is already inserted + if (neighbor_point_idx != i) + { + const Point &neighbor_pt = m_points[neighbor_point_idx]; + FT neighbor_weight = m_weights[neighbor_point_idx]; + + if (star_sphere_squared_radius + && k_sqdist(center_pt, neighbor_pt) + > *star_sphere_squared_radius + SQ_HALF_SPARSITY) + break; + + Tr_point proj_pt = project_point_and_compute_weight( + neighbor_pt, neighbor_weight, center_pt, m_tangent_spaces[i], + local_tr_traits); + + Tr_vertex_handle vh = local_tr.insert_if_in_star(proj_pt, center_vertex); + //Tr_vertex_handle vh = local_tr.insert(proj_pt); + if (vh != Tr_vertex_handle()) + { + // CJTODO TEMP TEST + /*if (star_sphere_squared_radius + && k_sqdist(center_pt, neighbor_pt) + > *star_sphere_squared_radius + SQ_HALF_SPARSITY) + std::cout << "ARGGGGGGGH" << std::endl;*/ + + vh->data() = neighbor_point_idx; + + // Let's recompute star_sphere_squared_radius + if (local_tr.current_dimension() >= Intrinsic_dimension) + { + star_sphere_squared_radius = 0; + // Get the incident cells and look for the biggest circumsphere + std::vector incident_cells; + local_tr.incident_full_cells( + center_vertex, + std::back_inserter(incident_cells)); + for (typename std::vector::iterator cit = + incident_cells.begin(); cit != incident_cells.end(); ++cit) + { + Tr_full_cell_handle cell = *cit; + if (local_tr.is_infinite(cell)) + { + star_sphere_squared_radius = boost::none; + break; + } + else + { + std::vector cell_pts; + cell_pts.reserve(Intrinsic_dimension + 1); + // For each point p + for (int ii = 0 ; ii <= Intrinsic_dimension ; ++ii) + cell_pts.push_back(cell->vertex(ii)->point()); + + Tr_point c = power_center(cell_pts.begin(), cell_pts.end()); + FT sq_power_sphere_diam = 4*point_weight(c); + + if (!star_sphere_squared_radius + || sq_power_sphere_diam > *star_sphere_squared_radius) + { + star_sphere_squared_radius = sq_power_sphere_diam; + } + } + } + } + } + //std::cerr << star_sphere_squared_radius << std::endl; + } + } + + + // CJTODO DEBUG + //std::cerr << "\nChecking topology and geometry..." + // << (local_tr.is_valid(true) ? "OK.\n" : "Error.\n"); + // DEBUG: output the local mesh into an OFF file + //std::stringstream sstr; + //sstr << "data/local_tri_" << i << ".off"; + //std::ofstream off_stream_tr(sstr.str()); + //CGAL::export_triangulation_to_off(off_stream_tr, local_tr); + } + + Tangent_space_basis compute_tangent_space(const Point &p +#ifdef CGAL_TC_EXPORT_NORMALS + , Vector *p_normal +#endif + ) const + { + // Kernel functors + typename Kernel::Construct_vector_d constr_vec = + m_k.construct_vector_d_object(); + typename Kernel::Compute_coordinate_d coord = + m_k.compute_coordinate_d_object(); + typename Kernel::Squared_length_d sqlen = + m_k.squared_length_d_object(); + typename Kernel::Scaled_vector_d scaled_vec = + m_k.scaled_vector_d_object(); + typename Kernel::Scalar_product_d inner_pdct = + m_k.scalar_product_d_object(); + typename Kernel::Difference_of_vectors_d diff_vec = + m_k.difference_of_vectors_d_object(); + + KNS_range kns_range = m_points_ds.query_ANN( + p, NUM_POINTS_FOR_PCA, false); + + //******************************* PCA ************************************* + + const int amb_dim = m_k.point_dimension_d_object()(p); + // One row = one point + Eigen::MatrixXd mat_points(NUM_POINTS_FOR_PCA, amb_dim); + KNS_iterator nn_it = kns_range.begin(); + for (int j = 0 ; + j < NUM_POINTS_FOR_PCA && nn_it != kns_range.end() ; + ++j, ++nn_it) + { + for (int i = 0 ; i < amb_dim ; ++i) + mat_points(j, i) = CGAL::to_double(coord(m_points[nn_it->first], i)); + } + Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean(); + Eigen::MatrixXd cov = centered.adjoint() * centered; + Eigen::SelfAdjointEigenSolver eig(cov); + + // The eigenvectors are sorted in increasing order of their corresponding + // eigenvalues + Tangent_space_basis ts; + for (int i = amb_dim - 1 ; i >= amb_dim - Intrinsic_dimension ; --i) + { + ts.push_back(constr_vec( + amb_dim, + eig.eigenvectors().col(i).data(), + eig.eigenvectors().col(i).data() + amb_dim)); + } +#ifdef CGAL_TC_EXPORT_NORMALS + *p_normal = constr_vec( + amb_dim, + eig.eigenvectors().col(amb_dim - Intrinsic_dimension - 1).data(), + eig.eigenvectors().col(amb_dim - Intrinsic_dimension - 1).data() + amb_dim); +#endif + + //************************************************************************* + + //Vector n = m_k.point_to_vector_d_object()(p); + //n = scaled_vec(n, FT(1)/sqrt(sqlen(n))); + //std::cerr << "IP = " << inner_pdct(n, ts[0]) << " & " << inner_pdct(n, ts[1]) << std::endl; + + return compute_gram_schmidt_basis(ts, m_k); + + /* + // CJTODO: this is only for a sphere in R^3 + Vector t1(-p[1] - p[2], p[0], p[0]); + Vector t2(p[1] * t1[2] - p[2] * t1[1], + p[2] * t1[0] - p[0] * t1[2], + p[0] * t1[1] - p[1] * t1[0]); + + // Normalize t1 and t2 + typename Kernel::Squared_length_d sqlen = m_k.squared_length_d_object(); + typename Kernel::Scaled_vector_d scaled_vec = m_k.scaled_vector_d_object(); + + Tangent_space_basis ts; + ts.reserve(Intrinsic_dimension); + ts.push_back(scaled_vec(t1, FT(1)/CGAL::sqrt(sqlen(t1)))); + ts.push_back(scaled_vec(t2, FT(1)/CGAL::sqrt(sqlen(t2)))); + + return ts;*/ + + /* + // Alternative code (to be used later) + //Vector n = m_k.point_to_vector_d_object()(p); + //n = scaled_vec(n, FT(1)/sqrt(sqlen(n))); + //Vector t1(12., 15., 65.); + //Vector t2(32., 5., 85.); + //Tangent_space_basis ts; + //ts.reserve(Intrinsic_dimension); + //ts.push_back(diff_vec(t1, scaled_vec(n, inner_pdct(t1, n)))); + //ts.push_back(diff_vec(t2, scaled_vec(n, inner_pdct(t2, n)))); + //return compute_gram_schmidt_basis(ts, m_k); + */ + } + + // Project the point in the tangent space + // The weight will be the squared distance between p and the projection of p + Tr_bare_point project_point(const Point &p, const Point &origin, + const Tangent_space_basis &ts) const + { + typename Kernel::Scalar_product_d inner_pdct = + m_k.scalar_product_d_object(); + typename Kernel::Difference_of_points_d diff_points = + m_k.difference_of_points_d_object(); + + std::vector coords; + // Ambiant-space coords of the projected point + coords.reserve(Intrinsic_dimension); + for (std::size_t i = 0 ; i < Intrinsic_dimension ; ++i) + { + // Compute the inner product p * ts[i] + Vector v = diff_points(p, origin); + FT coord = inner_pdct(v, ts[i]); + coords.push_back(coord); + } + + return Tr_bare_point(Intrinsic_dimension, coords.begin(), coords.end()); + } + + // Project the point in the tangent space + // The weight will be the squared distance between p and the projection of p + Tr_point project_point_and_compute_weight( + const Point &p, FT w, const Point &origin, const Tangent_space_basis &ts, + const Tr_traits &tr_traits) const + { + const int point_dim = m_k.point_dimension_d_object()(p); + typename Kernel::Scalar_product_d inner_pdct = + m_k.scalar_product_d_object(); + typename Kernel::Difference_of_points_d diff_points = + m_k.difference_of_points_d_object(); + typename Kernel::Construct_cartesian_const_iterator_d ccci = + m_k.construct_cartesian_const_iterator_d_object(); + + Vector v = diff_points(p, origin); + + std::vector coords; + // Ambiant-space coords of the projected point + std::vector p_proj(ccci(origin), ccci(origin, 0)); + coords.reserve(Intrinsic_dimension); + for (std::size_t i = 0 ; i < Intrinsic_dimension ; ++i) + { + // Compute the inner product p * ts[i] + FT coord = inner_pdct(v, ts[i]); + coords.push_back(coord); + + // p_proj += coord * v; + for (int j = 0 ; j < point_dim ; ++j) + p_proj[j] += coord * ts[i][j]; + } + + Point projected_pt(point_dim, p_proj.begin(), p_proj.end()); + + return tr_traits.construct_weighted_point_d_object() + ( + tr_traits.construct_point_d_object()( + Intrinsic_dimension, coords.begin(), coords.end()), + w - m_k.squared_distance_d_object()(p, projected_pt) + ); + } + + // A simplex here is a local tri's full cell handle + bool is_simplex_consistent(Tr_full_cell_handle fch) + { + std::set c; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + { + std::size_t data = fch->vertex(i)->data(); + c.insert(data); + } + return is_simplex_consistent(c); + } + + // A simplex here is a list of point indices + bool is_simplex_consistent(std::set const& simplex) + { + // Check if the simplex is in the stars of all its vertices + std::set::const_iterator it_point_idx = simplex.begin(); + // For each point p of the simplex, we parse the incidents cells of p + // and we check if "simplex" is among them + for ( ; it_point_idx != simplex.end() ; ++it_point_idx) + { + std::size_t point_idx = *it_point_idx; + if (point_idx == std::numeric_limits::max()) + continue; + Triangulation const& tr = m_triangulations[point_idx].tr(); + Tr_vertex_handle center_vh = m_triangulations[point_idx].center_vertex(); + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + typename std::vector::const_iterator it_c = incident_cells.begin(); + typename std::vector::const_iterator it_c_end= incident_cells.end(); + // For each cell + bool found = false; + for ( ; !found && it_c != it_c_end ; ++it_c) + { + std::set cell; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + cell.insert((*it_c)->vertex(i)->data()); + if (cell == simplex) + found = true; + } + + if (!found) + return false; + } + + return true; + } + +#ifdef CGAL_LINKED_WITH_TBB + // Functor for try_to_solve_inconsistencies_in_a_local_triangulation function + class Try_to_solve_inconsistencies_in_a_local_triangulation + { + Tangential_complex & m_tc; + + public: + // Constructor + Try_to_solve_inconsistencies_in_a_local_triangulation( + Tangential_complex &tc) + : m_tc(tc) + {} + + // Constructor + Try_to_solve_inconsistencies_in_a_local_triangulation( + const Compute_tangent_triangulation &ctt) + : m_tc(ctt.m_tc) + {} + + // operator() + void operator()( const tbb::blocked_range& r ) const + { + for( size_t i = r.begin() ; i != r.end() ; ++i) + m_tc.try_to_solve_inconsistencies_in_a_local_triangulation(i); + } + }; +#endif // CGAL_LINKED_WITH_TBB + + void try_to_solve_inconsistencies_in_a_local_triangulation( + std::size_t tr_index) + { +#ifdef CGAL_LINKED_WITH_TBB + //Tr_mutex::scoped_lock lock(m_tr_mutexes[tr_index]); +#endif + + Triangulation const& tr = m_triangulations[tr_index].tr(); + Tr_vertex_handle center_vh = m_triangulations[tr_index].center_vertex(); + const Tr_traits &local_tr_traits = tr.geom_traits(); + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + typename std::vector::const_iterator it_c = + incident_cells.begin(); + typename std::vector::const_iterator it_c_end= + incident_cells.end(); + // For each cell + for ( ; it_c != it_c_end ; ++it_c) + { + +#ifdef CGAL_TC_ONLY_CHANGE_SIMPLEX_WEIGHTS + + std::set c; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + { + std::size_t data = (*it_c)->vertex(i)->data(); + c.insert(data); + } + + // Inconsistent? + if (!is_simplex_consistent(c)) + { + CGAL::Random rng; + for (std::set::iterator it=c.begin(); it!=c.end(); ++it) + { + m_weights[*it] = rng.get_double(0., SQ_HALF_SPARSITY); + } + +#if !defined(CGAL_TC_GLOBAL_REFRESH) + refresh_tangential_complex(); +#endif + + // We will try the other cells next time (incident_cells is not + // valid anymore here) + break; + } +#else + // Inconsistent? + if (!is_simplex_consistent(*it_c)) + { + // Get the k + 2 closest points + + /*int point_dim = m_k.point_dimension_d_object()(*m_points.begin()); + std::vector center(point_dim, FT(0)); + for (int i = 0 ; i < point_dim ; ++i) + { + for (int j = 0 ; j < Intrinsic_dimension + 1 ; ++j) + { + std::size_t data = (*it_c)->vertex(j)->data(); + const Point &p = m_points[data]; + center[i] += p[i]; + } + center[i] /= (Intrinsic_dimension + 1); + } + Point global_center(center.begin(), center.end());*/ + + std::vector simplex_pts; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + simplex_pts.push_back((*it_c)->vertex(i)->point()); + + //typename Tr_traits::Power_center_d power_center = + // local_tr_traits.power_center_d_object(); // CJTODO + typename Get_functor::type power_center(local_tr_traits); + typename Tr_traits::Compute_coordinate_d coord = + local_tr_traits.compute_coordinate_d_object(); + + typename Kernel::Translated_point_d k_transl = + m_k.translated_point_d_object(); + typename Kernel::Scaled_vector_d k_scaled_vec = + m_k.scaled_vector_d_object(); + + Tr_point local_center = power_center(simplex_pts.begin(), simplex_pts.end()); + Point global_center = m_points[tr_index]; + const Tangent_space_basis &tsb = m_tangent_spaces[tr_index]; + for (int i = 0 ; i < Intrinsic_dimension ; ++i) + { + global_center = k_transl( + global_center, + k_scaled_vec(tsb[i], coord(local_center, i))); + } + + KNS_range kns_range = m_points_ds.query_ANN( + global_center, + Intrinsic_dimension + 1 + + CGAL_TC_NUMBER_OF_ADDITIONNAL_PERTURBED_POINTS); + std::vector neighbors; + for (KNS_iterator nn_it = kns_range.begin() ; + nn_it != kns_range.end() ; + ++nn_it) + { + neighbors.push_back(nn_it->first); + } + + CGAL::Random rng; + for (std::vector::iterator it = neighbors.begin(); + it != neighbors.end() ; + ++it) + { + m_weights[*it] = rng.get_double(0., SQ_HALF_SPARSITY); + } + +#if !defined(CGAL_TC_GLOBAL_REFRESH) + refresh_tangential_complex(); +#endif + + // We will try the other cells next time (incident_cells is not + // valid anymore here) + break; + } + +#endif + } + } + + std::ostream &export_vertices_to_off( + std::ostream & os, std::size_t &num_vertices) + { + if (m_points.empty()) + { + num_vertices = 0; + return os; + } + + // If Intrinsic_dimension = 1, we output each point two times + // to be able to export each segment as a flat triangle with 3 different + // indices (otherwise, Meshlab detects degenerated simplices) + const int N = (Intrinsic_dimension == 1 ? 2 : 1); + const int ambient_dim = m_k.point_dimension_d_object()(*m_points.begin()); + + // Kernel functors + typename Kernel::Compute_coordinate_d coord = + m_k.compute_coordinate_d_object(); + + int num_coords = min(ambient_dim, 3); +#ifdef CGAL_TC_EXPORT_NORMALS + Normals::const_iterator it_n = m_normals.begin(); +#endif + typename Points::const_iterator it_p = m_points.begin(); + typename Points::const_iterator it_p_end = m_points.end(); + // For each point p + for ( ; it_p != it_p_end ; ++it_p) + { + for (int ii = 0 ; ii < N ; ++ii) + { + int i = 0; + for ( ; i < num_coords ; ++i) + os << CGAL::to_double(coord(*it_p, i)) << " "; + if (i == 2) + os << "0"; + +#ifdef CGAL_TC_EXPORT_NORMALS + for (i = 0 ; i < num_coords ; ++i) + os << " " << CGAL::to_double(coord(*it_n, i)); +#endif + os << std::endl; + } +#ifdef CGAL_TC_EXPORT_NORMALS + ++it_n; +#endif + } + + num_vertices = N*m_points.size(); + return os; + } + + std::ostream &export_simplices_to_off( + std::ostream & os, std::size_t &num_simplices, + bool color_inconsistencies = false, + std::set > const* excluded_simplices = NULL, + bool show_excluded_vertices_in_color = false) + { + // If Intrinsic_dimension = 1, each point is output two times + // (see export_vertices_to_off) + int factor = (Intrinsic_dimension == 1 ? 2 : 1); + int OFF_simplices_dim = + (Intrinsic_dimension == 1 ? 3 : Intrinsic_dimension + 1); + num_simplices = 0; + std::size_t num_inconsistent_simplices = 0; + typename Tr_container::const_iterator it_tr = m_triangulations.begin(); + typename Tr_container::const_iterator it_tr_end = m_triangulations.end(); + // For each triangulation + for ( ; it_tr != it_tr_end ; ++it_tr) + { + Triangulation const& tr = it_tr->tr(); + Tr_vertex_handle center_vh = it_tr->center_vertex(); + + if (tr.current_dimension() < Intrinsic_dimension) + continue; + + // Color for this star + std::stringstream color; + //color << rand()%256 << " " << 100+rand()%156 << " " << 100+rand()%156; + color << 128 << " " << 128 << " " << 128; + + std::vector incident_cells; + tr.incident_full_cells(center_vh, std::back_inserter(incident_cells)); + + typename std::vector::const_iterator it_c = incident_cells.begin(); + typename std::vector::const_iterator it_c_end= incident_cells.end(); + // For each cell + for ( ; it_c != it_c_end ; ++it_c) + { + if (tr.is_infinite(*it_c)) // Don't export infinite cells + continue; + + if (color_inconsistencies || excluded_simplices) + { + std::set c; + std::stringstream sstr_c; + std::size_t data; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + { + data = (*it_c)->vertex(i)->data(); + sstr_c << data*factor << " "; + c.insert(data); + } + // See export_vertices_to_off + if (Intrinsic_dimension == 1) + sstr_c << (data*factor + 1) << " "; + + bool excluded = + (excluded_simplices + && excluded_simplices->find(c) != excluded_simplices->end()); + + if (!excluded) + { + os << OFF_simplices_dim << " " << sstr_c.str() << " "; + if (color_inconsistencies && is_simplex_consistent(c)) + os << color.str(); + else + { + os << "255 0 0"; + ++num_inconsistent_simplices; + } + ++num_simplices; + } + else if (show_excluded_vertices_in_color) + { + os << OFF_simplices_dim << " " + << sstr_c.str() << " " + << "0 0 255"; + ++num_simplices; + } + } + else + { + os << OFF_simplices_dim << " "; + std::size_t data; + for (int i = 0 ; i < Intrinsic_dimension + 1 ; ++i) + { + data = (*it_c)->vertex(i)->data(); + os << data*factor << " "; + } + // See export_vertices_to_off + if (Intrinsic_dimension == 1) + os << (data*factor + 1) << " "; + + ++num_simplices; + } + + os << std::endl; + } + } + +#ifdef CGAL_TC_VERBOSE + std::cerr << std::endl + << "================================================" << std::endl + << "Export to OFF:\n" + << " * Number of vertices: " << m_points.size() << std::endl + << " * Total number of simplices in stars (incl. duplicates): " + << num_simplices << std::endl + << " * Number of inconsistent simplices in stars (incl. duplicates): " + << num_inconsistent_simplices << std::endl + << " * Percentage of inconsistencies: " + << (num_simplices > 0 ? + 100. * num_inconsistent_simplices / num_simplices : 0.) << "%" + << std::endl + << "================================================" << std::endl; +#endif + + return os; + } + +private: + const Kernel m_k; + + Points m_points; + Weights m_weights; + + Points_ds m_points_ds; + TS_container m_tangent_spaces; + Tr_container m_triangulations; // Contains the triangulations + // and their center vertex +#ifdef CGAL_LINKED_WITH_TBB + //std::vector m_tr_mutexes; +#endif +#ifdef CGAL_TC_EXPORT_NORMALS + Normals m_normals; +#endif + +}; // /class Tangential_complex + +} // end namespace CGAL + +#endif // TANGENTIAL_COMPLEX_H diff --git a/Tangential_complex/include/CGAL/Tangential_complex/Point_cloud.h b/Tangential_complex/include/CGAL/Tangential_complex/Point_cloud.h new file mode 100644 index 00000000000..1e466d65049 --- /dev/null +++ b/Tangential_complex/include/CGAL/Tangential_complex/Point_cloud.h @@ -0,0 +1,340 @@ +// 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 POINT_CLOUD_H +#define POINT_CLOUD_H + +#include +#include +#include + +#ifdef CGAL_TC_USE_NANOFLANN + +#include "nanoflann.hpp" + +#include +#include +#include + +namespace CGAL { +namespace Tangential_complex_ { + +// "dataset to kd-tree" adaptor class +template +class Point_cloud_adaptator +{ +public: + typedef typename Point_container_::value_type Point; + typedef typename CGAL::Kernel_traits::type Kernel; + typedef typename Kernel::FT FT; + + /// The constructor that sets the data set source + Point_cloud_adaptator(Point_container_ &points, Kernel const& k) + : m_points(points), m_k(k) + {} + + /// CRTP helper method + inline Point_container_ const& points() const + { + return m_points; + } + inline Point_container_& points() + { + return m_points; + } + + // Must return the number of data points + inline size_t kdtree_get_point_count() const + { + return m_points.size(); + } + + // Returns the distance between the vector "p1[0:size-1]" + // and the data point with index "idx_p2" stored in the class: + inline FT kdtree_distance( + const FT *p1, const size_t idx_p2, size_t size) const + { + Point sp(p1, p1 + size); + return m_k.squared_distance_d_object()(sp, points()[idx_p2]); + } + + // Returns the dim'th component of the idx'th point in the class: + // Since this is inlined and the "dim" argument is typically an + // immediate value, the "if/else's" are actually solved at compile time. + inline FT kdtree_get_pt(const size_t idx, int dim) const + { + return m_k.compute_coordinate_d_object()(points()[idx], dim); + } + + // Optional bounding-box computation: return false to default to a standard + // bbox computation loop. + // Return true if the BBOX was already computed by the class and returned + // in "bb" so it can be avoided to redo it again. + // Look at bb.size() to find out the expected dimensionality + // (e.g. 2 or 3 for point clouds) + template + bool kdtree_get_bbox(Bbox &bb) const + { + return false; + } + + Kernel const& kernel() const + { + return m_k; + } + +protected: + Point_container_& m_points; //!< A ref to the data set origin + Kernel const& m_k; //!< A const ref to the kernel + +}; + +template +class Point_cloud_data_structure +{ +public: + typedef typename Point_container_::value_type Point; + typedef typename K Kernel; + typedef typename Kernel::FT FT; + + static const int AMB_DIM = Ambient_dimension::value; // CJTODO: use Point_dimension_d or similar + + /// Constructor + Point_cloud_data_structure(Point_container_ &points, Kernel const& k) + : m_adaptor(points, k), + m_kd_tree(AMB_DIM, + m_adaptor, + nanoflann::KDTreeSingleIndexAdaptorParams(10 /* max leaf */) ) + { + m_kd_tree.buildIndex(); + } + + /*Point_container_ &points() + { + return m_adaptor.points(); + } + + const Point_container_ &points() const + { + return m_adaptor.points(); + }*/ + + void query_ANN(const Point &sp, + std::size_t k, + size_t *neighbor_indices, + FT *squared_distance) const + { + /*std::vector sp_vec( + m_adaptor.kernel().construct_cartesian_const_iterator_d_object()(sp), + m_adaptor.kernel().construct_cartesian_const_iterator_d_object()(sp, 0));*/ // CJTODO remettre + std::vector sp_vec; + for (int i = 0 ; i < 4 ; ++i) + sp_vec.push_back(sp[i]); + nanoflann::KNNResultSet result_set(k); + result_set.init(neighbor_indices, squared_distance); + m_kd_tree.findNeighbors(result_set, + &sp_vec[0], + nanoflann::SearchParams()); + + /*std::cout << "knnSearch(nn="<< num_results <<"): \n"; + for (int i = 0 ; i < num_results ; ++i) + { + std::cout << " * neighbor_indices = " << neighbor_indices [i] + << " (out_dist_sqr = " << squared_distance[i] << ")" + << std::endl; + }*/ + } + + void query_ball(const Point &sp, + const FT radius, + std::vector > &neighbors, + bool sort_output = true) + { + /*std::vector sp_vec( + m_adaptor.kernel().construct_cartesian_const_iterator_d_object()(sp), + m_adaptor.kernel().construct_cartesian_const_iterator_d_object()(sp, 0));*/ // CJTODO remettre + std::vector sp_vec; + for (int i = 0 ; i < 4 ; ++i) + sp_vec.push_back(sp[i]); + m_kd_tree.radiusSearch(&sp_vec[0], + radius, + neighbors, + nanoflann::SearchParams(32, 0.f, sort_output)); + + /*std::cout << "radiusSearch(num="<< neighbors.size() <<"): \n"; + for (const auto idx_and_dist : neighbors) + { + std::cout << " * neighbor_indices = " << idx_and_dist.first + << " (out_dist_sqr = " << idx_and_dist.second << ")" + << std::endl; + }*/ + } + +protected: + typedef Point_cloud_adaptator Adaptor; + typedef nanoflann::KDTreeSingleIndexAdaptor< + nanoflann::L2_Simple_Adaptor , + Adaptor, + AMB_DIM // dim + > Kd_tree; + + Adaptor m_adaptor; + Kd_tree m_kd_tree; +}; + +} // namespace Tangential_complex_ +} //namespace CGAL + +#else // !CGAL_TC_USE_NANOFLANN => use CGAL Spatial searching + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace CGAL { +namespace Tangential_complex_ { + +template +class Point_cloud_data_structure +{ +public: + typedef typename Point_container_::value_type Point; + typedef K Kernel; + typedef typename Kernel::FT FT; + + typedef CGAL::Search_traits< + FT, Point, + typename Kernel::Cartesian_const_iterator_d, + typename Kernel::Construct_cartesian_const_iterator_d> Traits_base; + // using a pointer as a special property map type + typedef CGAL::Search_traits_adapter< + std::ptrdiff_t, Point*, Traits_base> STraits; + + typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; + typedef typename K_neighbor_search::Tree Tree; + typedef typename K_neighbor_search::Distance Distance; + typedef typename K_neighbor_search::iterator KNS_iterator; + typedef K_neighbor_search KNS_range; + + typedef CGAL::Orthogonal_incremental_neighbor_search< + STraits, Distance, CGAL::Sliding_midpoint, Tree> + Incremental_neighbor_search; + typedef typename Incremental_neighbor_search::iterator INS_iterator; + typedef Incremental_neighbor_search INS_range; + + static const int AMB_DIM = Ambient_dimension::value; // CJTODO: use Point_dimension_d or similar + + /// Constructor + Point_cloud_data_structure(Point_container_ const& points) + : m_points(points), + m_tree(boost::counting_iterator(0), + boost::counting_iterator(points.size()), + typename Tree::Splitter(), + STraits((Point*)&(points[0])) ) + { } + + /// Constructor + Point_cloud_data_structure( + Point_container_ const& points, + std::size_t begin_idx, std::size_t past_the_end_idx) + : m_points(points), + m_tree( + boost::counting_iterator(begin_idx), + boost::counting_iterator(past_the_end_idx), + typename Tree::Splitter(), + STraits((Point*)&(points[0])) ) + { + } + + /*Point_container_ &points() + { + return m_points; + } + + const Point_container_ &points() const + { + return m_points; + }*/ + + // Be careful, this function invalidates the tree, + // which will be recomputed at the next query + void insert(std::ptrdiff_t point_idx) + { + m_tree.insert(point_idx); + } + + KNS_range query_ANN(const + Point &sp, + unsigned int k, + bool sorted = true) const + { + // Initialize the search structure, and search all N points + // Note that we need to pass the Distance explicitly since it needs to + // know the property map + K_neighbor_search search( + m_tree, + sp, + k, + FT(0), + true, + Distance_adapter >( + (Point*)&(m_points[0])), + sorted); + + return search; + } + + INS_range query_incremental_ANN(const Point &sp) const + { + // Initialize the search structure, and search all N points + // Note that we need to pass the Distance explicitly since it needs to + // know the property map + Incremental_neighbor_search search( + m_tree, + sp, + FT(0), + true, + Distance_adapter >( + (Point*)&(m_points[0])) ); + + return search; + } + +protected: + Point_container_ const& m_points; + Tree m_tree; +}; + +} // namespace Tangential_complex_ +} //namespace CGAL + +#endif // CGAL_TC_USE_NANOFLANN + +#endif // POINT_CLOUD_H diff --git a/Tangential_complex/include/CGAL/Tangential_complex/config.h b/Tangential_complex/include/CGAL/Tangential_complex/config.h new file mode 100644 index 00000000000..87e009191b7 --- /dev/null +++ b/Tangential_complex/include/CGAL/Tangential_complex/config.h @@ -0,0 +1,59 @@ +// 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_TC_CONFIG_H +#define CGAL_TC_CONFIG_H + +#include + +// Without TBB_USE_THREADING_TOOL Intel Inspector XE will report false +// positives in Intel TBB +// (http://software.intel.com/en-us/articles/compiler-settings-for-threading-error-analysis-in-intel-inspector-xe/) +#ifdef _DEBUG +# define TBB_USE_THREADING_TOOL +#endif + + +//========================= Debugging & profiling ============================= +#define CGAL_TC_PROFILING +#define CGAL_TC_VERBOSE + +// Solving inconsistencies: only change the weights of the inconsistent simplex +// or more? +//#define CGAL_TC_ONLY_CHANGE_SIMPLEX_WEIGHTS + +// Otherwise... +const int CGAL_TC_NUMBER_OF_ADDITIONNAL_PERTURBED_POINTS = 1; + +//========================= Strategy ========================================== +//#define CGAL_TC_USE_NANOFLANN +//#define CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER +#define CGAL_TC_GLOBAL_REFRESH +//#define CGAL_TC_ON_DEMAND_REFRESH // CJTODO: not implemented yet + // The idea is to perform a global refresh + some local refreshes, just + // for local tri where there are some inconsistencies + // But be careful: refreshing the TC may invalidate cells, so the + // incident cells have to be recomputed again + +//========================= Parameters ======================================== +const std::size_t NUM_POINTS_FOR_PCA = 50; +const double INPUT_SPARSITY = 0.05; + +#endif // CGAL_TC_CONFIG_H diff --git a/Tangential_complex/include/CGAL/Tangential_complex/nanoflann.hpp b/Tangential_complex/include/CGAL/Tangential_complex/nanoflann.hpp new file mode 100644 index 00000000000..8bedc45a2d3 --- /dev/null +++ b/Tangential_complex/include/CGAL/Tangential_complex/nanoflann.hpp @@ -0,0 +1,1449 @@ +/*********************************************************************** + * Software License Agreement (BSD License) + * + * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved. + * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved. + * Copyright 2011-2013 Jose Luis Blanco (joseluisblancoc@gmail.com). + * All rights reserved. + * + * THE BSD LICENSE + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + *************************************************************************/ + +#ifndef NANOFLANN_HPP_ +#define NANOFLANN_HPP_ + +#include +#include +#include +#include +#include // for fwrite() +#include // for fabs(),... +#include + +// Avoid conflicting declaration of min/max macros in windows headers +#if !defined(NOMINMAX) && (defined(_WIN32) || defined(_WIN32_) || defined(WIN32) || defined(_WIN64)) +# define NOMINMAX +# ifdef max +# undef max +# undef min +# endif +#endif + +namespace nanoflann +{ +/** @addtogroup nanoflann_grp nanoflann C++ library for ANN + * @{ */ + + /** Library version: 0xMmP (M=Major,m=minor,P=path) */ + #define NANOFLANN_VERSION 0x117 + + /** @addtogroup result_sets_grp Result set classes + * @{ */ + template + class KNNResultSet + { + IndexType * indices; + DistanceType* dists; + CountType capacity; + CountType count; + + public: + inline KNNResultSet(CountType capacity_) : capacity(capacity_), count(0) + { + } + + inline void init(IndexType* indices_, DistanceType* dists_) + { + indices = indices_; + dists = dists_; + count = 0; + dists[capacity-1] = (std::numeric_limits::max)(); + } + + inline CountType size() const + { + return count; + } + + inline bool full() const + { + return count == capacity; + } + + + inline void addPoint(DistanceType dist, IndexType index) + { + CountType i; + for (i=count; i>0; --i) { +#ifdef NANOFLANN_FIRST_MATCH // If defined and two poins have the same distance, the one with the lowest-index will be returned first. + if ( (dists[i-1]>dist) || ((dist==dists[i-1])&&(indices[i-1]>index)) ) { +#else + if (dists[i-1]>dist) { +#endif + if (i + class RadiusResultSet + { + public: + const DistanceType radius; + + std::vector >& m_indices_dists; + + inline RadiusResultSet(DistanceType radius_, std::vector >& indices_dists) : radius(radius_), m_indices_dists(indices_dists) + { + init(); + } + + inline ~RadiusResultSet() { } + + inline void init() { clear(); } + inline void clear() { m_indices_dists.clear(); } + + inline size_t size() const { return m_indices_dists.size(); } + + inline bool full() const { return true; } + + inline void addPoint(DistanceType dist, IndexType index) + { + if (dist 0 + */ + std::pair worst_item() const + { + if (m_indices_dists.empty()) throw std::runtime_error("Cannot invoke RadiusResultSet::worst_item() on an empty list of results."); + typedef typename std::vector >::const_iterator DistIt; + DistIt it = std::max_element(m_indices_dists.begin(), m_indices_dists.end()); + return *it; + } + }; + + /** operator "<" for std::sort() */ + struct IndexDist_Sorter + { + /** PairType will be typically: std::pair */ + template + inline bool operator()(const PairType &p1, const PairType &p2) const { + return p1.second < p2.second; + } + }; + + /** @} */ + + + /** @addtogroup loadsave_grp Load/save auxiliary functions + * @{ */ + template + void save_value(FILE* stream, const T& value, size_t count = 1) + { + fwrite(&value, sizeof(value),count, stream); + } + + template + void save_value(FILE* stream, const std::vector& value) + { + size_t size = value.size(); + fwrite(&size, sizeof(size_t), 1, stream); + fwrite(&value[0], sizeof(T), size, stream); + } + + template + void load_value(FILE* stream, T& value, size_t count = 1) + { + size_t read_cnt = fread(&value, sizeof(value), count, stream); + if (read_cnt != count) { + throw std::runtime_error("Cannot read from file"); + } + } + + + template + void load_value(FILE* stream, std::vector& value) + { + size_t size; + size_t read_cnt = fread(&size, sizeof(size_t), 1, stream); + if (read_cnt!=1) { + throw std::runtime_error("Cannot read from file"); + } + value.resize(size); + read_cnt = fread(&value[0], sizeof(T), size, stream); + if (read_cnt!=size) { + throw std::runtime_error("Cannot read from file"); + } + } + /** @} */ + + + /** @addtogroup metric_grp Metric (distance) classes + * @{ */ + + template inline T abs(T x) { return (x<0) ? -x : x; } + template<> inline int abs(int x) { return ::abs(x); } + template<> inline float abs(float x) { return fabsf(x); } + template<> inline double abs(double x) { return fabs(x); } + template<> inline long double abs(long double x) { return fabsl(x); } + + /** Manhattan distance functor (generic version, optimized for high-dimensionality data sets). + * Corresponding distance traits: nanoflann::metric_L1 + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DistanceType Type of distance variables (must be signed) (e.g. float, double, int64_t) + */ + template + struct L1_Adaptor + { + typedef T ElementType; + typedef _DistanceType DistanceType; + + const DataSource &data_source; + + L1_Adaptor(const DataSource &_data_source) : data_source(_data_source) { } + + inline DistanceType operator()(const T* a, const size_t b_idx, size_t size, DistanceType worst_dist = -1) const + { + DistanceType result = DistanceType(); + const T* last = a + size; + const T* lastgroup = last - 3; + size_t d = 0; + + /* Process 4 items with each loop for efficiency. */ + while (a < lastgroup) { + const DistanceType diff0 = nanoflann::abs(a[0] - data_source.kdtree_get_pt(b_idx,d++)); + const DistanceType diff1 = nanoflann::abs(a[1] - data_source.kdtree_get_pt(b_idx,d++)); + const DistanceType diff2 = nanoflann::abs(a[2] - data_source.kdtree_get_pt(b_idx,d++)); + const DistanceType diff3 = nanoflann::abs(a[3] - data_source.kdtree_get_pt(b_idx,d++)); + result += diff0 + diff1 + diff2 + diff3; + a += 4; + if ((worst_dist>0)&&(result>worst_dist)) { + return result; + } + } + /* Process last 0-3 components. Not needed for standard vector lengths. */ + while (a < last) { + result += nanoflann::abs( *a++ - data_source.kdtree_get_pt(b_idx,d++) ); + } + return result; + } + + template + inline DistanceType accum_dist(const U a, const V b, int ) const + { + return nanoflann::abs(a-b); + } + }; + + /** Squared Euclidean distance functor (generic version, optimized for high-dimensionality data sets). + * Corresponding distance traits: nanoflann::metric_L2 + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DistanceType Type of distance variables (must be signed) (e.g. float, double, int64_t) + */ + template + struct L2_Adaptor + { + typedef T ElementType; + typedef _DistanceType DistanceType; + + const DataSource &data_source; + + L2_Adaptor(const DataSource &_data_source) : data_source(_data_source) { } + + inline DistanceType operator()(const T* a, const size_t b_idx, size_t size, DistanceType worst_dist = -1) const + { + DistanceType result = DistanceType(); + const T* last = a + size; + const T* lastgroup = last - 3; + size_t d = 0; + + /* Process 4 items with each loop for efficiency. */ + while (a < lastgroup) { + const DistanceType diff0 = a[0] - data_source.kdtree_get_pt(b_idx,d++); + const DistanceType diff1 = a[1] - data_source.kdtree_get_pt(b_idx,d++); + const DistanceType diff2 = a[2] - data_source.kdtree_get_pt(b_idx,d++); + const DistanceType diff3 = a[3] - data_source.kdtree_get_pt(b_idx,d++); + result += diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3; + a += 4; + if ((worst_dist>0)&&(result>worst_dist)) { + return result; + } + } + /* Process last 0-3 components. Not needed for standard vector lengths. */ + while (a < last) { + const DistanceType diff0 = *a++ - data_source.kdtree_get_pt(b_idx,d++); + result += diff0 * diff0; + } + return result; + } + + template + inline DistanceType accum_dist(const U a, const V b, int ) const + { + return (a-b)*(a-b); + } + }; + + /** Squared Euclidean distance functor (suitable for low-dimensionality datasets, like 2D or 3D point clouds) + * Corresponding distance traits: nanoflann::metric_L2_Simple + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DistanceType Type of distance variables (must be signed) (e.g. float, double, int64_t) + */ + template + struct L2_Simple_Adaptor + { + typedef T ElementType; + typedef _DistanceType DistanceType; + + const DataSource &data_source; + + L2_Simple_Adaptor(const DataSource &_data_source) : data_source(_data_source) { } + + inline DistanceType operator()(const T* a, const size_t b_idx, size_t size) const { + return data_source.kdtree_distance(a,b_idx,size); + } + + template + inline DistanceType accum_dist(const U a, const V b, int ) const + { + return (a-b)*(a-b); + } + }; + + /** Metaprogramming helper traits class for the L1 (Manhattan) metric */ + struct metric_L1 { + template + struct traits { + typedef L1_Adaptor distance_t; + }; + }; + /** Metaprogramming helper traits class for the L2 (Euclidean) metric */ + struct metric_L2 { + template + struct traits { + typedef L2_Adaptor distance_t; + }; + }; + /** Metaprogramming helper traits class for the L2_simple (Euclidean) metric */ + struct metric_L2_Simple { + template + struct traits { + typedef L2_Simple_Adaptor distance_t; + }; + }; + + /** @} */ + + + + /** @addtogroup param_grp Parameter structs + * @{ */ + + /** Parameters (see http://code.google.com/p/nanoflann/ for help choosing the parameters) + */ + struct KDTreeSingleIndexAdaptorParams + { + KDTreeSingleIndexAdaptorParams(size_t _leaf_max_size = 10, int dim_ = -1) : + leaf_max_size(_leaf_max_size), dim(dim_) + {} + + size_t leaf_max_size; + int dim; + }; + + /** Search options for KDTreeSingleIndexAdaptor::findNeighbors() */ + struct SearchParams + { + /** Note: The first argument (checks_IGNORED_) is ignored, but kept for compatibility with the FLANN interface */ + SearchParams(int checks_IGNORED_ = 32, float eps_ = 0, bool sorted_ = true ) : + checks(checks_IGNORED_), eps(eps_), sorted(sorted_) {} + + int checks; //!< Ignored parameter (Kept for compatibility with the FLANN interface). + float eps; //!< search for eps-approximate neighbours (default: 0) + bool sorted; //!< only for radius search, require neighbours sorted by distance (default: true) + }; + /** @} */ + + + /** @addtogroup memalloc_grp Memory allocation + * @{ */ + + /** + * Allocates (using C's malloc) a generic type T. + * + * Params: + * count = number of instances to allocate. + * Returns: pointer (of type T*) to memory buffer + */ + template + inline T* allocate(size_t count = 1) + { + T* mem = (T*) ::malloc(sizeof(T)*count); + return mem; + } + + + /** + * Pooled storage allocator + * + * The following routines allow for the efficient allocation of storage in + * small chunks from a specified pool. Rather than allowing each structure + * to be freed individually, an entire pool of storage is freed at once. + * This method has two advantages over just using malloc() and free(). First, + * it is far more efficient for allocating small objects, as there is + * no overhead for remembering all the information needed to free each + * object or consolidating fragmented memory. Second, the decision about + * how long to keep an object is made at the time of allocation, and there + * is no need to track down all the objects to free them. + * + */ + + const size_t WORDSIZE=16; + const size_t BLOCKSIZE=8192; + + class PooledAllocator + { + /* We maintain memory alignment to word boundaries by requiring that all + allocations be in multiples of the machine wordsize. */ + /* Size of machine word in bytes. Must be power of 2. */ + /* Minimum number of bytes requested at a time from the system. Must be multiple of WORDSIZE. */ + + + size_t remaining; /* Number of bytes left in current block of storage. */ + void* base; /* Pointer to base of current block of storage. */ + void* loc; /* Current location in block to next allocate memory. */ + size_t blocksize; + + void internal_init() + { + remaining = 0; + base = NULL; + usedMemory = 0; + wastedMemory = 0; + } + + public: + size_t usedMemory; + size_t wastedMemory; + + /** + Default constructor. Initializes a new pool. + */ + PooledAllocator(const size_t blocksize_ = BLOCKSIZE) : blocksize(blocksize_) { + internal_init(); + } + + /** + * Destructor. Frees all the memory allocated in this pool. + */ + ~PooledAllocator() { + free_all(); + } + + /** Frees all allocated memory chunks */ + void free_all() + { + while (base != NULL) { + void *prev = *((void**) base); /* Get pointer to prev block. */ + ::free(base); + base = prev; + } + internal_init(); + } + + /** + * Returns a pointer to a piece of new memory of the given size in bytes + * allocated from the pool. + */ + void* malloc(const size_t req_size) + { + /* Round size up to a multiple of wordsize. The following expression + only works for WORDSIZE that is a power of 2, by masking last bits of + incremented size to zero. + */ + const size_t size = (req_size + (WORDSIZE - 1)) & ~(WORDSIZE - 1); + + /* Check whether a new block must be allocated. Note that the first word + of a block is reserved for a pointer to the previous block. + */ + if (size > remaining) { + + wastedMemory += remaining; + + /* Allocate new storage. */ + const size_t blocksize = (size + sizeof(void*) + (WORDSIZE-1) > BLOCKSIZE) ? + size + sizeof(void*) + (WORDSIZE-1) : BLOCKSIZE; + + // use the standard C malloc to allocate memory + void* m = ::malloc(blocksize); + if (!m) { + fprintf(stderr,"Failed to allocate memory.\n"); + return NULL; + } + + /* Fill first word of new block with pointer to previous block. */ + ((void**) m)[0] = base; + base = m; + + size_t shift = 0; + //int size_t = (WORDSIZE - ( (((size_t)m) + sizeof(void*)) & (WORDSIZE-1))) & (WORDSIZE-1); + + remaining = blocksize - sizeof(void*) - shift; + loc = ((char*)m + sizeof(void*) + shift); + } + void* rloc = loc; + loc = (char*)loc + size; + remaining -= size; + + usedMemory += size; + + return rloc; + } + + /** + * Allocates (using this pool) a generic type T. + * + * Params: + * count = number of instances to allocate. + * Returns: pointer (of type T*) to memory buffer + */ + template + T* allocate(const size_t count = 1) + { + T* mem = (T*) this->malloc(sizeof(T)*count); + return mem; + } + + }; + /** @} */ + + /** @addtogroup nanoflann_metaprog_grp Auxiliary metaprogramming stuff + * @{ */ + + // ---------------- CArray ------------------------- + /** A STL container (as wrapper) for arrays of constant size defined at compile time (class imported from the MRPT project) + * This code is an adapted version from Boost, modifed for its integration + * within MRPT (JLBC, Dec/2009) (Renamed array -> CArray to avoid possible potential conflicts). + * See + * http://www.josuttis.com/cppcode + * for details and the latest version. + * See + * http://www.boost.org/libs/array for Documentation. + * for documentation. + * + * (C) Copyright Nicolai M. Josuttis 2001. + * Permission to copy, use, modify, sell and distribute this software + * is granted provided this copyright notice appears in all copies. + * This software is provided "as is" without express or implied + * warranty, and with no claim as to its suitability for any purpose. + * + * 29 Jan 2004 - minor fixes (Nico Josuttis) + * 04 Dec 2003 - update to synch with library TR1 (Alisdair Meredith) + * 23 Aug 2002 - fix for Non-MSVC compilers combined with MSVC libraries. + * 05 Aug 2001 - minor update (Nico Josuttis) + * 20 Jan 2001 - STLport fix (Beman Dawes) + * 29 Sep 2000 - Initial Revision (Nico Josuttis) + * + * Jan 30, 2004 + */ + template + class CArray { + public: + T elems[N]; // fixed-size array of elements of type T + + public: + // type definitions + typedef T value_type; + typedef T* iterator; + typedef const T* const_iterator; + typedef T& reference; + typedef const T& const_reference; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + + // iterator support + inline iterator begin() { return elems; } + inline const_iterator begin() const { return elems; } + inline iterator end() { return elems+N; } + inline const_iterator end() const { return elems+N; } + + // reverse iterator support +#if !defined(BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION) && !defined(BOOST_MSVC_STD_ITERATOR) && !defined(BOOST_NO_STD_ITERATOR_TRAITS) + typedef std::reverse_iterator reverse_iterator; + typedef std::reverse_iterator const_reverse_iterator; +#elif defined(_MSC_VER) && (_MSC_VER == 1300) && defined(BOOST_DINKUMWARE_STDLIB) && (BOOST_DINKUMWARE_STDLIB == 310) + // workaround for broken reverse_iterator in VC7 + typedef std::reverse_iterator > reverse_iterator; + typedef std::reverse_iterator > const_reverse_iterator; +#else + // workaround for broken reverse_iterator implementations + typedef std::reverse_iterator reverse_iterator; + typedef std::reverse_iterator const_reverse_iterator; +#endif + + reverse_iterator rbegin() { return reverse_iterator(end()); } + const_reverse_iterator rbegin() const { return const_reverse_iterator(end()); } + reverse_iterator rend() { return reverse_iterator(begin()); } + const_reverse_iterator rend() const { return const_reverse_iterator(begin()); } + // operator[] + inline reference operator[](size_type i) { return elems[i]; } + inline const_reference operator[](size_type i) const { return elems[i]; } + // at() with range check + reference at(size_type i) { rangecheck(i); return elems[i]; } + const_reference at(size_type i) const { rangecheck(i); return elems[i]; } + // front() and back() + reference front() { return elems[0]; } + const_reference front() const { return elems[0]; } + reference back() { return elems[N-1]; } + const_reference back() const { return elems[N-1]; } + // size is constant + static inline size_type size() { return N; } + static bool empty() { return false; } + static size_type max_size() { return N; } + enum { static_size = N }; + /** This method has no effects in this class, but raises an exception if the expected size does not match */ + inline void resize(const size_t nElements) { if (nElements!=N) throw std::logic_error("Try to change the size of a CArray."); } + // swap (note: linear complexity in N, constant for given instantiation) + void swap (CArray& y) { std::swap_ranges(begin(),end(),y.begin()); } + // direct access to data (read-only) + const T* data() const { return elems; } + // use array as C array (direct read/write access to data) + T* data() { return elems; } + // assignment with type conversion + template CArray& operator= (const CArray& rhs) { + std::copy(rhs.begin(),rhs.end(), begin()); + return *this; + } + // assign one value to all elements + inline void assign (const T& value) { for (size_t i=0;i= size()) { throw std::out_of_range("CArray<>: index out of range"); } } + }; // end of CArray + + /** Used to declare fixed-size arrays when DIM>0, dynamically-allocated vectors when DIM=-1. + * Fixed size version for a generic DIM: + */ + template + struct array_or_vector_selector + { + typedef CArray container_t; + }; + /** Dynamic size version */ + template + struct array_or_vector_selector<-1,T> { + typedef std::vector container_t; + }; + /** @} */ + + /** @addtogroup kdtrees_grp KD-tree classes and adaptors + * @{ */ + + /** kd-tree index + * + * Contains the k-d trees and other information for indexing a set of points + * for nearest-neighbor matching. + * + * The class "DatasetAdaptor" must provide the following interface (can be non-virtual, inlined methods): + * + * \code + * // Must return the number of data points + * inline size_t kdtree_get_point_count() const { ... } + * + * // Must return the Euclidean (L2) distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class: + * inline DistanceType kdtree_distance(const T *p1, const size_t idx_p2,size_t size) const { ... } + * + * // Must return the dim'th component of the idx'th point in the class: + * inline T kdtree_get_pt(const size_t idx, int dim) const { ... } + * + * // Optional bounding-box computation: return false to default to a standard bbox computation loop. + * // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it again. + * // Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds) + * template + * bool kdtree_get_bbox(BBOX &bb) const + * { + * bb[0].low = ...; bb[0].high = ...; // 0th dimension limits + * bb[1].low = ...; bb[1].high = ...; // 1st dimension limits + * ... + * return true; + * } + * + * \endcode + * + * \tparam IndexType Will be typically size_t or int + */ + template + class KDTreeSingleIndexAdaptor + { + public: + typedef typename Distance::ElementType ElementType; + typedef typename Distance::DistanceType DistanceType; + protected: + + /** + * Array of indices to vectors in the dataset. + */ + std::vector vind; + + size_t m_leaf_max_size; + + + /** + * The dataset used by this index + */ + const DatasetAdaptor &dataset; //!< The source of our data + + const KDTreeSingleIndexAdaptorParams index_params; + + size_t m_size; + int dim; //!< Dimensionality of each data point + + + /*--------------------- Internal Data Structures --------------------------*/ + struct Node + { + union { + struct + { + /** + * Indices of points in leaf node + */ + IndexType left, right; + } lr; + struct + { + /** + * Dimension used for subdivision. + */ + int divfeat; + /** + * The values used for subdivision. + */ + DistanceType divlow, divhigh; + } sub; + }; + /** + * The child nodes. + */ + Node* child1, * child2; + }; + typedef Node* NodePtr; + + + struct Interval + { + ElementType low, high; + }; + + /** Define "BoundingBox" as a fixed-size or variable-size container depending on "DIM" */ + typedef typename array_or_vector_selector::container_t BoundingBox; + + /** Define "distance_vector_t" as a fixed-size or variable-size container depending on "DIM" */ + typedef typename array_or_vector_selector::container_t distance_vector_t; + + /** This record represents a branch point when finding neighbors in + the tree. It contains a record of the minimum distance to the query + point, as well as the node at which the search resumes. + */ + template + struct BranchStruct + { + T node; /* Tree node at which search resumes */ + DistanceType mindist; /* Minimum distance to query for all nodes below. */ + + BranchStruct() {} + BranchStruct(const T& aNode, DistanceType dist) : node(aNode), mindist(dist) {} + + inline bool operator<(const BranchStruct& rhs) const + { + return mindist BranchSt; + typedef BranchSt* Branch; + + BoundingBox root_bbox; + + /** + * Pooled memory allocator. + * + * Using a pooled memory allocator is more efficient + * than allocating memory directly when there is a large + * number small of memory allocations. + */ + PooledAllocator pool; + + public: + + Distance distance; + + /** + * KDTree constructor + * + * Params: + * inputData = dataset with the input features + * params = parameters passed to the kdtree algorithm (see http://code.google.com/p/nanoflann/ for help choosing the parameters) + */ + KDTreeSingleIndexAdaptor(const int dimensionality, const DatasetAdaptor& inputData, const KDTreeSingleIndexAdaptorParams& params = KDTreeSingleIndexAdaptorParams() ) : + dataset(inputData), index_params(params), root_node(NULL), distance(inputData) + { + m_size = dataset.kdtree_get_point_count(); + dim = dimensionality; + if (DIM>0) dim=DIM; + else { + if (params.dim>0) dim = params.dim; + } + m_leaf_max_size = params.leaf_max_size; + + // Create a permutable array of indices to the input vectors. + init_vind(); + } + + /** + * Standard destructor + */ + ~KDTreeSingleIndexAdaptor() + { + } + + /** Frees the previously-built index. Automatically called within buildIndex(). */ + void freeIndex() + { + pool.free_all(); + root_node=NULL; + } + + /** + * Builds the index + */ + void buildIndex() + { + init_vind(); + computeBoundingBox(root_bbox); + freeIndex(); + root_node = divideTree(0, m_size, root_bbox ); // construct the tree + } + + /** + * Returns size of index. + */ + size_t size() const + { + return m_size; + } + + /** + * Returns the length of an index feature. + */ + size_t veclen() const + { + return static_cast(DIM>0 ? DIM : dim); + } + + /** + * Computes the inde memory usage + * Returns: memory used by the index + */ + size_t usedMemory() const + { + return pool.usedMemory+pool.wastedMemory+dataset.kdtree_get_point_count()*sizeof(IndexType); // pool memory and vind array memory + } + + /** \name Query methods + * @{ */ + + /** + * Find set of nearest neighbors to vec[0:dim-1]. Their indices are stored inside + * the result object. + * + * Params: + * result = the result object in which the indices of the nearest-neighbors are stored + * vec = the vector for which to search the nearest neighbors + * + * \tparam RESULTSET Should be any ResultSet + * \sa knnSearch, radiusSearch + */ + template + void findNeighbors(RESULTSET& result, const ElementType* vec, const SearchParams& searchParams) const + { + assert(vec); + if (!root_node) throw std::runtime_error("[nanoflann] findNeighbors() called before building the index."); + float epsError = 1+searchParams.eps; + + distance_vector_t dists; // fixed or variable-sized container (depending on DIM) + dists.assign((DIM>0 ? DIM : dim) ,0); // Fill it with zeros. + DistanceType distsq = computeInitialDistances(vec, dists); + searchLevel(result, vec, root_node, distsq, dists, epsError); // "count_leaf" parameter removed since was neither used nor returned to the user. + } + + /** + * Find the "num_closest" nearest neighbors to the \a query_point[0:dim-1]. Their indices are stored inside + * the result object. + * \sa radiusSearch, findNeighbors + * \note nChecks_IGNORED is ignored but kept for compatibility with the original FLANN interface. + */ + inline void knnSearch(const ElementType *query_point, const size_t num_closest, IndexType *out_indices, DistanceType *out_distances_sq, const int nChecks_IGNORED = 10) const + { + nanoflann::KNNResultSet resultSet(num_closest); + resultSet.init(out_indices, out_distances_sq); + this->findNeighbors(resultSet, query_point, nanoflann::SearchParams()); + } + + /** + * Find all the neighbors to \a query_point[0:dim-1] within a maximum radius. + * The output is given as a vector of pairs, of which the first element is a point index and the second the corresponding distance. + * Previous contents of \a IndicesDists are cleared. + * + * If searchParams.sorted==true, the output list is sorted by ascending distances. + * + * For a better performance, it is advisable to do a .reserve() on the vector if you have any wild guess about the number of expected matches. + * + * \sa knnSearch, findNeighbors + * \return The number of points within the given radius (i.e. indices.size() or dists.size() ) + */ + size_t radiusSearch(const ElementType *query_point,const DistanceType radius, std::vector >& IndicesDists, const SearchParams& searchParams) const + { + RadiusResultSet resultSet(radius,IndicesDists); + this->findNeighbors(resultSet, query_point, searchParams); + + if (searchParams.sorted) + std::sort(IndicesDists.begin(),IndicesDists.end(), IndexDist_Sorter() ); + + return resultSet.size(); + } + + /** @} */ + + private: + /** Make sure the auxiliary list \a vind has the same size than the current dataset, and re-generate if size has changed. */ + void init_vind() + { + // Create a permutable array of indices to the input vectors. + m_size = dataset.kdtree_get_point_count(); + if (vind.size()!=m_size) vind.resize(m_size); + for (size_t i = 0; i < m_size; i++) vind[i] = i; + } + + /// Helper accessor to the dataset points: + inline ElementType dataset_get(size_t idx, int component) const { + return dataset.kdtree_get_pt(idx,component); + } + + + void save_tree(FILE* stream, NodePtr tree) + { + save_value(stream, *tree); + if (tree->child1!=NULL) { + save_tree(stream, tree->child1); + } + if (tree->child2!=NULL) { + save_tree(stream, tree->child2); + } + } + + + void load_tree(FILE* stream, NodePtr& tree) + { + tree = pool.allocate(); + load_value(stream, *tree); + if (tree->child1!=NULL) { + load_tree(stream, tree->child1); + } + if (tree->child2!=NULL) { + load_tree(stream, tree->child2); + } + } + + + void computeBoundingBox(BoundingBox& bbox) + { + bbox.resize((DIM>0 ? DIM : dim)); + if (dataset.kdtree_get_bbox(bbox)) + { + // Done! It was implemented in derived class + } + else + { + for (int i=0; i<(DIM>0 ? DIM : dim); ++i) { + bbox[i].low = + bbox[i].high = dataset_get(0,i); + } + const size_t N = dataset.kdtree_get_point_count(); + for (size_t k=1; k0 ? DIM : dim); ++i) { + if (dataset_get(k,i)bbox[i].high) bbox[i].high = dataset_get(k,i); + } + } + } + } + + + /** + * Create a tree node that subdivides the list of vecs from vind[first] + * to vind[last]. The routine is called recursively on each sublist. + * Place a pointer to this new tree node in the location pTree. + * + * Params: pTree = the new node to create + * first = index of the first vector + * last = index of the last vector + */ + NodePtr divideTree(const IndexType left, const IndexType right, BoundingBox& bbox) + { + NodePtr node = pool.allocate(); // allocate memory + + /* If too few exemplars remain, then make this a leaf node. */ + if ( (right-left) <= m_leaf_max_size) { + node->child1 = node->child2 = NULL; /* Mark as leaf node. */ + node->lr.left = left; + node->lr.right = right; + + // compute bounding-box of leaf points + for (int i=0; i<(DIM>0 ? DIM : dim); ++i) { + bbox[i].low = dataset_get(vind[left],i); + bbox[i].high = dataset_get(vind[left],i); + } + for (IndexType k=left+1; k0 ? DIM : dim); ++i) { + if (bbox[i].low>dataset_get(vind[k],i)) bbox[i].low=dataset_get(vind[k],i); + if (bbox[i].highsub.divfeat = cutfeat; + + BoundingBox left_bbox(bbox); + left_bbox[cutfeat].high = cutval; + node->child1 = divideTree(left, left+idx, left_bbox); + + BoundingBox right_bbox(bbox); + right_bbox[cutfeat].low = cutval; + node->child2 = divideTree(left+idx, right, right_bbox); + + node->sub.divlow = left_bbox[cutfeat].high; + node->sub.divhigh = right_bbox[cutfeat].low; + + for (int i=0; i<(DIM>0 ? DIM : dim); ++i) { + bbox[i].low = std::min(left_bbox[i].low, right_bbox[i].low); + bbox[i].high = std::max(left_bbox[i].high, right_bbox[i].high); + } + } + + return node; + } + + void computeMinMax(IndexType* ind, IndexType count, int element, ElementType& min_elem, ElementType& max_elem) + { + min_elem = dataset_get(ind[0],element); + max_elem = dataset_get(ind[0],element); + for (IndexType i=1; imax_elem) max_elem = val; + } + } + + void middleSplit(IndexType* ind, IndexType count, IndexType& index, int& cutfeat, DistanceType& cutval, const BoundingBox& bbox) + { + // find the largest span from the approximate bounding box + ElementType max_span = bbox[0].high-bbox[0].low; + cutfeat = 0; + cutval = (bbox[0].high+bbox[0].low)/2; + for (int i=1; i<(DIM>0 ? DIM : dim); ++i) { + ElementType span = bbox[i].low-bbox[i].low; + if (span>max_span) { + max_span = span; + cutfeat = i; + cutval = (bbox[i].high+bbox[i].low)/2; + } + } + + // compute exact span on the found dimension + ElementType min_elem, max_elem; + computeMinMax(ind, count, cutfeat, min_elem, max_elem); + cutval = (min_elem+max_elem)/2; + max_span = max_elem - min_elem; + + // check if a dimension of a largest span exists + size_t k = cutfeat; + for (size_t i=0; i<(DIM>0 ? DIM : dim); ++i) { + if (i==k) continue; + ElementType span = bbox[i].high-bbox[i].low; + if (span>max_span) { + computeMinMax(ind, count, i, min_elem, max_elem); + span = max_elem - min_elem; + if (span>max_span) { + max_span = span; + cutfeat = i; + cutval = (min_elem+max_elem)/2; + } + } + } + IndexType lim1, lim2; + planeSplit(ind, count, cutfeat, cutval, lim1, lim2); + + if (lim1>count/2) index = lim1; + else if (lim2(0.00001); + ElementType max_span = bbox[0].high-bbox[0].low; + for (int i=1; i<(DIM>0 ? DIM : dim); ++i) { + ElementType span = bbox[i].high-bbox[i].low; + if (span>max_span) { + max_span = span; + } + } + ElementType max_spread = -1; + cutfeat = 0; + for (int i=0; i<(DIM>0 ? DIM : dim); ++i) { + ElementType span = bbox[i].high-bbox[i].low; + if (span>(1-EPS)*max_span) { + ElementType min_elem, max_elem; + computeMinMax(ind, count, cutfeat, min_elem, max_elem); + ElementType spread = max_elem-min_elem;; + if (spread>max_spread) { + cutfeat = i; + max_spread = spread; + } + } + } + // split in the middle + DistanceType split_val = (bbox[cutfeat].low+bbox[cutfeat].high)/2; + ElementType min_elem, max_elem; + computeMinMax(ind, count, cutfeat, min_elem, max_elem); + + if (split_valmax_elem) cutval = max_elem; + else cutval = split_val; + + IndexType lim1, lim2; + planeSplit(ind, count, cutfeat, cutval, lim1, lim2); + + if (lim1>count/2) index = lim1; + else if (lim2cutval + */ + void planeSplit(IndexType* ind, const IndexType count, int cutfeat, DistanceType cutval, IndexType& lim1, IndexType& lim2) + { + /* Move vector indices for left subtree to front of list. */ + IndexType left = 0; + IndexType right = count-1; + for (;; ) { + while (left<=right && dataset_get(ind[left],cutfeat)=cutval) --right; + if (left>right || !right) break; // "!right" was added to support unsigned Index types + std::swap(ind[left], ind[right]); + ++left; + --right; + } + /* If either list is empty, it means that all remaining features + * are identical. Split in the middle to maintain a balanced tree. + */ + lim1 = left; + right = count-1; + for (;; ) { + while (left<=right && dataset_get(ind[left],cutfeat)<=cutval) ++left; + while (right && left<=right && dataset_get(ind[right],cutfeat)>cutval) --right; + if (left>right || !right) break; // "!right" was added to support unsigned Index types + std::swap(ind[left], ind[right]); + ++left; + --right; + } + lim2 = left; + } + + DistanceType computeInitialDistances(const ElementType* vec, distance_vector_t& dists) const + { + assert(vec); + DistanceType distsq = 0.0; + + for (int i = 0; i < (DIM>0 ? DIM : dim); ++i) { + if (vec[i] < root_bbox[i].low) { + dists[i] = distance.accum_dist(vec[i], root_bbox[i].low, i); + distsq += dists[i]; + } + if (vec[i] > root_bbox[i].high) { + dists[i] = distance.accum_dist(vec[i], root_bbox[i].high, i); + distsq += dists[i]; + } + } + + return distsq; + } + + /** + * Performs an exact search in the tree starting from a node. + * \tparam RESULTSET Should be any ResultSet + */ + template + void searchLevel(RESULTSET& result_set, const ElementType* vec, const NodePtr node, DistanceType mindistsq, + distance_vector_t& dists, const float epsError) const + { + /* If this is a leaf node, then do check and return. */ + if ((node->child1 == NULL)&&(node->child2 == NULL)) { + //count_leaf += (node->lr.right-node->lr.left); // Removed since was neither used nor returned to the user. + DistanceType worst_dist = result_set.worstDist(); + for (IndexType i=node->lr.left; ilr.right; ++i) { + const IndexType index = vind[i];// reorder... : i; + DistanceType dist = distance(vec, index, (DIM>0 ? DIM : dim)); + if (distsub.divfeat; + ElementType val = vec[idx]; + DistanceType diff1 = val - node->sub.divlow; + DistanceType diff2 = val - node->sub.divhigh; + + NodePtr bestChild; + NodePtr otherChild; + DistanceType cut_dist; + if ((diff1+diff2)<0) { + bestChild = node->child1; + otherChild = node->child2; + cut_dist = distance.accum_dist(val, node->sub.divhigh, idx); + } + else { + bestChild = node->child2; + otherChild = node->child1; + cut_dist = distance.accum_dist( val, node->sub.divlow, idx); + } + + /* Call recursively to search next level down. */ + searchLevel(result_set, vec, bestChild, mindistsq, dists, epsError); + + DistanceType dst = dists[idx]; + mindistsq = mindistsq + cut_dist - dst; + dists[idx] = cut_dist; + if (mindistsq*epsError<=result_set.worstDist()) { + searchLevel(result_set, vec, otherChild, mindistsq, dists, epsError); + } + dists[idx] = dst; + } + + public: + /** Stores the index in a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so when loading the index object it must be constructed associated to the same source of data points used while building it. + * See the example: examples/saveload_example.cpp + * \sa loadIndex */ + void saveIndex(FILE* stream) + { + save_value(stream, m_size); + save_value(stream, dim); + save_value(stream, root_bbox); + save_value(stream, m_leaf_max_size); + save_value(stream, vind); + save_tree(stream, root_node); + } + + /** Loads a previous index from a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so the index object must be constructed associated to the same source of data points used while building the index. + * See the example: examples/saveload_example.cpp + * \sa loadIndex */ + void loadIndex(FILE* stream) + { + load_value(stream, m_size); + load_value(stream, dim); + load_value(stream, root_bbox); + load_value(stream, m_leaf_max_size); + load_value(stream, vind); + load_tree(stream, root_node); + } + + }; // class KDTree + + + /** A simple KD-tree adaptor for working with data directly stored in an Eigen Matrix, without duplicating the data storage. + * Each row in the matrix represents a point in the state space. + * + * Example of usage: + * \code + * Eigen::Matrix mat; + * // Fill out "mat"... + * + * typedef KDTreeEigenMatrixAdaptor< Eigen::Matrix > my_kd_tree_t; + * const int max_leaf = 10; + * my_kd_tree_t mat_index(dimdim, mat, max_leaf ); + * mat_index.index->buildIndex(); + * mat_index.index->... + * \endcode + * + * \tparam DIM If set to >0, it specifies a compile-time fixed dimensionality for the points in the data set, allowing more compiler optimizations. + * \tparam Distance The distance metric to use: nanoflann::metric_L1, nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc. + * \tparam IndexType The type for indices in the KD-tree index (typically, size_t of int) + */ + template + struct KDTreeEigenMatrixAdaptor + { + typedef KDTreeEigenMatrixAdaptor self_t; + typedef typename MatrixType::Scalar num_t; + typedef typename Distance::template traits::distance_t metric_t; + typedef KDTreeSingleIndexAdaptor< metric_t,self_t,DIM,IndexType> index_t; + + index_t* index; //! The kd-tree index for the user to call its methods as usual with any other FLANN index. + + /// Constructor: takes a const ref to the matrix object with the data points + KDTreeEigenMatrixAdaptor(const int dimensionality, const MatrixType &mat, const int leaf_max_size = 10) : m_data_matrix(mat) + { + const size_t dims = mat.cols(); + if (DIM>0 && static_cast(dims)!=DIM) + throw std::runtime_error("Data set dimensionality does not match the 'DIM' template argument"); + index = new index_t( dims, *this /* adaptor */, nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dims ) ); + index->buildIndex(); + } + + ~KDTreeEigenMatrixAdaptor() { + delete index; + } + + const MatrixType &m_data_matrix; + + /** Query for the \a num_closest closest points to a given point (entered as query_point[0:dim-1]). + * Note that this is a short-cut method for index->findNeighbors(). + * The user can also call index->... methods as desired. + * \note nChecks_IGNORED is ignored but kept for compatibility with the original FLANN interface. + */ + inline void query(const num_t *query_point, const size_t num_closest, IndexType *out_indices, num_t *out_distances_sq, const int nChecks_IGNORED = 10) const + { + nanoflann::KNNResultSet resultSet(num_closest); + resultSet.init(out_indices, out_distances_sq); + index->findNeighbors(resultSet, query_point, nanoflann::SearchParams()); + } + + /** @name Interface expected by KDTreeSingleIndexAdaptor + * @{ */ + + const self_t & derived() const { + return *this; + } + self_t & derived() { + return *this; + } + + // Must return the number of data points + inline size_t kdtree_get_point_count() const { + return m_data_matrix.rows(); + } + + // Returns the distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class: + inline num_t kdtree_distance(const num_t *p1, const size_t idx_p2,size_t size) const + { + num_t s=0; + for (size_t i=0; i + bool kdtree_get_bbox(BBOX &bb) const { + return false; + } + + /** @} */ + + }; // end of KDTreeEigenMatrixAdaptor + /** @} */ + +/** @} */ // end of grouping +} // end of NS + + +#endif /* NANOFLANN_HPP_ */ diff --git a/Tangential_complex/include/CGAL/Tangential_complex/utilities.h b/Tangential_complex/include/CGAL/Tangential_complex/utilities.h new file mode 100644 index 00000000000..3797268d6ae --- /dev/null +++ b/Tangential_complex/include/CGAL/Tangential_complex/utilities.h @@ -0,0 +1,76 @@ +// 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_TC_UTILITIES_H +#define CGAL_TC_UTILITIES_H + +#include +#include + +#include + +namespace CGAL { +namespace Tangential_complex_ { + + template + std::vector + compute_gram_schmidt_basis( + std::vector const& input_basis, + K const& kernel) + { + typedef typename K::FT FT; + typedef typename K::Vector_d Vector; + typedef std::vector Basis; + + const int D = Ambient_dimension::value; // CJTODO: use Point_dimension_d or similar + + // Kernel functors + typename K::Squared_length_d sqlen = kernel.squared_length_d_object(); + typename K::Scaled_vector_d scaled_vec = kernel.scaled_vector_d_object(); + typename K::Scalar_product_d inner_pdct = kernel.scalar_product_d_object(); + typename K::Difference_of_vectors_d diff_vec = kernel.difference_of_vectors_d_object(); + + Basis output_basis; + + typename Basis::const_iterator inb_it = input_basis.begin(); + typename Basis::const_iterator inb_it_end = input_basis.end(); + for (int i = 0 ; inb_it != inb_it_end ; ++inb_it, ++i) + { + Vector u = *inb_it; + + typename Basis::iterator outb_it = output_basis.begin(); + for (int j = 0 ; j < i ; ++j) + { + Vector const& ej = *outb_it; + Vector u_proj = scaled_vec(ej, inner_pdct(u, ej)); + u = diff_vec(u, u_proj); + } + + output_basis.push_back( + scaled_vec(u, FT(1)/CGAL::sqrt(sqlen(u)))); + } + + return output_basis; + } + +} // namespace Tangential_complex_ +} //namespace CGAL + +#endif // CGAL_TC_UTILITIES_H diff --git a/Tangential_complex/package_info/Tangential_complex/copyright b/Tangential_complex/package_info/Tangential_complex/copyright new file mode 100644 index 00000000000..8932b3233d2 --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/copyright @@ -0,0 +1,2 @@ +INRIA Sophia-Antipolis (France) + diff --git a/Tangential_complex/package_info/Tangential_complex/copyright.txt b/Tangential_complex/package_info/Tangential_complex/copyright.txt new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Tangential_complex/package_info/Tangential_complex/description.txt b/Tangential_complex/package_info/Tangential_complex/description.txt new file mode 100644 index 00000000000..b8969f7fc3e --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/description.txt @@ -0,0 +1,3 @@ +Tangential complex + +This CGAL component provides an implementation of the tangential complex. diff --git a/Tangential_complex/package_info/Tangential_complex/license.txt b/Tangential_complex/package_info/Tangential_complex/license.txt new file mode 100644 index 00000000000..8bb8efcb72b --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Tangential_complex/package_info/Tangential_complex/maintainer b/Tangential_complex/package_info/Tangential_complex/maintainer new file mode 100644 index 00000000000..769c1668e20 --- /dev/null +++ b/Tangential_complex/package_info/Tangential_complex/maintainer @@ -0,0 +1 @@ +ClĂ©ment Jamin \ No newline at end of file diff --git a/Tangential_complex/test/Tangential_complex/CMakeLists.txt b/Tangential_complex/test/Tangential_complex/CMakeLists.txt new file mode 100644 index 00000000000..c1b10d3b85c --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/CMakeLists.txt @@ -0,0 +1,77 @@ +# Created by the script cgal_create_cmake_script +# This is the CMake script for compiling a CGAL application. + + +project( Tangential_complex_test ) + +cmake_minimum_required(VERSION 2.6.2) +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6) + if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3) + cmake_policy(VERSION 2.8.4) + else() + cmake_policy(VERSION 2.6) + endif() +endif() + +# Creates a new CMake option, turned ON by default +option(ACTIVATE_MSVC_PRECOMPILED_HEADERS + "Activate precompiled headers in MSVC" + OFF) + +# Macro to add precompiled headers for MSVC +# This function does two things: +# 1. Enable precompiled headers on each file which is listed in "SourcesVar". +# 2. Add the content of "PrecompiledSource" (e.g. "StdAfx.cpp") to "SourcesVar". +MACRO(ADD_MSVC_PRECOMPILED_HEADER PrecompiledHeader PrecompiledSource SourcesVar) + IF(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS) + GET_FILENAME_COMPONENT(PrecompiledBasename ${PrecompiledHeader} NAME_WE) + SET(Sources ${${SourcesVar}}) + + SET_SOURCE_FILES_PROPERTIES(${PrecompiledSource} + PROPERTIES COMPILE_FLAGS "/Yc\"${PrecompiledHeader}\"") + SET_SOURCE_FILES_PROPERTIES(${Sources} + PROPERTIES COMPILE_FLAGS "/Yu\"${PrecompiledHeaders}\" /FI\"${PrecompiledHeader}\"") + # Add precompiled header to SourcesVar + LIST(APPEND ${SourcesVar} ${PrecompiledSource}) + ENDIF(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS) +ENDMACRO(ADD_MSVC_PRECOMPILED_HEADER) + +# The compiler might need more memory because of precompiled headers +if(MSVC AND ACTIVATE_MSVC_PRECOMPILED_HEADERS AND NOT(MSVC_VERSION LESS 1310)) + set(CGAL_C_FLAGS "${CGAL_C_FLAGS} /Zm1000") + set(CGAL_CXX_FLAGS "${CGAL_CXX_FLAGS} /Zm1000") +endif() + +find_package(CGAL QUIET COMPONENTS Core ) + +if ( CGAL_FOUND ) + include( ${CGAL_USE_FILE} ) + + find_package( TBB QUIET ) + + if( TBB_FOUND ) + include(${TBB_USE_FILE}) + list(APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES}) + endif() + + + include( CGAL_CreateSingleSourceCGALProgram ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + include_directories (BEFORE "include") + + set (SOURCE_FILES "test_tangential_complex.cpp") + ADD_MSVC_PRECOMPILED_HEADER("StdAfx.h" "StdAfx.cpp" SOURCE_FILES) + create_single_source_cgal_program( ${SOURCE_FILES} ) + + 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/Tangential_complex/test/Tangential_complex/StdAfx.cpp b/Tangential_complex/test/Tangential_complex/StdAfx.cpp new file mode 100644 index 00000000000..15668dcadef --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/StdAfx.cpp @@ -0,0 +1,2 @@ +// Build the precompiled headers. +#include "StdAfx.h" \ No newline at end of file diff --git a/Tangential_complex/test/Tangential_complex/StdAfx.h b/Tangential_complex/test/Tangential_complex/StdAfx.h new file mode 100644 index 00000000000..573345d0022 --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/StdAfx.h @@ -0,0 +1,14 @@ +#ifndef STDAFX_H +#define STDAFX_H + +// CGAL +#include +#include +#include +#include +#include +#include +#include +#include + +#endif //STDAFX_H \ No newline at end of file diff --git a/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp b/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp new file mode 100644 index 00000000000..77102284fbd --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/test_tangential_complex.cpp @@ -0,0 +1,124 @@ +//#undef CGAL_LINKED_WITH_TBB // CJTODO TEMP + +// Without TBB_USE_THREADING_TOOL Intel Inspector XE will report false positives in Intel TBB +// (http://software.intel.com/en-us/articles/compiler-settings-for-threading-error-analysis-in-intel-inspector-xe/) +#ifdef _DEBUG +# define TBB_USE_THREADING_TOOL +#endif + +#include +#include +#include +#include + +#include "test_utilities.h" + +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +# include +#endif + +#ifdef _DEBUG + const int NUM_POINTS = 50; +#else + const int NUM_POINTS = 10000; +#endif + +int main() +{ +#ifdef CGAL_LINKED_WITH_TBB +# ifdef _DEBUG + tbb::task_scheduler_init init(1); +# else + tbb::task_scheduler_init init(10); +# endif +#endif + + const int INTRINSIC_DIMENSION = 2; + const int AMBIENT_DIMENSION = 4; + + typedef CGAL::Epick_d > Kernel; + typedef Kernel::FT FT; + typedef Kernel::Point_d Point; + + int i = 0; + bool stop = false; + //for ( ; !stop ; ++i) + { + Kernel k; + Wall_clock_timer t; + CGAL::default_random = CGAL::Random(i); + std::cerr << "Random seed = " << i << std::endl; + +#ifdef CGAL_TC_PROFILING + Wall_clock_timer t_gen; +#endif + std::vector points = + //generate_points_on_circle_2(NUM_POINTS, 3.); + //generate_points_on_moment_curve(NUM_POINTS, AMBIENT_DIMENSION, 0., 1.); + //generate_points_on_plane(NUM_POINTS); + //generate_points_on_sphere_3(NUM_POINTS, 3.0); + //generate_points_on_sphere_d(NUM_POINTS, AMBIENT_DIMENSION, 3.0); + //generate_points_on_klein_bottle_3D(NUM_POINTS, 4., 3.); + generate_points_on_klein_bottle_4D(NUM_POINTS, 4., 3.); + //generate_points_on_klein_bottle_variant_5D(NUM_POINTS, 4., 3.); + +#ifdef CGAL_TC_PROFILING + std::cerr << "Point set generated in " << t_gen.elapsed() + << " seconds." << std::endl; +#endif + + points = sparsify_point_set( + k, points, FT(INPUT_SPARSITY)*FT(INPUT_SPARSITY)); + std::cerr << "Number of points after sparsification: " + << points.size() << std::endl; + + CGAL::Tangential_complex< + Kernel, + INTRINSIC_DIMENSION, + CGAL::Parallel_tag> tc(points.begin(), points.end(), k); + double init_time = t.elapsed(); t.reset(); + + tc.compute_tangential_complex(); + double computation_time = t.elapsed(); t.reset(); + + std::set > incorrect_simplices; + //stop = !tc.check_if_all_simplices_are_in_the_ambient_delaunay(&incorrect_simplices); + + t.reset(); + tc.fix_inconsistencies(); + double fix_time = t.elapsed(); t.reset(); + + double export_time = -1.; + if (INTRINSIC_DIMENSION <= 3) + { + t.reset(); + std::stringstream output_filename; + output_filename << "data/test_tc_" << INTRINSIC_DIMENSION + << "_in_R" << AMBIENT_DIMENSION << ".off"; + std::ofstream off_stream(output_filename.str().c_str()); + tc.export_to_off(off_stream, true, &incorrect_simplices, true); + double export_time = t.elapsed(); t.reset(); + } + /*else + tc.number_of_inconsistent_simplices();*/ + + + std::cerr << std::endl + << "================================================" << std::endl + << "Number of vertices: " << tc.number_of_vertices() << std::endl + << "Computation times (seconds): " << std::endl + << " * Tangential complex: " << init_time + computation_time + << std::endl + << " - Init + kd-tree = " << init_time << std::endl + << " - TC computation = " << computation_time << std::endl + << " * Fix inconsistencies: " << fix_time << std::endl + << " * Export to OFF: " << export_time << std::endl + << "================================================" << std::endl + << std::endl; + } + + return 0; +} diff --git a/Tangential_complex/test/Tangential_complex/test_utilities.h b/Tangential_complex/test/Tangential_complex/test_utilities.h new file mode 100644 index 00000000000..95a4a9d20e7 --- /dev/null +++ b/Tangential_complex/test/Tangential_complex/test_utilities.h @@ -0,0 +1,460 @@ +// 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 +// +//****************************************************************************** +// File Description : +// +//****************************************************************************** + +#ifndef CGAL_TC_TEST_TEST_UTILITIES_H +#define CGAL_TC_TEST_TEST_UTILITIES_H + +#include +#include +#include +#include +#include + +// Actually, this is very slow because the "m_points_ds->insert" +// cleans the tree, which is thus built at each query_ANN call +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER +template +class Point_sparsifier +{ +public: + typedef typename Kernel::FT FT; + typedef typename Point_container::value_type Point; + typedef typename CGAL::Tangential_complex_::Point_cloud_data_structure< + Kernel, Point_container> Points_ds; + typedef typename Points_ds::KNS_range KNS_range; + + // We can't instantiate m_points_ds right now since it requires that + // points is not empty (which be the case here) + Point_sparsifier(Point_container &points, + FT sparsity = FT(INPUT_SPARSITY*INPUT_SPARSITY)) + : m_points(points), m_sparsity(sparsity), m_points_ds(NULL) + {} + + bool try_to_insert_point(const Point &p) + { + if (m_points_ds == NULL) + { + m_points.push_back(p); + m_points_ds = new Points_ds(m_points); + m_points_ds->insert(0); + return true; + } + else + { + KNS_range kns_range = m_points_ds->query_ANN(p, 1, false); + if (kns_range.begin()->second >= m_sparsity) + { + m_points.push_back(p); + m_points_ds->insert(m_points.size() - 1); + return true; + } + } + + return false; + } + +private: + Point_container & m_points; + Points_ds * m_points_ds; + FT m_sparsity; +}; +#endif + +template +std::vector +sparsify_point_set( + const Kernel &k, Point_container const& input_pts, + typename Kernel::FT min_squared_dist) +{ + typedef typename Point_container::value_type Point; + typedef typename CGAL::Tangential_complex_::Point_cloud_data_structure< + Kernel, Point_container> Points_ds; + typedef typename Points_ds::INS_iterator INS_iterator; + typedef typename Points_ds::INS_range INS_range; + + typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); + +#ifdef CGAL_TC_PROFILING + Wall_clock_timer t; +#endif + + // Create the output container + std::vector output; + + Points_ds points_ds(input_pts); + + // Parse the following points, and add them if they are not too close to + // the other points + std::size_t pt_idx = 0; + for (typename Point_container::const_iterator it_pt = input_pts.begin() ; + it_pt != input_pts.end(); + ++it_pt, ++pt_idx) + { + INS_range ins_range = points_ds.query_incremental_ANN(*it_pt); + + // Drop it if there is another point that: + // - is closer that min_squared_dist + // - and has a higher index + for (INS_iterator nn_it = ins_range.begin() ; + nn_it != ins_range.end() ; + ++nn_it) + { + std::size_t neighbor_point_idx = nn_it->first; + typename Kernel::FT sq_dist = nn_it->second; + // The neighbor is further, we keep the point + if (sq_dist >= min_squared_dist) + { + output.push_back(*it_pt); + break; + } + // The neighbor is close and it has a higher index + else if (neighbor_point_idx > pt_idx) + break; // We drop the point + // Otherwise, we go the next closest point + } + } + +#ifdef CGAL_TC_PROFILING + std::cerr << "Point set sparsified in " << t.elapsed() + << " seconds." << std::endl; +#endif + + return output; +} + +template +bool load_points_from_file( + const std::string &filename, OutputIterator points, + std::size_t only_first_n_points = std::numeric_limits::max()) +{ + std::ifstream in(filename); + if (!in.is_open()) + { + std::cerr << "Could not open '" << filename << "'" << std::endl; + return false; + } + + Point p; + int dim_from_file; + in >> dim_from_file; + + std::size_t i = 0; + while(i < only_first_n_points && in >> p) + { + *points++ = p; + ++i; + } + + return true; +} + +template +std::vector generate_points_on_plane(std::size_t num_points) +{ + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + CGAL::Random rng; + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + FT x = rng.get_double(0, 5); + FT y = rng.get_double(0, 5); + Point p = Kernel().construct_point_d_object()(x, y, 0); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +template +std::vector generate_points_on_moment_curve( + std::size_t num_points, int dim, + typename Kernel::FT min_x , typename Kernel::FT max_x) +{ + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + CGAL::Random rng; + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + FT x = rng.get_double(min_x, max_x); + std::vector coords; + coords.reserve(dim); + for (int p = 1 ; p <= dim ; ++p) + coords.push_back(std::pow(CGAL::to_double(x), p)); + Point p = Kernel().construct_point_d_object()( + dim, coords.begin(), coords.end()); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +template +std::vector generate_points_on_circle_2( + std::size_t num_points, double radius) +{ + typedef typename Kernel::Point_d Point; + CGAL::Random_points_on_circle_2 generator(radius); + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + Point p = *generator++; +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +template +std::vector generate_points_on_sphere_3( + std::size_t num_points, double radius) +{ + typedef typename Kernel::Point_d Point; + CGAL::Random_points_on_sphere_3 generator(radius); + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + Point p = *generator++; +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +template +std::vector generate_points_on_sphere_d( + std::size_t num_points, int dim, double radius) +{ + typedef typename Kernel::Point_d Point; + CGAL::Random_points_on_sphere_d generator(dim, radius); + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + Point p = *generator++; +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +// a = big radius, b = small radius +template +std::vector generate_points_on_klein_bottle_3D( + std::size_t num_points, double a, double b, bool uniform = false) +{ + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + std::size_t num_cols = num_points/num_lines + 1; + + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + FT u, v; + if (uniform) + { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } + else + { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + double tmp = cos(u/2)*sin(v) - sin(u/2)*sin(2.*v); + Point p = Kernel().construct_point_d_object()( + (a + b*tmp)*cos(u), + (a + b*tmp)*sin(u), + b*(sin(u/2)*sin(v) + cos(u/2)*sin(2.*v))); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +// a = big radius, b = small radius +template +std::vector generate_points_on_klein_bottle_4D( + std::size_t num_points, double a, double b, double noise = 0., bool uniform = false) +{ + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + std::size_t num_cols = num_points/num_lines + 1; + + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + FT u, v; + if (uniform) + { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } + else + { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + Point p = Kernel().construct_point_d_object()( + Kernel().construct_point_d_object()( + (a + b*cos(v))*cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)), + (a + b*cos(v))*sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)), + b*sin(v)*cos(u/2) + (noise == 0. ? 0. : rng.get_double(0, noise)), + b*sin(v)*sin(u/2) + (noise == 0. ? 0. : rng.get_double(0, noise)))); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + + +// a = big radius, b = small radius +template +std::vector +generate_points_on_klein_bottle_variant_5D( + std::size_t num_points, double a, double b, bool uniform = false) +{ + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + std::size_t num_cols = num_points/num_lines + 1; + + std::vector points; + points.reserve(num_points); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + Point_sparsifier > sparsifier(points); +#endif + for (std::size_t i = 0 ; i < num_points ; ) + { + FT u, v; + if (uniform) + { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } + else + { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + FT x1 = (a + b*cos(v))*cos(u); + FT x2 = (a + b*cos(v))*sin(u); + FT x3 = b*sin(v)*cos(u/2); + FT x4 = b*sin(v)*sin(u/2); + FT x5 = x1 + x2 + x3 + x4; + + Point p = Kernel().construct_point_d_object()(x1, x2, x3, x4, x5); +#ifdef CGAL_TC_USE_SLOW_BUT_ACCURATE_SPARSIFIER + if (sparsifier.try_to_insert_point(p)) + ++i; +#else + points.push_back(p); + ++i; +#endif + } + return points; +} + +#endif // CGAL_MESH_3_TEST_TEST_UTILITIES_H diff --git a/Triangulation/include/CGAL/Delaunay_triangulation.h b/Triangulation/include/CGAL/Delaunay_triangulation.h index 75eabc243e0..a1c8015bb3e 100644 --- a/Triangulation/include/CGAL/Delaunay_triangulation.h +++ b/Triangulation/include/CGAL/Delaunay_triangulation.h @@ -75,6 +75,8 @@ public: // PUBLIC NESTED TYPES 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; @@ -118,7 +120,7 @@ public: using Base::vertices_begin; using Base::vertices_end; // using Base:: - + private: //*** Side_of_oriented_subsphere_d *** typedef typename Base::Flat_orientation_d Flat_orientation_d; @@ -155,12 +157,12 @@ public: } // With this constructor, - // the user can specify a Flat_orientation_d object to be used for - // orienting simplices of a specific dimension + // 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 Delaunay_triangulation( - int dim, + int dim, const std::pair &preset_flat_orientation, const Geom_traits &k = Geom_traits()) : Base(dim, preset_flat_orientation, k) @@ -175,8 +177,8 @@ public: Side_of_oriented_subsphere_d side_of_oriented_subsphere_predicate() const { return Side_of_oriented_subsphere_d ( - flat_orientation_, - geom_traits().construct_flat_orientation_d_object(), + flat_orientation_, + geom_traits().construct_flat_orientation_d_object(), geom_traits().in_flat_side_of_oriented_sphere_d_object() ); } @@ -288,7 +290,7 @@ public: Orientation o = ori_( boost::make_transform_iterator(s->vertices_begin(), spivi), - boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1, + boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1, spivi)); if( POSITIVE == o ) @@ -317,9 +319,9 @@ 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: @@ -336,7 +338,7 @@ private: Conflict_traversal_pred_in_fullspace; }; -// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = +// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = // FUNCTIONS THAT ARE MEMBER METHODS: // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS @@ -527,7 +529,7 @@ Delaunay_triangulation 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(); @@ -625,7 +627,7 @@ Delaunay_triangulation 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); @@ -855,27 +857,27 @@ 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 ) + for (Finite_full_cell_const_iterator cit = this->finite_full_cells_begin() ; + cit != this->finite_full_cells_end() ; ++cit ) { Full_cell_const_handle ch = cit.base(); - for(int i = 0; i < dim+1 ; ++i ) + for(int i = 0; i < dim+1 ; ++i ) { // If the i-th neighbor is not an infinite cell - Vertex_handle opposite_vh = + Vertex_handle opposite_vh = ch->neighbor(i)->vertex(ch->neighbor(i)->index(ch)); if (!is_infinite(opposite_vh)) { - Side_of_oriented_sphere_d side = + Side_of_oriented_sphere_d side = geom_traits().side_of_oriented_sphere_d_object(); - if (side(Point_const_iterator(ch->vertices_begin()), + if (side(Point_const_iterator(ch->vertices_begin()), Point_const_iterator(ch->vertices_end()), opposite_vh->point()) == ON_BOUNDED_SIDE) { diff --git a/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h b/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h index 590c2a6b0e6..1ca28265e21 100644 --- a/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h +++ b/Triangulation/include/CGAL/IO/Triangulation_off_ostream.h @@ -53,7 +53,7 @@ input_point(std::istream & is, const Traits &traits, P & p) { typedef typename Traits::FT FT; std::vector coords; - + std::string line; for(;;) { @@ -84,7 +84,7 @@ 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(;;) { @@ -124,8 +124,8 @@ operator>>(std::istream &is, typename Wrap::Weighted_point_d & wp) std::vector coords; while (line_sstr >> temp) coords.push_back(temp); - - std::vector::iterator last = coords.end() - 1; + + typename std::vector::iterator last = coords.end() - 1; P p = P(coords.begin(), last); wp = WP(p, *last); @@ -134,7 +134,7 @@ operator>>(std::istream &is, typename Wrap::Weighted_point_d & wp) template < class GT, class TDS > std::ostream & -export_triangulation_to_off(std::ostream & os, +export_triangulation_to_off(std::ostream & os, const Triangulation & tr, bool in_3D_export_surface_only = false) { @@ -147,7 +147,7 @@ export_triangulation_to_off(std::ostream & os, 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."; @@ -158,11 +158,11 @@ export_triangulation_to_off(std::ostream & 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(); + 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()); @@ -172,7 +172,7 @@ export_triangulation_to_off(std::ostream & os, index_of_vertex[it.base()] = i; } CGAL_assertion( i == n ); - + size_t number_of_triangles = 0; if (tr.maximal_dimension() == 2) { @@ -207,7 +207,7 @@ export_triangulation_to_off(std::ostream & os, output << index_of_vertex[*vit] << " "; } output << std::endl; - ++number_of_triangles; + ++number_of_triangles; } } } @@ -243,7 +243,7 @@ export_triangulation_to_off(std::ostream & os, } os << "OFF \n" - << n << " " + << n << " " << number_of_triangles << " 0\n" << output.str(); diff --git a/Triangulation/include/CGAL/Regular_triangulation.h b/Triangulation/include/CGAL/Regular_triangulation.h index 9d8853629c8..2755add7ab1 100644 --- a/Triangulation/include/CGAL/Regular_triangulation.h +++ b/Triangulation/include/CGAL/Regular_triangulation.h @@ -42,15 +42,15 @@ class Regular_triangulation { typedef typename RTTraits::Dimension Maximal_dimension_; typedef typename Default::Get< - TDS_, + TDS_, Triangulation_data_structure< Maximal_dimension_, Triangulation_vertex, - Triangulation_full_cell + 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; @@ -68,8 +68,9 @@ public: // PUBLIC NESTED TYPES typedef typename Base::Full_cell Full_cell; typedef typename Base::Facet Facet; typedef typename Base::Face Face; - + typedef Maximal_dimension_ Maximal_dimension; + typedef typename Base::Point Point; typedef typename RTTraits::Bare_point Bare_point; typedef typename RTTraits::Weighted_point Weighted_point; @@ -153,7 +154,7 @@ private: } }; public: - + // - - - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS Regular_triangulation(int dim, const Geom_traits &k = Geom_traits()) @@ -162,12 +163,12 @@ public: } // With this constructor, - // the user can specify a Flat_orientation_d object to be used for - // orienting simplices of a specific dimension + // 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, + int dim, const std::pair &preset_flat_orientation, const Geom_traits &k = Geom_traits()) : Base(dim, preset_flat_orientation, k) @@ -182,8 +183,8 @@ public: 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(), + flat_orientation_, + geom_traits().construct_flat_orientation_d_object(), geom_traits().in_flat_power_test_d_object() ); } @@ -248,13 +249,13 @@ public: return number_of_vertices() - n; } - Vertex_handle insert(const Weighted_point &, - const Locate_type, - const Face &, - const Facet &, + Vertex_handle insert(const Weighted_point &, + const Locate_type, + const Face &, + const Facet &, const Full_cell_handle); - Vertex_handle insert(const Weighted_point & p, + Vertex_handle insert(const Weighted_point & p, const Full_cell_handle start = Full_cell_handle()) { Locate_type lt; @@ -274,14 +275,14 @@ public: 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 &, + + 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()) @@ -294,7 +295,7 @@ public: } Vertex_handle insert_if_in_star( - const Weighted_point & p, const Vertex_handle star_center, + const Weighted_point & p, const Vertex_handle star_center, const Vertex_handle hint) { CGAL_assertion( Vertex_handle() != hint ); @@ -350,7 +351,7 @@ public: Orientation o = ori_( boost::make_transform_iterator(s->vertices_begin(), spivi), - boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1, + boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1, spivi)); if( POSITIVE == o ) @@ -379,16 +380,16 @@ public: 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, + does_cell_range_contain_vertex(InputIterator cz_begin, InputIterator cz_end, Vertex_handle vh) const { // Check all vertices @@ -407,7 +408,7 @@ private: template void - process_conflict_zone(InputIterator cz_begin, InputIterator cz_end, + process_conflict_zone(InputIterator cz_begin, InputIterator cz_end, OutputIterator vertices_out) const { // Get all vertices @@ -427,10 +428,10 @@ private: } } - + template void - process_cz_vertices_after_insertion(InputIterator vertices_begin, + process_cz_vertices_after_insertion(InputIterator vertices_begin, InputIterator vertices_end) { // Get all vertices @@ -458,14 +459,14 @@ private: 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 @@ -540,8 +541,8 @@ Regular_triangulation 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; @@ -660,7 +661,7 @@ Regular_triangulation 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(); @@ -685,8 +686,8 @@ Regular_triangulation 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(); + for(typename Dark_full_cells::iterator it = dark_incident_s.begin(); + it != dark_incident_s.end(); ++it) { (*it)->data().count_ += 1; @@ -760,7 +761,7 @@ Regular_triangulation 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); @@ -870,9 +871,9 @@ Regular_triangulation template< typename RTTraits, typename TDS > typename Regular_triangulation::Vertex_handle Regular_triangulation -::insert_if_in_star(const Weighted_point & p, +::insert_if_in_star(const Weighted_point & p, const Vertex_handle star_center, - const Locate_type lt, const Face & f, const Facet & ft, + const Locate_type lt, const Face & f, const Facet & ft, const Full_cell_handle s) { switch( lt ) @@ -914,7 +915,7 @@ Regular_triangulation template< typename RTTraits, typename TDS > typename Regular_triangulation::Vertex_handle Regular_triangulation -::insert_in_conflicting_cell(const Weighted_point & p, +::insert_in_conflicting_cell(const Weighted_point & p, const Full_cell_handle s, const Vertex_handle only_if_this_vertex_is_in_the_cz) { @@ -950,10 +951,10 @@ Regular_triangulation 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(), + && !does_cell_range_contain_vertex(cs.begin(), cs.end(), only_if_this_vertex_is_in_the_cz)) { return Vertex_handle(); @@ -962,7 +963,7 @@ Regular_triangulation // Otherwise, proceed with the insertion std::vector cz_vertices; cz_vertices.reserve(64); - process_conflict_zone(cs.begin(), cs.end(), + process_conflict_zone(cs.begin(), cs.end(), std::back_inserter(cz_vertices)); Vertex_handle ret = insert_in_hole(p, cs.begin(), cs.end(), ft); @@ -1043,8 +1044,8 @@ Regular_triangulation if( current_dimension() < maximal_dimension() ) { Conflict_pred_in_subspace c( - *this, p, - coaffine_orientation_predicate(), + *this, p, + coaffine_orientation_predicate(), power_test_in_flat_predicate()); return c(s); } @@ -1067,8 +1068,8 @@ Regular_triangulation if( current_dimension() < maximal_dimension() ) { Conflict_pred_in_subspace c( - *this, p, - coaffine_orientation_predicate(), + *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); @@ -1089,7 +1090,7 @@ template< typename RTTraits, typename TDS > bool Regular_triangulation ::is_valid(bool verbose, int level) const -{ +{ if (!Base::is_valid(verbose, level)) return false; @@ -1100,16 +1101,16 @@ Regular_triangulation cit != finite_full_cells_end() ; ++cit ) { Full_cell_const_handle ch = cit.base(); - for(int i = 0; i < dim+1 ; ++i ) + for(int i = 0; i < dim+1 ; ++i ) { // If the i-th neighbor is not an infinite cell - Vertex_handle opposite_vh = + Vertex_handle opposite_vh = ch->neighbor(i)->vertex(ch->neighbor(i)->index(ch)); if (!is_infinite(opposite_vh)) { - Power_test_d side = + Power_test_d side = geom_traits().power_test_d_object(); - if (side(Point_const_iterator(ch->vertices_begin()), + if (side(Point_const_iterator(ch->vertices_begin()), Point_const_iterator(ch->vertices_end()), opposite_vh->point()) == ON_POSITIVE_SIDE) { diff --git a/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h b/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h index 3f5a9630a85..c5eb0248859 100644 --- a/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h +++ b/Triangulation/include/CGAL/Regular_triangulation_euclidean_traits.h @@ -24,7 +24,6 @@ #include #include #include -#include #include diff --git a/Triangulation/include/CGAL/Triangulation_data_structure.h b/Triangulation/include/CGAL/Triangulation_data_structure.h index f6097292c13..94818cc4627 100644 --- a/Triangulation/include/CGAL/Triangulation_data_structure.h +++ b/Triangulation/include/CGAL/Triangulation_data_structure.h @@ -713,7 +713,9 @@ Triangulation_data_structure } } } - clear_visited_marks(start); + clear_visited_marks(start); // CJTODO: couldn't we use what is in "out" + // to make ot faster? (would require to + // replace the output iterator by a container) return ft; } @@ -999,7 +1001,7 @@ Triangulation_data_structure associate_vertex_with_full_cell(new_s, facet_index, v); set_neighbors(new_s, facet_index, - neighbor(old_s, facet_index), + outside_neighbor, mirror_index(old_s, facet_index)); // add the new full_cell to the list of new full_cells diff --git a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h index 427bb375010..08461059f56 100644 --- a/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h @@ -281,7 +281,7 @@ public: #endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO { #ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif size_type n = number_of_vertices(); @@ -622,7 +622,7 @@ public: size_type n = number_of_vertices(); #ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif // Parallel diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index cd47e98642b..8ffbe734f11 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -208,7 +208,7 @@ namespace CGAL { #endif #ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif size_type n = number_of_vertices(); @@ -618,7 +618,7 @@ namespace CGAL { size_type n = number_of_vertices(); #ifdef CGAL_TRIANGULATION_3_PROFILING - WallClockTimer t; + Wall_clock_timer t; #endif // Parallel