23 #ifndef TANGENTIAL_COMPLEX_H_ 24 #define TANGENTIAL_COMPLEX_H_ 26 #include <gudhi/Tangential_complex/config.h> 27 #include <gudhi/Tangential_complex/Simplicial_complex.h> 28 #include <gudhi/Tangential_complex/utilities.h> 29 #include <gudhi/Kd_tree_search.h> 30 #include <gudhi/console_color.h> 31 #include <gudhi/Clock.h> 32 #include <gudhi/Simplex_tree.h> 34 #include <CGAL/Default.h> 35 #include <CGAL/Dimension.h> 36 #include <CGAL/function_objects.h> 37 #include <CGAL/Epick_d.h> 38 #include <CGAL/Regular_triangulation_traits_adapter.h> 39 #include <CGAL/Regular_triangulation.h> 40 #include <CGAL/Delaunay_triangulation.h> 41 #include <CGAL/Combination_enumerator.h> 42 #include <CGAL/point_generators_d.h> 45 #include <Eigen/Eigen> 47 #include <boost/optional.hpp> 48 #include <boost/iterator/transform_iterator.hpp> 49 #include <boost/range/adaptor/transformed.hpp> 50 #include <boost/range/counting_range.hpp> 51 #include <boost/math/special_functions/factorials.hpp> 52 #include <boost/container/flat_set.hpp> 69 #include <tbb/parallel_for.h> 70 #include <tbb/combinable.h> 71 #include <tbb/mutex.h> 80 namespace tangential_complex {
82 using namespace internal;
86 Vertex_data(std::size_t data = (std::numeric_limits<std::size_t>::max)())
89 operator std::size_t() {
93 operator std::size_t()
const {
128 typename DimensionTag,
129 typename Concurrency_tag = CGAL::Parallel_tag,
130 typename Triangulation_ = CGAL::Default
134 typedef typename K::FT FT;
135 typedef typename K::Point_d Point;
136 typedef typename K::Weighted_point_d Weighted_point;
137 typedef typename K::Vector_d Vector;
139 typedef typename CGAL::Default::Get
142 CGAL::Regular_triangulation
144 CGAL::Epick_d<DimensionTag>,
145 CGAL::Triangulation_data_structure
147 typename CGAL::Epick_d<DimensionTag>::Dimension,
148 CGAL::Triangulation_vertex
150 CGAL::Regular_triangulation_traits_adapter< CGAL::Epick_d<DimensionTag> >, Vertex_data
152 CGAL::Triangulation_full_cell<CGAL::Regular_triangulation_traits_adapter< CGAL::Epick_d<DimensionTag> > >
155 >::type Triangulation;
156 typedef typename Triangulation::Geom_traits Tr_traits;
157 typedef typename Triangulation::Weighted_point Tr_point;
158 typedef typename Tr_traits::Base::Point_d Tr_bare_point;
159 typedef typename Triangulation::Vertex_handle Tr_vertex_handle;
160 typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle;
161 typedef typename Tr_traits::Vector_d Tr_vector;
163 #if defined(GUDHI_USE_TBB) 164 typedef tbb::mutex Mutex_for_perturb;
165 typedef Vector Translation_for_perturb;
166 typedef std::vector<Atomic_wrapper<FT> > Weights;
168 typedef Vector Translation_for_perturb;
169 typedef std::vector<FT> Weights;
171 typedef std::vector<Translation_for_perturb> Translations_for_perturb;
181 : m_tr(
new Triangulation(dim)) { }
184 destroy_triangulation();
187 Triangulation & construct_triangulation(
int dim) {
189 m_tr =
new Triangulation(dim);
193 void destroy_triangulation() {
198 Triangulation & tr() {
202 Triangulation
const& tr()
const {
206 Tr_vertex_handle
const& center_vertex()
const {
207 return m_center_vertex;
210 Tr_vertex_handle & center_vertex() {
211 return m_center_vertex;
216 Tr_vertex_handle m_center_vertex;
220 typedef Basis<K> Tangent_space_basis;
221 typedef Basis<K> Orthogonal_space_basis;
222 typedef std::vector<Tangent_space_basis> TS_container;
223 typedef std::vector<Orthogonal_space_basis> OS_container;
225 typedef std::vector<Point> Points;
227 typedef boost::container::flat_set<std::size_t> Simplex;
228 typedef std::set<Simplex> Simplex_set;
235 typedef std::vector<Tr_and_VH> Tr_container;
236 typedef std::vector<Vector> Vectors;
240 typedef boost::container::flat_set<std::size_t> Incident_simplex;
241 typedef std::vector<Incident_simplex> Star;
242 typedef std::vector<Star> Stars_container;
246 static const Tr_point &vertex_handle_to_point(Tr_vertex_handle vh) {
250 template <
typename P,
typename VH>
251 static const P &vertex_handle_to_point(VH vh) {
256 typedef internal::Simplicial_complex Simplicial_complex;
267 template <
typename Po
int_range>
269 int intrinsic_dimension,
270 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
271 InputIterator first_for_tse, InputIterator last_for_tse,
276 m_intrinsic_dim(intrinsic_dimension),
277 m_ambient_dim(points.empty() ? 0 : k.point_dimension_d_object()(*points.begin())),
278 m_points(points.begin(), points.end()),
279 m_weights(m_points.size(), FT(0))
280 #if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION)
281 , m_p_perturb_mutexes(NULL)
283 , m_points_ds(m_points)
284 , m_last_max_perturb(0.)
285 , m_are_tangent_spaces_computed(m_points.size(), false)
286 , m_tangent_spaces(m_points.size(), Tangent_space_basis())
287 #ifdef GUDHI_TC_EXPORT_NORMALS
288 , m_orth_spaces(m_points.size(), Orthogonal_space_basis())
290 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
291 , m_points_for_tse(first_for_tse, last_for_tse)
292 , m_points_ds_for_tse(m_points_for_tse)
298 #if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION) 299 delete [] m_p_perturb_mutexes;
305 return m_intrinsic_dim;
310 return m_ambient_dim;
313 Points
const& points()
const {
323 return m_points[vertex];
332 return compute_perturbed_point(vertex);
338 return m_points.size();
341 void set_weights(
const Weights& weights) {
345 void set_tangent_planes(
const TS_container& tangent_spaces
346 #ifdef GUDHI_TC_EXPORT_NORMALS
347 ,
const OS_container& orthogonal_spaces
350 #ifdef GUDHI_TC_EXPORT_NORMALS 352 m_points.size() == tangent_spaces.size()
353 && m_points.size() == orthogonal_spaces.size(),
354 std::logic_error(
"Wrong sizes"));
357 m_points.size() == tangent_spaces.size(),
358 std::logic_error(
"Wrong sizes"));
360 m_tangent_spaces = tangent_spaces;
361 #ifdef GUDHI_TC_EXPORT_NORMALS 362 m_orth_spaces = orthogonal_spaces;
364 for (std::size_t i = 0; i < m_points.size(); ++i)
365 m_are_tangent_spaces_computed[i] =
true;
370 #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS 371 std::cerr << red <<
"WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. " 372 <<
"Computation might be slower than usual.\n" << white;
375 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB) 382 m_triangulations.resize(m_points.size());
383 m_stars.resize(m_points.size());
384 m_squared_star_spheres_radii_incl_margin.resize(m_points.size(), FT(-1));
385 #ifdef GUDHI_TC_PERTURB_POSITION 386 if (m_points.empty())
387 m_translations.clear();
389 m_translations.resize(m_points.size(),
390 m_k.construct_vector_d_object()(m_ambient_dim));
391 #if defined(GUDHI_USE_TBB) 392 delete [] m_p_perturb_mutexes;
393 m_p_perturb_mutexes =
new Mutex_for_perturb[m_points.size()];
399 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
400 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()),
401 Compute_tangent_triangulation(*
this));
403 #endif // GUDHI_USE_TBB 405 for (std::size_t i = 0; i < m_points.size(); ++i)
406 compute_tangent_triangulation(i);
409 #endif // GUDHI_USE_TBB 411 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB) 413 std::cerr <<
"Tangential complex computed in " << t.num_seconds()
421 bool success =
false;
423 unsigned int num_steps = 0;
425 std::size_t initial_num_inconsistent_stars = 0;
427 std::size_t best_num_inconsistent_stars = 0;
429 std::size_t final_num_inconsistent_stars = 0;
440 if (time_limit == 0.)
445 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 446 std::tuple<std::size_t, std::size_t, std::size_t> stats_before =
447 number_of_inconsistent_simplices(
false);
449 if (std::get<1>(stats_before) == 0) {
451 std::cerr <<
"Nothing to fix.\n";
456 #endif // GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 458 m_last_max_perturb = max_perturb;
464 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 466 <<
"\nBefore fix step:\n" 467 <<
" * Total number of simplices in stars (incl. duplicates): " 468 << std::get<0>(stats_before) <<
"\n" 469 <<
" * Num inconsistent simplices in stars (incl. duplicates): " 470 << red << std::get<1>(stats_before) << white <<
" (" 471 << 100. * std::get<1>(stats_before) / std::get<0>(stats_before) <<
"%)\n" 472 <<
" * Number of stars containing inconsistent simplices: " 473 << red << std::get<2>(stats_before) << white <<
" (" 474 << 100. * std::get<2>(stats_before) / m_points.size() <<
"%)\n";
477 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 479 <<
"\nAttempt to fix inconsistencies using perturbations - step #" 480 << info.
num_steps + 1 <<
"... " << white;
483 std::size_t num_inconsistent_stars = 0;
484 std::vector<std::size_t> updated_points;
486 #ifdef GUDHI_TC_PROFILING 487 Gudhi::Clock t_fix_step;
491 #if defined(GUDHI_USE_TBB) 492 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
493 tbb::combinable<std::size_t> num_inconsistencies;
494 tbb::combinable<std::vector<std::size_t> > tls_updated_points;
496 tbb::blocked_range<size_t>(0, m_triangulations.size()),
497 Try_to_solve_inconsistencies_in_a_local_triangulation(*
this, max_perturb,
499 tls_updated_points));
500 num_inconsistent_stars =
501 num_inconsistencies.combine(std::plus<std::size_t>());
502 updated_points = tls_updated_points.combine(
503 [](std::vector<std::size_t>
const& x,
504 std::vector<std::size_t>
const& y) {
505 std::vector<std::size_t> res;
506 res.reserve(x.size() + y.size());
507 res.insert(res.end(), x.begin(), x.end());
508 res.insert(res.end(), y.begin(), y.end());
512 #endif // GUDHI_USE_TBB 514 for (std::size_t i = 0; i < m_triangulations.size(); ++i) {
515 num_inconsistent_stars +=
516 try_to_solve_inconsistencies_in_a_local_triangulation(i, max_perturb,
517 std::back_inserter(updated_points));
519 #if defined(GUDHI_USE_TBB) 521 #endif // GUDHI_USE_TBB 523 #ifdef GUDHI_TC_PROFILING 527 #if defined(GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES) || defined(DEBUG_TRACES) 529 <<
"\nEncountered during fix:\n" 530 <<
" * Num stars containing inconsistent simplices: " 531 << red << num_inconsistent_stars << white
532 <<
" (" << 100. * num_inconsistent_stars / m_points.size() <<
"%)\n";
535 #ifdef GUDHI_TC_PROFILING 536 std::cerr << yellow <<
"done in " << t_fix_step.num_seconds()
537 <<
" seconds.\n" << white;
538 #elif defined(DEBUG_TRACES) 539 std::cerr << yellow <<
"done.\n" << white;
542 if (num_inconsistent_stars > 0)
543 refresh_tangential_complex(updated_points);
545 #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS 547 std::size_t num_inc_1 =
548 std::get<1>(number_of_inconsistent_simplices(
false));
549 refresh_tangential_complex();
550 std::size_t num_inc_2 =
551 std::get<1>(number_of_inconsistent_simplices(
false));
552 if (num_inc_1 != num_inc_2)
553 std::cerr << red <<
"REFRESHMENT CHECK: FAILED. (" 554 << num_inc_1 <<
" vs " << num_inc_2 <<
")\n" << white;
556 std::cerr << green <<
"REFRESHMENT CHECK: PASSED.\n" << white;
559 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 560 std::tuple<std::size_t, std::size_t, std::size_t> stats_after =
561 number_of_inconsistent_simplices(
false);
565 <<
" * Total number of simplices in stars (incl. duplicates): " 566 << std::get<0>(stats_after) <<
"\n" 567 <<
" * Num inconsistent simplices in stars (incl. duplicates): " 568 << red << std::get<1>(stats_after) << white <<
" (" 569 << 100. * std::get<1>(stats_after) / std::get<0>(stats_after) <<
"%)\n" 570 <<
" * Number of stars containing inconsistent simplices: " 571 << red << std::get<2>(stats_after) << white <<
" (" 572 << 100. * std::get<2>(stats_after) / m_points.size() <<
"%)\n";
574 stats_before = stats_after;
585 done = (num_inconsistent_stars == 0);
588 if (time_limit > 0. && t.num_seconds() > time_limit) {
590 std::cerr << red <<
"Time limit reached.\n" << white;
599 std::cerr << green <<
"Fixed!\n" << white;
608 std::size_t num_simplices = 0;
610 std::size_t num_inconsistent_simplices = 0;
612 std::size_t num_inconsistent_stars = 0;
629 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
630 bool is_star_inconsistent =
false;
633 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
634 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
635 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
637 if (is_infinite(*it_inc_simplex))
640 Simplex c = *it_inc_simplex;
643 if (!is_simplex_consistent(c)) {
645 is_star_inconsistent =
true;
655 <<
"\n==========================================================\n" 656 <<
"Inconsistencies:\n" 657 <<
" * Total number of simplices in stars (incl. duplicates): " 659 <<
" * Number of inconsistent simplices in stars (incl. duplicates): " 662 <<
" * Number of stars containing inconsistent simplices: " 665 <<
"==========================================================\n";
681 template <
typename Simplex_tree_>
683 ,
bool export_inconsistent_simplices =
true 685 ,
bool export_infinite_simplices =
false 686 , Simplex_set *p_inconsistent_simplices = NULL
689 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 691 <<
"\nExporting the TC as a Simplex_tree... " << white;
693 #ifdef GUDHI_TC_PROFILING 700 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
702 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
703 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
704 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
705 Simplex c = *it_inc_simplex;
708 if (!export_infinite_simplices && is_infinite(c))
711 if (!export_inconsistent_simplices && !is_simplex_consistent(c))
714 if (static_cast<int> (c.size()) > max_dim)
715 max_dim =
static_cast<int> (c.size());
720 bool inserted = tree.insert_simplex_and_subfaces(c).second;
723 if (p_inconsistent_simplices && inserted && !is_simplex_consistent(c)) {
724 p_inconsistent_simplices->insert(c);
729 #ifdef GUDHI_TC_PROFILING 731 std::cerr << yellow <<
"done in " << t.num_seconds()
732 <<
" seconds.\n" << white;
733 #elif defined(DEBUG_TRACES) 734 std::cerr << yellow <<
"done.\n" << white;
750 int create_complex(Simplicial_complex &complex,
751 bool export_inconsistent_simplices =
true,
752 bool export_infinite_simplices =
false,
753 int check_lower_and_higher_dim_simplices = 2,
754 Simplex_set *p_inconsistent_simplices = NULL)
const {
755 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 757 <<
"\nExporting the TC as a Simplicial_complex... " << white;
759 #ifdef GUDHI_TC_PROFILING 767 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
769 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
770 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
771 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
772 Simplex c = *it_inc_simplex;
775 if (!export_infinite_simplices && is_infinite(c))
778 if (!export_inconsistent_simplices && !is_simplex_consistent(c))
782 if (check_lower_and_higher_dim_simplices == 2
784 && static_cast<int> (c.size()) != max_dim) {
787 "Info: check_lower_and_higher_dim_simplices ACTIVATED. " 788 "Export might be take some time...\n" << white;
789 check_lower_and_higher_dim_simplices = 1;
792 if (static_cast<int> (c.size()) > max_dim)
793 max_dim =
static_cast<int> (c.size());
799 complex.add_simplex(c, check_lower_and_higher_dim_simplices == 1);
802 if (p_inconsistent_simplices && added && !is_simplex_consistent(c)) {
803 p_inconsistent_simplices->insert(c);
808 #ifdef GUDHI_TC_PROFILING 810 std::cerr << yellow <<
"done in " << t.num_seconds()
811 <<
" seconds.\n" << white;
812 #elif defined(DEBUG_TRACES) 813 std::cerr << yellow <<
"done.\n" << white;
819 template<
typename ProjectionFunctor = CGAL::Identity<Po
int> >
820 std::ostream &export_to_off(
821 const Simplicial_complex &complex, std::ostream & os,
822 Simplex_set
const *p_simpl_to_color_in_red = NULL,
823 Simplex_set
const *p_simpl_to_color_in_green = NULL,
824 Simplex_set
const *p_simpl_to_color_in_blue = NULL,
825 ProjectionFunctor
const& point_projection = ProjectionFunctor())
827 return export_to_off(
828 os,
false, p_simpl_to_color_in_red, p_simpl_to_color_in_green,
829 p_simpl_to_color_in_blue, &complex, point_projection);
832 template<
typename ProjectionFunctor = CGAL::Identity<Po
int> >
833 std::ostream &export_to_off(
834 std::ostream & os,
bool color_inconsistencies =
false,
835 Simplex_set
const *p_simpl_to_color_in_red = NULL,
836 Simplex_set
const *p_simpl_to_color_in_green = NULL,
837 Simplex_set
const *p_simpl_to_color_in_blue = NULL,
838 const Simplicial_complex *p_complex = NULL,
839 ProjectionFunctor
const& point_projection = ProjectionFunctor())
const {
840 if (m_points.empty())
843 if (m_ambient_dim < 2) {
844 std::cerr <<
"Error: export_to_off => ambient dimension should be >= 2.\n";
845 os <<
"Error: export_to_off => ambient dimension should be >= 2.\n";
848 if (m_ambient_dim > 3) {
849 std::cerr <<
"Warning: export_to_off => ambient dimension should be " 850 "<= 3. Only the first 3 coordinates will be exported.\n";
853 if (m_intrinsic_dim < 1 || m_intrinsic_dim > 3) {
854 std::cerr <<
"Error: export_to_off => intrinsic dimension should be " 855 "between 1 and 3.\n";
856 os <<
"Error: export_to_off => intrinsic dimension should be " 857 "between 1 and 3.\n";
861 std::stringstream output;
862 std::size_t num_simplices, num_vertices;
863 export_vertices_to_off(output, num_vertices,
false, point_projection);
865 export_simplices_to_off(
866 *p_complex, output, num_simplices, p_simpl_to_color_in_red,
867 p_simpl_to_color_in_green, p_simpl_to_color_in_blue);
869 export_simplices_to_off(
870 output, num_simplices, color_inconsistencies, p_simpl_to_color_in_red,
871 p_simpl_to_color_in_green, p_simpl_to_color_in_blue);
874 #ifdef GUDHI_TC_EXPORT_NORMALS 879 << num_vertices <<
" " 880 << num_simplices <<
" " 888 void refresh_tangential_complex() {
889 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 890 std::cerr << yellow <<
"\nRefreshing TC... " << white;
893 #ifdef GUDHI_TC_PROFILING 898 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
899 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()),
900 Compute_tangent_triangulation(*
this));
902 #endif // GUDHI_USE_TBB 904 for (std::size_t i = 0; i < m_points.size(); ++i)
905 compute_tangent_triangulation(i);
908 #endif // GUDHI_USE_TBB 910 #ifdef GUDHI_TC_PROFILING 912 std::cerr << yellow <<
"done in " << t.num_seconds()
913 <<
" seconds.\n" << white;
914 #elif defined(DEBUG_TRACES) 915 std::cerr << yellow <<
"done.\n" << white;
920 template <
typename Po
int_indices_range>
921 void refresh_tangential_complex(
922 Point_indices_range
const& perturbed_points_indices) {
923 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 924 std::cerr << yellow <<
"\nRefreshing TC... " << white;
927 #ifdef GUDHI_TC_PROFILING 932 Points_ds updated_pts_ds(m_points, perturbed_points_indices);
936 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
937 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()),
938 Refresh_tangent_triangulation(*
this, updated_pts_ds));
940 #endif // GUDHI_USE_TBB 942 for (std::size_t i = 0; i < m_points.size(); ++i)
943 refresh_tangent_triangulation(i, updated_pts_ds);
946 #endif // GUDHI_USE_TBB 948 #ifdef GUDHI_TC_PROFILING 950 std::cerr << yellow <<
"done in " << t.num_seconds()
951 <<
" seconds.\n" << white;
952 #elif defined(DEBUG_TRACES) 953 std::cerr << yellow <<
"done.\n" << white;
957 void export_inconsistent_stars_to_OFF_files(std::string
const& filename_base)
const {
959 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
961 Simplicial_complex sc;
963 bool is_inconsistent =
false;
964 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
965 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
966 for (; it_inc_simplex != it_inc_simplex_end;
969 if (is_infinite(*it_inc_simplex))
972 Simplex c = *it_inc_simplex;
978 if (!is_inconsistent && !is_simplex_consistent(c))
979 is_inconsistent =
true;
982 if (is_inconsistent) {
984 std::stringstream output_filename;
985 output_filename << filename_base <<
"_" << idx <<
".off";
986 std::ofstream off_stream(output_filename.str().c_str());
987 export_to_off(sc, off_stream);
992 class Compare_distance_to_ref_point {
994 Compare_distance_to_ref_point(Point
const& ref, K
const& k)
995 : m_ref(ref), m_k(k) { }
997 bool operator()(Point
const& p1, Point
const& p2) {
998 typename K::Squared_distance_d sqdist =
999 m_k.squared_distance_d_object();
1000 return sqdist(p1, m_ref) < sqdist(p2, m_ref);
1008 #ifdef GUDHI_USE_TBB 1010 class Compute_tangent_triangulation {
1019 Compute_tangent_triangulation(
const Compute_tangent_triangulation &ctt)
1020 : m_tc(ctt.m_tc) { }
1023 void operator()(
const tbb::blocked_range<size_t>& r)
const {
1024 for (
size_t i = r.begin(); i != r.end(); ++i)
1025 m_tc.compute_tangent_triangulation(i);
1030 class Refresh_tangent_triangulation {
1032 Points_ds
const& m_updated_pts_ds;
1036 Refresh_tangent_triangulation(
Tangential_complex &tc, Points_ds
const& updated_pts_ds)
1037 : m_tc(tc), m_updated_pts_ds(updated_pts_ds) { }
1040 Refresh_tangent_triangulation(
const Refresh_tangent_triangulation &ctt)
1041 : m_tc(ctt.m_tc), m_updated_pts_ds(ctt.m_updated_pts_ds) { }
1044 void operator()(
const tbb::blocked_range<size_t>& r)
const {
1045 for (
size_t i = r.begin(); i != r.end(); ++i)
1046 m_tc.refresh_tangent_triangulation(i, m_updated_pts_ds);
1049 #endif // GUDHI_USE_TBB 1051 bool is_infinite(Simplex
const& s)
const {
1052 return *s.rbegin() == (std::numeric_limits<std::size_t>::max)();
1058 Tr_vertex_handle compute_star(std::size_t i,
const Point ¢er_pt,
const Tangent_space_basis &tsb,
1059 Triangulation &triangulation,
bool verbose =
false) {
1060 int tangent_space_dim = tsb.dimension();
1061 const Tr_traits &local_tr_traits = triangulation.geom_traits();
1064 typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
1067 typename Tr_traits::Compute_weight_d point_weight = local_tr_traits.compute_weight_d_object();
1068 typename Tr_traits::Power_center_d power_center = local_tr_traits.power_center_d_object();
1077 if (i == tsb.origin()) {
1079 proj_wp = local_tr_traits.construct_weighted_point_d_object()(local_tr_traits.construct_point_d_object()(tangent_space_dim, CGAL::ORIGIN),
1082 const Weighted_point& wp = compute_perturbed_weighted_point(i);
1083 proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
1086 Tr_vertex_handle center_vertex = triangulation.insert(proj_wp);
1087 center_vertex->data() = i;
1089 std::cerr <<
"* Inserted point #" << i <<
"\n";
1091 #ifdef GUDHI_TC_VERY_VERBOSE 1092 std::size_t num_attempts_to_insert_points = 1;
1093 std::size_t num_inserted_points = 1;
1097 INS_range ins_range = m_points_ds.incremental_nearest_neighbors(center_pt);
1103 boost::optional<FT> squared_star_sphere_radius_plus_margin;
1106 for (
auto nn_it = ins_range.begin();
1107 nn_it != ins_range.end();
1109 std::size_t neighbor_point_idx = nn_it->first;
1112 if (neighbor_point_idx != i) {
1117 compute_perturbed_weighted_point(neighbor_point_idx, neighbor_pt, neighbor_weight);
1119 if (squared_star_sphere_radius_plus_margin &&
1120 k_sqdist(center_pt, neighbor_pt) > *squared_star_sphere_radius_plus_margin)
1123 Tr_point proj_pt = project_point_and_compute_weight(neighbor_pt, neighbor_weight, tsb,
1126 #ifdef GUDHI_TC_VERY_VERBOSE 1127 ++num_attempts_to_insert_points;
1131 Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex);
1133 if (vh != Tr_vertex_handle() && vh->data() == (std::numeric_limits<std::size_t>::max)()) {
1134 #ifdef GUDHI_TC_VERY_VERBOSE 1135 ++num_inserted_points;
1138 std::cerr <<
"* Inserted point #" << neighbor_point_idx <<
"\n";
1140 vh->data() = neighbor_point_idx;
1143 if (triangulation.current_dimension() >= tangent_space_dim) {
1144 squared_star_sphere_radius_plus_margin = boost::none;
1146 std::vector<Tr_full_cell_handle> incident_cells;
1147 triangulation.incident_full_cells(
1149 std::back_inserter(incident_cells));
1150 for (
typename std::vector<Tr_full_cell_handle>::iterator cit =
1151 incident_cells.begin(); cit != incident_cells.end(); ++cit) {
1152 Tr_full_cell_handle cell = *cit;
1153 if (triangulation.is_infinite(cell)) {
1154 squared_star_sphere_radius_plus_margin = boost::none;
1159 Tr_point c = power_center(boost::make_transform_iterator(cell->vertices_begin(),
1160 vertex_handle_to_point<Tr_point,
1162 boost::make_transform_iterator(cell->vertices_end(),
1163 vertex_handle_to_point<Tr_point,
1164 Tr_vertex_handle>));
1166 FT sq_power_sphere_diam = 4 * point_weight(c);
1168 if (!squared_star_sphere_radius_plus_margin ||
1169 sq_power_sphere_diam > *squared_star_sphere_radius_plus_margin) {
1170 squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
1177 if (squared_star_sphere_radius_plus_margin) {
1179 squared_star_sphere_radius_plus_margin = CGAL::square(std::sqrt(*squared_star_sphere_radius_plus_margin)
1180 + 2 * m_last_max_perturb);
1183 m_squared_star_spheres_radii_incl_margin[i] =
1184 *squared_star_sphere_radius_plus_margin;
1186 m_squared_star_spheres_radii_incl_margin[i] = FT(-1);
1193 return center_vertex;
1196 void refresh_tangent_triangulation(std::size_t i, Points_ds
const& updated_pts_ds,
bool verbose =
false) {
1198 std::cerr <<
"** Refreshing tangent tri #" << i <<
" **\n";
1200 if (m_squared_star_spheres_radii_incl_margin[i] == FT(-1))
1201 return compute_tangent_triangulation(i, verbose);
1203 Point center_point = compute_perturbed_point(i);
1205 std::size_t closest_pt_index =
1208 typename K::Construct_weighted_point_d k_constr_wp =
1209 m_k.construct_weighted_point_d_object();
1210 typename K::Power_distance_d k_power_dist = m_k.power_distance_d_object();
1213 Weighted_point star_sphere = k_constr_wp(compute_perturbed_point(i),
1214 m_squared_star_spheres_radii_incl_margin[i]);
1215 Weighted_point closest_updated_point =
1216 compute_perturbed_weighted_point(closest_pt_index);
1219 if (k_power_dist(star_sphere, closest_updated_point) <= FT(0))
1220 compute_tangent_triangulation(i, verbose);
1223 void compute_tangent_triangulation(std::size_t i,
bool verbose =
false) {
1225 std::cerr <<
"** Computing tangent tri #" << i <<
" **\n";
1230 const Point center_pt = compute_perturbed_point(i);
1231 Tangent_space_basis &tsb = m_tangent_spaces[i];
1234 if (!m_are_tangent_spaces_computed[i]) {
1235 #ifdef GUDHI_TC_EXPORT_NORMALS 1236 tsb = compute_tangent_space(center_pt, i,
true , &m_orth_spaces[i]);
1238 tsb = compute_tangent_space(center_pt, i);
1242 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1245 int tangent_space_dim = tangent_basis_dim(i);
1246 Triangulation &local_tr =
1247 m_triangulations[i].construct_triangulation(tangent_space_dim);
1249 m_triangulations[i].center_vertex() =
1250 compute_star(i, center_pt, tsb, local_tr, verbose);
1252 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1254 std::cerr <<
" - triangulation construction: " << t.num_seconds() <<
" s.\n";
1258 #ifdef GUDHI_TC_VERY_VERBOSE 1259 std::cerr <<
"Inserted " << num_inserted_points <<
" points / " 1260 << num_attempts_to_insert_points <<
" attemps to compute the star\n";
1265 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1267 std::cerr <<
" - update_star: " << t.num_seconds() <<
" s.\n";
1273 void update_star(std::size_t i) {
1274 Star &star = m_stars[i];
1276 Triangulation &local_tr = m_triangulations[i].tr();
1277 Tr_vertex_handle center_vertex = m_triangulations[i].center_vertex();
1278 int cur_dim_plus_1 = local_tr.current_dimension() + 1;
1280 std::vector<Tr_full_cell_handle> incident_cells;
1281 local_tr.incident_full_cells(
1282 center_vertex, std::back_inserter(incident_cells));
1284 typename std::vector<Tr_full_cell_handle>::const_iterator it_c = incident_cells.begin();
1285 typename std::vector<Tr_full_cell_handle>::const_iterator it_c_end = incident_cells.end();
1287 for (; it_c != it_c_end; ++it_c) {
1289 Incident_simplex incident_simplex;
1290 for (
int j = 0; j < cur_dim_plus_1; ++j) {
1291 std::size_t index = (*it_c)->vertex(j)->data();
1293 incident_simplex.insert(index);
1295 GUDHI_CHECK(incident_simplex.size() == cur_dim_plus_1 - 1,
1296 std::logic_error(
"update_star: wrong size of incident simplex"));
1297 star.push_back(incident_simplex);
1303 Tangent_space_basis compute_tangent_space(
const Point &p
1304 ,
const std::size_t i
1305 ,
bool normalize_basis =
true 1306 , Orthogonal_space_basis *p_orth_space_basis = NULL
1308 unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1309 static_cast<unsigned int> (m_points.size()));
1312 typename K::Construct_vector_d constr_vec =
1313 m_k.construct_vector_d_object();
1314 typename K::Compute_coordinate_d coord =
1315 m_k.compute_coordinate_d_object();
1317 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 1318 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca,
false);
1319 const Points &points_for_pca = m_points_for_tse;
1321 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca,
false);
1322 const Points &points_for_pca = m_points;
1326 Eigen::MatrixXd mat_points(num_pts_for_pca, m_ambient_dim);
1327 auto nn_it = kns_range.begin();
1328 for (
unsigned int j = 0;
1329 j < num_pts_for_pca && nn_it != kns_range.end();
1331 for (
int i = 0; i < m_ambient_dim; ++i) {
1332 mat_points(j, i) = CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1335 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1336 Eigen::MatrixXd cov = centered.adjoint() * centered;
1337 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1339 Tangent_space_basis tsb(i);
1343 for (
int j = m_ambient_dim - 1;
1344 j >= m_ambient_dim - m_intrinsic_dim;
1346 if (normalize_basis) {
1347 Vector v = constr_vec(m_ambient_dim,
1348 eig.eigenvectors().col(j).data(),
1349 eig.eigenvectors().col(j).data() + m_ambient_dim);
1350 tsb.push_back(normalize_vector(v, m_k));
1352 tsb.push_back(constr_vec(
1354 eig.eigenvectors().col(j).data(),
1355 eig.eigenvectors().col(j).data() + m_ambient_dim));
1359 if (p_orth_space_basis) {
1360 p_orth_space_basis->set_origin(i);
1361 for (
int j = m_ambient_dim - m_intrinsic_dim - 1;
1364 if (normalize_basis) {
1365 Vector v = constr_vec(m_ambient_dim,
1366 eig.eigenvectors().col(j).data(),
1367 eig.eigenvectors().col(j).data() + m_ambient_dim);
1368 p_orth_space_basis->push_back(normalize_vector(v, m_k));
1370 p_orth_space_basis->push_back(constr_vec(
1372 eig.eigenvectors().col(j).data(),
1373 eig.eigenvectors().col(j).data() + m_ambient_dim));
1378 m_are_tangent_spaces_computed[i] =
true;
1388 Tangent_space_basis compute_tangent_space(
const Simplex &s,
bool normalize_basis =
true) {
1389 unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1390 static_cast<unsigned int> (m_points.size()));
1393 typename K::Construct_vector_d constr_vec =
1394 m_k.construct_vector_d_object();
1395 typename K::Compute_coordinate_d coord =
1396 m_k.compute_coordinate_d_object();
1397 typename K::Squared_length_d sqlen =
1398 m_k.squared_length_d_object();
1399 typename K::Scaled_vector_d scaled_vec =
1400 m_k.scaled_vector_d_object();
1401 typename K::Scalar_product_d scalar_pdct =
1402 m_k.scalar_product_d_object();
1403 typename K::Difference_of_vectors_d diff_vec =
1404 m_k.difference_of_vectors_d_object();
1407 Eigen::MatrixXd mat_points(s.size() * num_pts_for_pca, m_ambient_dim);
1408 unsigned int current_row = 0;
1410 for (Simplex::const_iterator it_index = s.begin();
1411 it_index != s.end(); ++it_index) {
1412 const Point &p = m_points[*it_index];
1414 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 1415 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca,
false);
1416 const Points &points_for_pca = m_points_for_tse;
1418 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca,
false);
1419 const Points &points_for_pca = m_points;
1422 auto nn_it = kns_range.begin();
1424 current_row < num_pts_for_pca && nn_it != kns_range.end();
1425 ++current_row, ++nn_it) {
1426 for (
int i = 0; i < m_ambient_dim; ++i) {
1427 mat_points(current_row, i) =
1428 CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1432 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1433 Eigen::MatrixXd cov = centered.adjoint() * centered;
1434 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1436 Tangent_space_basis tsb;
1440 for (
int j = m_ambient_dim - 1;
1441 j >= m_ambient_dim - m_intrinsic_dim;
1443 if (normalize_basis) {
1444 Vector v = constr_vec(m_ambient_dim,
1445 eig.eigenvectors().col(j).data(),
1446 eig.eigenvectors().col(j).data() + m_ambient_dim);
1447 tsb.push_back(normalize_vector(v, m_k));
1449 tsb.push_back(constr_vec(
1451 eig.eigenvectors().col(j).data(),
1452 eig.eigenvectors().col(j).data() + m_ambient_dim));
1461 int tangent_basis_dim(std::size_t i)
const {
1462 return m_tangent_spaces[i].dimension();
1465 Point compute_perturbed_point(std::size_t pt_idx)
const {
1466 #ifdef GUDHI_TC_PERTURB_POSITION 1467 return m_k.translated_point_d_object()(
1468 m_points[pt_idx], m_translations[pt_idx]);
1470 return m_points[pt_idx];
1474 void compute_perturbed_weighted_point(std::size_t pt_idx, Point &p, FT &w)
const {
1475 #ifdef GUDHI_TC_PERTURB_POSITION 1476 p = m_k.translated_point_d_object()(
1477 m_points[pt_idx], m_translations[pt_idx]);
1479 p = m_points[pt_idx];
1481 w = m_weights[pt_idx];
1484 Weighted_point compute_perturbed_weighted_point(std::size_t pt_idx)
const {
1485 typename K::Construct_weighted_point_d k_constr_wp =
1486 m_k.construct_weighted_point_d_object();
1488 Weighted_point wp = k_constr_wp(
1489 #ifdef GUDHI_TC_PERTURB_POSITION
1490 m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]),
1499 Point unproject_point(
const Tr_point &p,
1500 const Tangent_space_basis &tsb,
1501 const Tr_traits &tr_traits)
const {
1502 typename K::Translated_point_d k_transl =
1503 m_k.translated_point_d_object();
1504 typename K::Scaled_vector_d k_scaled_vec =
1505 m_k.scaled_vector_d_object();
1506 typename Tr_traits::Compute_coordinate_d coord =
1507 tr_traits.compute_coordinate_d_object();
1509 Point global_point = compute_perturbed_point(tsb.origin());
1510 for (
int i = 0; i < m_intrinsic_dim; ++i)
1511 global_point = k_transl(global_point,
1512 k_scaled_vec(tsb[i], coord(p, i)));
1514 return global_point;
1519 Tr_bare_point project_point(
const Point &p,
1520 const Tangent_space_basis &tsb,
1521 const Tr_traits &tr_traits)
const {
1522 typename K::Scalar_product_d scalar_pdct =
1523 m_k.scalar_product_d_object();
1524 typename K::Difference_of_points_d diff_points =
1525 m_k.difference_of_points_d_object();
1527 Vector v = diff_points(p, compute_perturbed_point(tsb.origin()));
1529 std::vector<FT> coords;
1531 coords.reserve(tsb.dimension());
1532 for (std::size_t i = 0; i < m_intrinsic_dim; ++i) {
1534 FT coord = scalar_pdct(v, tsb[i]);
1535 coords.push_back(coord);
1538 return tr_traits.construct_point_d_object()(
1539 static_cast<int> (coords.size()), coords.begin(), coords.end());
1546 Tr_point project_point_and_compute_weight(
const Weighted_point &wp,
1547 const Tangent_space_basis &tsb,
1548 const Tr_traits &tr_traits)
const {
1549 typename K::Point_drop_weight_d k_drop_w =
1550 m_k.point_drop_weight_d_object();
1551 typename K::Compute_weight_d k_point_weight =
1552 m_k.compute_weight_d_object();
1553 return project_point_and_compute_weight(
1554 k_drop_w(wp), k_point_weight(wp), tsb, tr_traits);
1558 Tr_point project_point_and_compute_weight(
const Point &p,
const FT w,
1559 const Tangent_space_basis &tsb,
1560 const Tr_traits &tr_traits)
const {
1561 const int point_dim = m_k.point_dimension_d_object()(p);
1563 typename K::Construct_point_d constr_pt =
1564 m_k.construct_point_d_object();
1565 typename K::Scalar_product_d scalar_pdct =
1566 m_k.scalar_product_d_object();
1567 typename K::Difference_of_points_d diff_points =
1568 m_k.difference_of_points_d_object();
1569 typename K::Compute_coordinate_d coord =
1570 m_k.compute_coordinate_d_object();
1571 typename K::Construct_cartesian_const_iterator_d ccci =
1572 m_k.construct_cartesian_const_iterator_d_object();
1574 Point origin = compute_perturbed_point(tsb.origin());
1575 Vector v = diff_points(p, origin);
1578 bool same_dim = (point_dim == tsb.dimension());
1580 std::vector<FT> coords;
1582 std::vector<FT> p_proj(ccci(origin), ccci(origin, 0));
1583 coords.reserve(tsb.dimension());
1584 for (
int i = 0; i < tsb.dimension(); ++i) {
1586 FT c = scalar_pdct(v, tsb[i]);
1587 coords.push_back(c);
1591 for (
int j = 0; j < point_dim; ++j)
1592 p_proj[j] += c * coord(tsb[i], j);
1597 FT sq_dist_to_proj_pt = 0;
1599 Point projected_pt = constr_pt(point_dim, p_proj.begin(), p_proj.end());
1600 sq_dist_to_proj_pt = m_k.squared_distance_d_object()(p, projected_pt);
1603 return tr_traits.construct_weighted_point_d_object()
1604 (tr_traits.construct_point_d_object()(
static_cast<int> (coords.size()), coords.begin(), coords.end()),
1605 w - sq_dist_to_proj_pt);
1610 template <
typename Indexed_po
int_range>
1611 std::vector<Tr_point> project_points_and_compute_weights(
1612 const Indexed_point_range &point_indices,
1613 const Tangent_space_basis &tsb,
1614 const Tr_traits &tr_traits)
const {
1615 std::vector<Tr_point> ret;
1616 for (
typename Indexed_point_range::const_iterator
1617 it = point_indices.begin(), it_end = point_indices.end();
1618 it != it_end; ++it) {
1619 ret.push_back(project_point_and_compute_weight(
1620 compute_perturbed_weighted_point(*it), tsb, tr_traits));
1627 bool is_simplex_consistent(Tr_full_cell_handle fch,
int cur_dim)
const {
1629 for (
int i = 0; i < cur_dim + 1; ++i) {
1630 std::size_t data = fch->vertex(i)->data();
1633 return is_simplex_consistent(c);
1639 bool is_simplex_consistent(Simplex
const& simplex)
const {
1641 Simplex::const_iterator it_point_idx = simplex.begin();
1644 for (; it_point_idx != simplex.end(); ++it_point_idx) {
1645 std::size_t point_idx = *it_point_idx;
1647 if (point_idx == (std::numeric_limits<std::size_t>::max)())
1650 Star
const& star = m_stars[point_idx];
1653 Incident_simplex is_to_find = simplex;
1654 is_to_find.erase(point_idx);
1657 if (std::find(star.begin(), star.end(), is_to_find) == star.end())
1670 template <
typename OutputIterator>
1671 bool is_simplex_consistent(
1672 std::size_t center_point,
1673 Incident_simplex
const& s,
1674 OutputIterator points_whose_star_does_not_contain_s,
1675 bool check_also_in_non_maximal_faces =
false)
const {
1676 Simplex full_simplex = s;
1677 full_simplex.insert(center_point);
1680 Incident_simplex::const_iterator it_point_idx = s.begin();
1683 for (; it_point_idx != s.end(); ++it_point_idx) {
1684 std::size_t point_idx = *it_point_idx;
1686 if (point_idx == (std::numeric_limits<std::size_t>::max)())
1689 Star
const& star = m_stars[point_idx];
1692 Incident_simplex is_to_find = full_simplex;
1693 is_to_find.erase(point_idx);
1695 if (check_also_in_non_maximal_faces) {
1699 for (Star::const_iterator is = star.begin(), is_end = star.end();
1700 !found && is != is_end; ++is) {
1701 if (std::includes(is->begin(), is->end(),
1702 is_to_find.begin(), is_to_find.end()))
1707 *points_whose_star_does_not_contain_s++ = point_idx;
1710 if (std::find(star.begin(), star.end(), is_to_find) == star.end())
1711 *points_whose_star_does_not_contain_s++ = point_idx;
1721 bool is_simplex_in_star(std::size_t p,
1722 Incident_simplex
const& s,
1723 bool check_also_in_non_maximal_faces =
true)
const {
1724 Star
const& star = m_stars[p];
1726 if (check_also_in_non_maximal_faces) {
1730 for (Star::const_iterator is = star.begin(), is_end = star.end();
1731 !found && is != is_end; ++is) {
1732 if (std::includes(is->begin(), is->end(), s.begin(), s.end()))
1738 return !(std::find(star.begin(), star.end(), s) == star.end());
1742 #ifdef GUDHI_USE_TBB 1744 class Try_to_solve_inconsistencies_in_a_local_triangulation {
1746 double m_max_perturb;
1747 tbb::combinable<std::size_t> &m_num_inconsistencies;
1748 tbb::combinable<std::vector<std::size_t> > &m_updated_points;
1754 tbb::combinable<std::size_t> &num_inconsistencies,
1755 tbb::combinable<std::vector<std::size_t> > &updated_points)
1757 m_max_perturb(max_perturb),
1758 m_num_inconsistencies(num_inconsistencies),
1759 m_updated_points(updated_points) { }
1762 Try_to_solve_inconsistencies_in_a_local_triangulation(
const Try_to_solve_inconsistencies_in_a_local_triangulation&
1765 m_max_perturb(tsilt.m_max_perturb),
1766 m_num_inconsistencies(tsilt.m_num_inconsistencies),
1767 m_updated_points(tsilt.m_updated_points) { }
1770 void operator()(
const tbb::blocked_range<size_t>& r)
const {
1771 for (
size_t i = r.begin(); i != r.end(); ++i) {
1772 m_num_inconsistencies.local() +=
1773 m_tc.try_to_solve_inconsistencies_in_a_local_triangulation(i, m_max_perturb,
1774 std::back_inserter(m_updated_points.local()));
1778 #endif // GUDHI_USE_TBB 1780 void perturb(std::size_t point_idx,
double max_perturb) {
1781 const Tr_traits &local_tr_traits =
1782 m_triangulations[point_idx].tr().geom_traits();
1783 typename Tr_traits::Compute_coordinate_d coord =
1784 local_tr_traits.compute_coordinate_d_object();
1785 typename K::Translated_point_d k_transl =
1786 m_k.translated_point_d_object();
1787 typename K::Construct_vector_d k_constr_vec =
1788 m_k.construct_vector_d_object();
1789 typename K::Scaled_vector_d k_scaled_vec =
1790 m_k.scaled_vector_d_object();
1792 CGAL::Random_points_in_ball_d<Tr_bare_point>
1793 tr_point_in_ball_generator(m_intrinsic_dim,
1794 m_random_generator.get_double(0., max_perturb));
1796 Tr_point local_random_transl =
1797 local_tr_traits.construct_weighted_point_d_object()(*tr_point_in_ball_generator++, 0);
1798 Translation_for_perturb global_transl = k_constr_vec(m_ambient_dim);
1799 const Tangent_space_basis &tsb = m_tangent_spaces[point_idx];
1800 for (
int i = 0; i < m_intrinsic_dim; ++i) {
1801 global_transl = k_transl(global_transl,
1802 k_scaled_vec(tsb[i], coord(local_random_transl, i)));
1805 #if defined(GUDHI_USE_TBB) 1806 m_p_perturb_mutexes[point_idx].lock();
1807 m_translations[point_idx] = global_transl;
1808 m_p_perturb_mutexes[point_idx].unlock();
1811 m_translations[point_idx] = global_transl;
1816 template <
typename OutputIt>
1817 bool try_to_solve_inconsistencies_in_a_local_triangulation(std::size_t tr_index,
1819 OutputIt perturbed_pts_indices = CGAL::Emptyset_iterator()) {
1820 bool is_inconsistent =
false;
1822 Star
const& star = m_stars[tr_index];
1825 Star::const_iterator it_inc_simplex = star.begin();
1826 Star::const_iterator it_inc_simplex_end = star.end();
1827 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1828 const Incident_simplex &incident_simplex = *it_inc_simplex;
1831 if (is_infinite(incident_simplex))
1834 Simplex c = incident_simplex;
1838 if (!is_simplex_consistent(c)) {
1839 is_inconsistent =
true;
1841 std::size_t idx = tr_index;
1843 perturb(tr_index, max_perturb);
1844 *perturbed_pts_indices++ = idx;
1851 return is_inconsistent;
1857 std::ostream &export_point_set(std::ostream & os,
1858 bool use_perturbed_points =
false,
1859 const char *coord_separator =
" ")
const {
1860 if (use_perturbed_points) {
1861 std::vector<Point> perturbed_points;
1862 perturbed_points.reserve(m_points.size());
1863 for (std::size_t i = 0; i < m_points.size(); ++i)
1864 perturbed_points.push_back(compute_perturbed_point(i));
1866 return export_point_set(
1867 m_k, perturbed_points, os, coord_separator);
1869 return export_point_set(
1870 m_k, m_points, os, coord_separator);
1874 template<
typename ProjectionFunctor = CGAL::Identity<Po
int> >
1875 std::ostream &export_vertices_to_off(
1876 std::ostream & os, std::size_t &num_vertices,
1877 bool use_perturbed_points =
false,
1878 ProjectionFunctor
const& point_projection = ProjectionFunctor())
const {
1879 if (m_points.empty()) {
1887 const int N = (m_intrinsic_dim == 1 ? 2 : 1);
1890 typename K::Compute_coordinate_d coord =
1891 m_k.compute_coordinate_d_object();
1893 #ifdef GUDHI_TC_EXPORT_ALL_COORDS_IN_OFF 1894 int num_coords = m_ambient_dim;
1896 int num_coords = (std::min)(m_ambient_dim, 3);
1899 #ifdef GUDHI_TC_EXPORT_NORMALS 1900 OS_container::const_iterator it_os = m_orth_spaces.begin();
1902 typename Points::const_iterator it_p = m_points.begin();
1903 typename Points::const_iterator it_p_end = m_points.end();
1905 for (std::size_t i = 0; it_p != it_p_end; ++it_p, ++i) {
1906 Point p = point_projection(
1907 use_perturbed_points ? compute_perturbed_point(i) : *it_p);
1908 for (
int ii = 0; ii < N; ++ii) {
1910 for (; j < num_coords; ++j)
1911 os << CGAL::to_double(coord(p, j)) <<
" ";
1915 #ifdef GUDHI_TC_EXPORT_NORMALS 1916 for (j = 0; j < num_coords; ++j)
1917 os <<
" " << CGAL::to_double(coord(*it_os->begin(), j));
1921 #ifdef GUDHI_TC_EXPORT_NORMALS 1926 num_vertices = N * m_points.size();
1930 std::ostream &export_simplices_to_off(std::ostream & os, std::size_t &num_OFF_simplices,
1931 bool color_inconsistencies =
false,
1932 Simplex_set
const *p_simpl_to_color_in_red = NULL,
1933 Simplex_set
const *p_simpl_to_color_in_green = NULL,
1934 Simplex_set
const *p_simpl_to_color_in_blue = NULL)
1938 num_OFF_simplices = 0;
1939 std::size_t num_maximal_simplices = 0;
1940 std::size_t num_inconsistent_maximal_simplices = 0;
1941 std::size_t num_inconsistent_stars = 0;
1942 typename Tr_container::const_iterator it_tr = m_triangulations.begin();
1943 typename Tr_container::const_iterator it_tr_end = m_triangulations.end();
1945 for (std::size_t idx = 0; it_tr != it_tr_end; ++it_tr, ++idx) {
1946 bool is_star_inconsistent =
false;
1948 Triangulation
const& tr = it_tr->tr();
1950 if (tr.current_dimension() < m_intrinsic_dim)
1954 std::stringstream color;
1956 color << 128 <<
" " << 128 <<
" " << 128;
1959 typedef std::vector<std::pair<Simplex, int> > Star_using_triangles;
1960 Star_using_triangles star_using_triangles;
1963 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
1964 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
1965 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1966 Simplex c = *it_inc_simplex;
1968 std::size_t num_vertices = c.size();
1969 ++num_maximal_simplices;
1971 int color_simplex = -1;
1972 if (color_inconsistencies && !is_simplex_consistent(c)) {
1973 ++num_inconsistent_maximal_simplices;
1975 is_star_inconsistent =
true;
1977 if (p_simpl_to_color_in_red &&
1979 p_simpl_to_color_in_red->begin(),
1980 p_simpl_to_color_in_red->end(),
1981 c) != p_simpl_to_color_in_red->end()) {
1983 }
else if (p_simpl_to_color_in_green &&
1985 p_simpl_to_color_in_green->begin(),
1986 p_simpl_to_color_in_green->end(),
1987 c) != p_simpl_to_color_in_green->end()) {
1989 }
else if (p_simpl_to_color_in_blue &&
1991 p_simpl_to_color_in_blue->begin(),
1992 p_simpl_to_color_in_blue->end(),
1993 c) != p_simpl_to_color_in_blue->end()) {
2002 if (m_intrinsic_dim == 1) {
2004 Simplex::iterator it = c.begin();
2005 for (; it != c.end(); ++it)
2006 tmp_c.insert(*it * 2);
2007 if (num_vertices == 2)
2008 tmp_c.insert(*tmp_c.rbegin() + 1);
2013 if (num_vertices <= 3) {
2014 star_using_triangles.push_back(std::make_pair(c, color_simplex));
2017 std::vector<bool> booleans(num_vertices,
false);
2018 std::fill(booleans.begin() + num_vertices - 3, booleans.end(),
true);
2021 Simplex::iterator it = c.begin();
2022 for (
int i = 0; it != c.end(); ++i, ++it) {
2024 triangle.insert(*it);
2026 star_using_triangles.push_back(
2027 std::make_pair(triangle, color_simplex));
2028 }
while (std::next_permutation(booleans.begin(), booleans.end()));
2033 Star_using_triangles::const_iterator it_simplex =
2034 star_using_triangles.begin();
2035 Star_using_triangles::const_iterator it_simplex_end =
2036 star_using_triangles.end();
2037 for (; it_simplex != it_simplex_end; ++it_simplex) {
2038 const Simplex &c = it_simplex->first;
2044 int color_simplex = it_simplex->second;
2046 std::stringstream sstr_c;
2048 Simplex::const_iterator it_point_idx = c.begin();
2049 for (; it_point_idx != c.end(); ++it_point_idx) {
2050 sstr_c << *it_point_idx <<
" ";
2053 os << 3 <<
" " << sstr_c.str();
2054 if (color_inconsistencies || p_simpl_to_color_in_red
2055 || p_simpl_to_color_in_green || p_simpl_to_color_in_blue) {
2056 switch (color_simplex) {
2057 case 0: os <<
" 255 255 0";
2059 case 1: os <<
" 255 0 0";
2061 case 2: os <<
" 0 255 0";
2063 case 3: os <<
" 0 0 255";
2065 default: os <<
" " << color.str();
2069 ++num_OFF_simplices;
2072 if (is_star_inconsistent)
2073 ++num_inconsistent_stars;
2078 <<
"\n==========================================================\n" 2079 <<
"Export from list of stars to OFF:\n" 2080 <<
" * Number of vertices: " << m_points.size() <<
"\n" 2081 <<
" * Total number of maximal simplices: " << num_maximal_simplices
2083 if (color_inconsistencies) {
2085 <<
" * Number of inconsistent stars: " 2086 << num_inconsistent_stars <<
" (" 2087 << (m_points.size() > 0 ?
2088 100. * num_inconsistent_stars / m_points.size() : 0.) <<
"%)\n" 2089 <<
" * Number of inconsistent maximal simplices: " 2090 << num_inconsistent_maximal_simplices <<
" (" 2091 << (num_maximal_simplices > 0 ?
2092 100. * num_inconsistent_maximal_simplices / num_maximal_simplices
2095 std::cerr <<
"==========================================================\n";
2102 std::ostream &export_simplices_to_off(
2103 const Simplicial_complex &complex,
2104 std::ostream & os, std::size_t &num_OFF_simplices,
2105 Simplex_set
const *p_simpl_to_color_in_red = NULL,
2106 Simplex_set
const *p_simpl_to_color_in_green = NULL,
2107 Simplex_set
const *p_simpl_to_color_in_blue = NULL)
2109 typedef Simplicial_complex::Simplex Simplex;
2110 typedef Simplicial_complex::Simplex_set Simplex_set;
2114 num_OFF_simplices = 0;
2115 std::size_t num_maximal_simplices = 0;
2117 typename Simplex_set::const_iterator it_s =
2118 complex.simplex_range().begin();
2119 typename Simplex_set::const_iterator it_s_end =
2120 complex.simplex_range().end();
2122 for (; it_s != it_s_end; ++it_s) {
2124 ++num_maximal_simplices;
2126 int color_simplex = -1;
2127 if (p_simpl_to_color_in_red &&
2129 p_simpl_to_color_in_red->begin(),
2130 p_simpl_to_color_in_red->end(),
2131 c) != p_simpl_to_color_in_red->end()) {
2133 }
else if (p_simpl_to_color_in_green &&
2134 std::find(p_simpl_to_color_in_green->begin(),
2135 p_simpl_to_color_in_green->end(),
2136 c) != p_simpl_to_color_in_green->end()) {
2138 }
else if (p_simpl_to_color_in_blue &&
2139 std::find(p_simpl_to_color_in_blue->begin(),
2140 p_simpl_to_color_in_blue->end(),
2141 c) != p_simpl_to_color_in_blue->end()) {
2146 typedef std::vector<Simplex> Triangles;
2147 Triangles triangles;
2149 int num_vertices =
static_cast<int>(c.size());
2151 if (num_vertices < m_intrinsic_dim + 1)
2158 if (m_intrinsic_dim == 1) {
2160 Simplex::iterator it = c.begin();
2161 for (; it != c.end(); ++it)
2162 tmp_c.insert(*it * 2);
2163 if (num_vertices == 2)
2164 tmp_c.insert(*tmp_c.rbegin() + 1);
2169 if (num_vertices <= 3) {
2170 triangles.push_back(c);
2173 std::vector<bool> booleans(num_vertices,
false);
2174 std::fill(booleans.begin() + num_vertices - 3, booleans.end(),
true);
2177 Simplex::iterator it = c.begin();
2178 for (
int i = 0; it != c.end(); ++i, ++it) {
2180 triangle.insert(*it);
2182 triangles.push_back(triangle);
2183 }
while (std::next_permutation(booleans.begin(), booleans.end()));
2187 Triangles::const_iterator it_tri = triangles.begin();
2188 Triangles::const_iterator it_tri_end = triangles.end();
2189 for (; it_tri != it_tri_end; ++it_tri) {
2191 if (is_infinite(*it_tri))
2195 Simplex::const_iterator it_point_idx = it_tri->begin();
2196 for (; it_point_idx != it_tri->end(); ++it_point_idx) {
2197 os << *it_point_idx <<
" ";
2200 if (p_simpl_to_color_in_red || p_simpl_to_color_in_green
2201 || p_simpl_to_color_in_blue) {
2202 switch (color_simplex) {
2203 case 0: os <<
" 255 255 0";
2205 case 1: os <<
" 255 0 0";
2207 case 2: os <<
" 0 255 0";
2209 case 3: os <<
" 0 0 255";
2211 default: os <<
" 128 128 128";
2216 ++num_OFF_simplices;
2223 <<
"\n==========================================================\n" 2224 <<
"Export from complex to OFF:\n" 2225 <<
" * Number of vertices: " << m_points.size() <<
"\n" 2226 <<
" * Total number of maximal simplices: " << num_maximal_simplices
2228 <<
"==========================================================\n";
2236 const int m_intrinsic_dim;
2237 const int m_ambient_dim;
2241 #ifdef GUDHI_TC_PERTURB_POSITION 2242 Translations_for_perturb m_translations;
2243 #if defined(GUDHI_USE_TBB) 2244 Mutex_for_perturb *m_p_perturb_mutexes;
2248 Points_ds m_points_ds;
2249 double m_last_max_perturb;
2250 std::vector<bool> m_are_tangent_spaces_computed;
2251 TS_container m_tangent_spaces;
2252 #ifdef GUDHI_TC_EXPORT_NORMALS 2253 OS_container m_orth_spaces;
2255 Tr_container m_triangulations;
2257 Stars_container m_stars;
2258 std::vector<FT> m_squared_star_spheres_radii_incl_margin;
2260 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 2261 Points m_points_for_tse;
2262 Points_ds m_points_ds_for_tse;
2265 mutable CGAL::Random m_random_generator;
2271 #endif // TANGENTIAL_COMPLEX_H_
std::size_t num_inconsistent_simplices
Number of inconsistent simplices.
Definition: Tangential_complex.h:610
Type returned by Tangential_complex::fix_inconsistencies_using_perturbation.
Definition: Tangential_complex.h:419
Point get_point(std::size_t vertex) const
Returns the point corresponding to the vertex given as parameter.
Definition: Tangential_complex.h:322
Tangential complex data structure.
Definition: Tangential_complex.h:132
Incremental_neighbor_search INS_range
The range returned by an incremental nearest or furthest neighbor search. Its value type is std::pair...
Definition: Kd_tree_search.h:110
K_neighbor_search KNS_range
The range returned by a k-nearest or k-furthest neighbor search. Its value type is std::pair<std::siz...
Definition: Kd_tree_search.h:102
KNS_range k_nearest_neighbors(Point const &p, unsigned int k, bool sorted=true, FT eps=FT(0)) const
Search for the k-nearest neighbors from a query point.
Definition: Kd_tree_search.h:174
std::size_t final_num_inconsistent_stars
final number of inconsistent stars
Definition: Tangential_complex.h:429
Definition: SimplicialComplexForAlpha.h:26
~Tangential_complex()
Destructor.
Definition: Tangential_complex.h:297
std::size_t num_simplices
Total number of simplices in stars (including duplicates that appear in several stars) ...
Definition: Tangential_complex.h:608
std::size_t best_num_inconsistent_stars
best number of inconsistent stars during the process
Definition: Tangential_complex.h:427
std::size_t num_inconsistent_stars
Number of stars containing at least one inconsistent simplex.
Definition: Tangential_complex.h:612
int create_complex(Simplex_tree_ &tree, bool export_inconsistent_simplices=true) const
Exports the complex into a Simplex_tree.
Definition: Tangential_complex.h:682
Fix_inconsistencies_info fix_inconsistencies_using_perturbation(double max_perturb, double time_limit=-1.)
Attempts to fix inconsistencies by perturbing the point positions.
Definition: Tangential_complex.h:437
std::size_t initial_num_inconsistent_stars
initial number of inconsistent stars
Definition: Tangential_complex.h:425
Tangential_complex(Point_range points, int intrinsic_dimension, const K &k=K())
Constructor from a range of points. Points are copied into the instance, and a search data structure ...
Definition: Tangential_complex.h:268
Num_inconsistencies number_of_inconsistent_simplices(bool verbose=false) const
Definition: Tangential_complex.h:619
bool success
true if all inconsistencies could be removed, false if the time limit has been reached before ...
Definition: Tangential_complex.h:421
int intrinsic_dimension() const
Returns the intrinsic dimension of the manifold.
Definition: Tangential_complex.h:304
Definition: Intro_spatial_searching.h:28
void compute_tangential_complex()
Computes the tangential complex.
Definition: Tangential_complex.h:369
int ambient_dimension() const
Returns the ambient dimension.
Definition: Tangential_complex.h:309
Type returned by Tangential_complex::number_of_inconsistent_simplices.
Definition: Tangential_complex.h:606
std::size_t number_of_vertices() const
Returns the number of vertices.
Definition: Tangential_complex.h:337
unsigned int num_steps
number of steps performed
Definition: Tangential_complex.h:423
Point get_perturbed_point(std::size_t vertex) const
Returns the perturbed position of the point corresponding to the vertex given as parameter.
Definition: Tangential_complex.h:331