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)()) : m_data(data) {}
88 operator std::size_t() {
return m_data; }
90 operator std::size_t()
const {
return m_data; }
121 template <
typename Kernel_,
122 typename DimensionTag,
123 typename Concurrency_tag = CGAL::Parallel_tag,
typename Triangulation_ = CGAL::Default>
126 typedef typename K::FT FT;
127 typedef typename K::Point_d Point;
128 typedef typename K::Weighted_point_d Weighted_point;
129 typedef typename K::Vector_d Vector;
131 typedef typename CGAL::Default::Get<
133 CGAL::Regular_triangulation<
134 CGAL::Epick_d<DimensionTag>,
135 CGAL::Triangulation_data_structure<
136 typename CGAL::Epick_d<DimensionTag>::Dimension,
137 CGAL::Triangulation_vertex<CGAL::Regular_triangulation_traits_adapter<CGAL::Epick_d<DimensionTag> >,
139 CGAL::Triangulation_full_cell<
140 CGAL::Regular_triangulation_traits_adapter<CGAL::Epick_d<DimensionTag> > > > > >::type Triangulation;
141 typedef typename Triangulation::Geom_traits Tr_traits;
142 typedef typename Triangulation::Weighted_point Tr_point;
143 typedef typename Tr_traits::Base::Point_d Tr_bare_point;
144 typedef typename Triangulation::Vertex_handle Tr_vertex_handle;
145 typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle;
146 typedef typename Tr_traits::Vector_d Tr_vector;
148 #if defined(GUDHI_USE_TBB) 149 typedef tbb::mutex Mutex_for_perturb;
150 typedef Vector Translation_for_perturb;
151 typedef std::vector<Atomic_wrapper<FT> > Weights;
153 typedef Vector Translation_for_perturb;
154 typedef std::vector<FT> Weights;
156 typedef std::vector<Translation_for_perturb> Translations_for_perturb;
162 Tr_and_VH() : m_tr(NULL) {}
164 Tr_and_VH(
int dim) : m_tr(
new Triangulation(dim)) {}
166 ~Tr_and_VH() { destroy_triangulation(); }
168 Triangulation &construct_triangulation(
int dim) {
170 m_tr =
new Triangulation(dim);
174 void destroy_triangulation() {
179 Triangulation &tr() {
return *m_tr; }
181 Triangulation
const &tr()
const {
return *m_tr; }
183 Tr_vertex_handle
const ¢er_vertex()
const {
return m_center_vertex; }
185 Tr_vertex_handle ¢er_vertex() {
return m_center_vertex; }
189 Tr_vertex_handle m_center_vertex;
193 typedef Basis<K> Tangent_space_basis;
194 typedef Basis<K> Orthogonal_space_basis;
195 typedef std::vector<Tangent_space_basis> TS_container;
196 typedef std::vector<Orthogonal_space_basis> OS_container;
198 typedef std::vector<Point> Points;
200 typedef boost::container::flat_set<std::size_t> Simplex;
201 typedef std::set<Simplex> Simplex_set;
208 typedef std::vector<Tr_and_VH> Tr_container;
209 typedef std::vector<Vector> Vectors;
213 typedef boost::container::flat_set<std::size_t> Incident_simplex;
214 typedef std::vector<Incident_simplex> Star;
215 typedef std::vector<Star> Stars_container;
219 static const Tr_point &vertex_handle_to_point(Tr_vertex_handle vh) {
return vh->point(); }
221 template <
typename P,
typename VH>
222 static const P &vertex_handle_to_point(VH vh) {
227 typedef internal::Simplicial_complex Simplicial_complex;
238 template <
typename Po
int_range>
240 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
241 InputIterator first_for_tse, InputIterator last_for_tse,
245 m_intrinsic_dim(intrinsic_dimension),
246 m_ambient_dim(points.empty() ? 0 : k.point_dimension_d_object()(*points.begin())),
247 m_points(points.begin(), points.end()),
248 m_weights(m_points.size(), FT(0))
249 #if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION)
251 m_p_perturb_mutexes(NULL)
254 m_points_ds(m_points),
255 m_last_max_perturb(0.),
256 m_are_tangent_spaces_computed(m_points.size(), false),
257 m_tangent_spaces(m_points.size(), Tangent_space_basis())
258 #ifdef GUDHI_TC_EXPORT_NORMALS
260 m_orth_spaces(m_points.size(), Orthogonal_space_basis())
262 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
264 m_points_for_tse(first_for_tse, last_for_tse),
265 m_points_ds_for_tse(m_points_for_tse)
272 #if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION) 273 delete[] m_p_perturb_mutexes;
283 Points
const &points()
const {
return m_points; }
290 Point
get_point(std::size_t vertex)
const {
return m_points[vertex]; }
303 void set_weights(
const Weights &weights) { m_weights = weights; }
305 void set_tangent_planes(
const TS_container &tangent_spaces
306 #ifdef GUDHI_TC_EXPORT_NORMALS
308 const OS_container &orthogonal_spaces
311 #ifdef GUDHI_TC_EXPORT_NORMALS 312 GUDHI_CHECK(m_points.size() == tangent_spaces.size() && m_points.size() == orthogonal_spaces.size(),
313 std::logic_error(
"Wrong sizes"));
315 GUDHI_CHECK(m_points.size() == tangent_spaces.size(), std::logic_error(
"Wrong sizes"));
317 m_tangent_spaces = tangent_spaces;
318 #ifdef GUDHI_TC_EXPORT_NORMALS 319 m_orth_spaces = orthogonal_spaces;
321 for (std::size_t i = 0; i < m_points.size(); ++i) m_are_tangent_spaces_computed[i] =
true;
326 #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS 327 std::cerr << red <<
"WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. " 328 <<
"Computation might be slower than usual.\n" 332 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB) 339 m_triangulations.resize(m_points.size());
340 m_stars.resize(m_points.size());
341 m_squared_star_spheres_radii_incl_margin.resize(m_points.size(), FT(-1));
342 #ifdef GUDHI_TC_PERTURB_POSITION 343 if (m_points.empty())
344 m_translations.clear();
346 m_translations.resize(m_points.size(), m_k.construct_vector_d_object()(m_ambient_dim));
347 #if defined(GUDHI_USE_TBB) 348 delete[] m_p_perturb_mutexes;
349 m_p_perturb_mutexes =
new Mutex_for_perturb[m_points.size()];
355 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
356 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*
this));
358 #endif // GUDHI_USE_TBB 360 for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
363 #endif // GUDHI_USE_TBB 365 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB) 367 std::cerr <<
"Tangential complex computed in " << t.num_seconds() <<
" seconds.\n";
374 bool success =
false;
376 unsigned int num_steps = 0;
378 std::size_t initial_num_inconsistent_stars = 0;
380 std::size_t best_num_inconsistent_stars = 0;
382 std::size_t final_num_inconsistent_stars = 0;
393 if (time_limit == 0.)
return info;
397 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 398 std::tuple<std::size_t, std::size_t, std::size_t> stats_before = number_of_inconsistent_simplices(
false);
400 if (std::get<1>(stats_before) == 0) {
402 std::cerr <<
"Nothing to fix.\n";
407 #endif // GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 409 m_last_max_perturb = max_perturb;
415 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 416 std::cerr <<
"\nBefore fix step:\n" 417 <<
" * Total number of simplices in stars (incl. duplicates): " << std::get<0>(stats_before) <<
"\n" 418 <<
" * Num inconsistent simplices in stars (incl. duplicates): " << red << std::get<1>(stats_before)
419 << white <<
" (" << 100. * std::get<1>(stats_before) / std::get<0>(stats_before) <<
"%)\n" 420 <<
" * Number of stars containing inconsistent simplices: " << red << std::get<2>(stats_before)
421 << white <<
" (" << 100. * std::get<2>(stats_before) / m_points.size() <<
"%)\n";
424 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 425 std::cerr << yellow <<
"\nAttempt to fix inconsistencies using perturbations - step #" << info.
num_steps + 1
429 std::size_t num_inconsistent_stars = 0;
430 std::vector<std::size_t> updated_points;
432 #ifdef GUDHI_TC_PROFILING 433 Gudhi::Clock t_fix_step;
437 #if defined(GUDHI_USE_TBB) 438 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
439 tbb::combinable<std::size_t> num_inconsistencies;
440 tbb::combinable<std::vector<std::size_t> > tls_updated_points;
441 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_triangulations.size()),
442 Try_to_solve_inconsistencies_in_a_local_triangulation(*
this, max_perturb, num_inconsistencies,
443 tls_updated_points));
444 num_inconsistent_stars = num_inconsistencies.combine(std::plus<std::size_t>());
446 tls_updated_points.combine([](std::vector<std::size_t>
const &x, std::vector<std::size_t>
const &y) {
447 std::vector<std::size_t> res;
448 res.reserve(x.size() + y.size());
449 res.insert(res.end(), x.begin(), x.end());
450 res.insert(res.end(), y.begin(), y.end());
454 #endif // GUDHI_USE_TBB 456 for (std::size_t i = 0; i < m_triangulations.size(); ++i) {
457 num_inconsistent_stars +=
458 try_to_solve_inconsistencies_in_a_local_triangulation(i, max_perturb, std::back_inserter(updated_points));
460 #if defined(GUDHI_USE_TBB) 462 #endif // GUDHI_USE_TBB 464 #ifdef GUDHI_TC_PROFILING 468 #if defined(GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES) || defined(DEBUG_TRACES) 469 std::cerr <<
"\nEncountered during fix:\n" 470 <<
" * Num stars containing inconsistent simplices: " << red << num_inconsistent_stars << white <<
" (" 471 << 100. * num_inconsistent_stars / m_points.size() <<
"%)\n";
474 #ifdef GUDHI_TC_PROFILING 475 std::cerr << yellow <<
"done in " << t_fix_step.num_seconds() <<
" seconds.\n" << white;
476 #elif defined(DEBUG_TRACES) 477 std::cerr << yellow <<
"done.\n" << white;
480 if (num_inconsistent_stars > 0) refresh_tangential_complex(updated_points);
482 #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS 484 std::size_t num_inc_1 = std::get<1>(number_of_inconsistent_simplices(
false));
485 refresh_tangential_complex();
486 std::size_t num_inc_2 = std::get<1>(number_of_inconsistent_simplices(
false));
487 if (num_inc_1 != num_inc_2)
488 std::cerr << red <<
"REFRESHMENT CHECK: FAILED. (" << num_inc_1 <<
" vs " << num_inc_2 <<
")\n" << white;
490 std::cerr << green <<
"REFRESHMENT CHECK: PASSED.\n" << white;
493 #ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES 494 std::tuple<std::size_t, std::size_t, std::size_t> stats_after = number_of_inconsistent_simplices(
false);
496 std::cerr <<
"\nAfter fix:\n" 497 <<
" * Total number of simplices in stars (incl. duplicates): " << std::get<0>(stats_after) <<
"\n" 498 <<
" * Num inconsistent simplices in stars (incl. duplicates): " << red << std::get<1>(stats_after)
499 << white <<
" (" << 100. * std::get<1>(stats_after) / std::get<0>(stats_after) <<
"%)\n" 500 <<
" * Number of stars containing inconsistent simplices: " << red << std::get<2>(stats_after) << white
501 <<
" (" << 100. * std::get<2>(stats_after) / m_points.size() <<
"%)\n";
503 stats_before = stats_after;
513 done = (num_inconsistent_stars == 0);
516 if (time_limit > 0. && t.num_seconds() > time_limit) {
518 std::cerr << red <<
"Time limit reached.\n" << white;
527 std::cerr << green <<
"Fixed!\n" << white;
536 std::size_t num_simplices = 0;
538 std::size_t num_inconsistent_simplices = 0;
540 std::size_t num_inconsistent_stars = 0;
556 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
557 bool is_star_inconsistent =
false;
560 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
561 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
562 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
564 if (is_infinite(*it_inc_simplex))
continue;
566 Simplex c = *it_inc_simplex;
569 if (!is_simplex_consistent(c)) {
571 is_star_inconsistent =
true;
580 std::cerr <<
"\n==========================================================\n" 581 <<
"Inconsistencies:\n" 582 <<
" * Total number of simplices in stars (incl. duplicates): " << stats.
num_simplices <<
"\n" 583 <<
" * Number of inconsistent simplices in stars (incl. duplicates): " 588 <<
"==========================================================\n";
604 template <
typename Simplex_tree_>
606 bool export_inconsistent_simplices =
true 609 bool export_infinite_simplices =
false, Simplex_set *p_inconsistent_simplices = NULL
612 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 613 std::cerr << yellow <<
"\nExporting the TC as a Simplex_tree... " << white;
615 #ifdef GUDHI_TC_PROFILING 622 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
624 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
625 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
626 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
627 Simplex c = *it_inc_simplex;
630 if (!export_infinite_simplices && is_infinite(c))
continue;
632 if (!export_inconsistent_simplices && !is_simplex_consistent(c))
continue;
634 if (static_cast<int>(c.size()) > max_dim) max_dim =
static_cast<int>(c.size());
639 bool inserted = tree.insert_simplex_and_subfaces(c).second;
642 if (p_inconsistent_simplices && inserted && !is_simplex_consistent(c)) {
643 p_inconsistent_simplices->insert(c);
648 #ifdef GUDHI_TC_PROFILING 650 std::cerr << yellow <<
"done in " << t.num_seconds() <<
" seconds.\n" << white;
651 #elif defined(DEBUG_TRACES) 652 std::cerr << yellow <<
"done.\n" << white;
668 int create_complex(Simplicial_complex &complex,
bool export_inconsistent_simplices =
true,
669 bool export_infinite_simplices =
false,
int check_lower_and_higher_dim_simplices = 2,
670 Simplex_set *p_inconsistent_simplices = NULL)
const {
671 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 672 std::cerr << yellow <<
"\nExporting the TC as a Simplicial_complex... " << white;
674 #ifdef GUDHI_TC_PROFILING 682 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
684 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
685 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
686 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
687 Simplex c = *it_inc_simplex;
690 if (!export_infinite_simplices && is_infinite(c))
continue;
692 if (!export_inconsistent_simplices && !is_simplex_consistent(c))
continue;
695 if (check_lower_and_higher_dim_simplices == 2 && max_dim != -1 && static_cast<int>(c.size()) != max_dim) {
698 <<
"Info: check_lower_and_higher_dim_simplices ACTIVATED. " 699 "Export might be take some time...\n" 701 check_lower_and_higher_dim_simplices = 1;
704 if (static_cast<int>(c.size()) > max_dim) max_dim =
static_cast<int>(c.size());
709 bool added = complex.add_simplex(c, check_lower_and_higher_dim_simplices == 1);
712 if (p_inconsistent_simplices && added && !is_simplex_consistent(c)) {
713 p_inconsistent_simplices->insert(c);
718 #ifdef GUDHI_TC_PROFILING 720 std::cerr << yellow <<
"done in " << t.num_seconds() <<
" seconds.\n" << white;
721 #elif defined(DEBUG_TRACES) 722 std::cerr << yellow <<
"done.\n" << white;
728 template <
typename ProjectionFunctor = CGAL::Identity<Po
int> >
729 std::ostream &export_to_off(
const Simplicial_complex &complex, std::ostream &os,
730 Simplex_set
const *p_simpl_to_color_in_red = NULL,
731 Simplex_set
const *p_simpl_to_color_in_green = NULL,
732 Simplex_set
const *p_simpl_to_color_in_blue = NULL,
733 ProjectionFunctor
const &point_projection = ProjectionFunctor())
const {
734 return export_to_off(os,
false, p_simpl_to_color_in_red, p_simpl_to_color_in_green, p_simpl_to_color_in_blue,
735 &complex, point_projection);
738 template <
typename ProjectionFunctor = CGAL::Identity<Po
int> >
739 std::ostream &export_to_off(std::ostream &os,
bool color_inconsistencies =
false,
740 Simplex_set
const *p_simpl_to_color_in_red = NULL,
741 Simplex_set
const *p_simpl_to_color_in_green = NULL,
742 Simplex_set
const *p_simpl_to_color_in_blue = NULL,
743 const Simplicial_complex *p_complex = NULL,
744 ProjectionFunctor
const &point_projection = ProjectionFunctor())
const {
745 if (m_points.empty())
return os;
747 if (m_ambient_dim < 2) {
748 std::cerr <<
"Error: export_to_off => ambient dimension should be >= 2.\n";
749 os <<
"Error: export_to_off => ambient dimension should be >= 2.\n";
752 if (m_ambient_dim > 3) {
753 std::cerr <<
"Warning: export_to_off => ambient dimension should be " 754 "<= 3. Only the first 3 coordinates will be exported.\n";
757 if (m_intrinsic_dim < 1 || m_intrinsic_dim > 3) {
758 std::cerr <<
"Error: export_to_off => intrinsic dimension should be " 759 "between 1 and 3.\n";
760 os <<
"Error: export_to_off => intrinsic dimension should be " 761 "between 1 and 3.\n";
765 std::stringstream output;
766 std::size_t num_simplices, num_vertices;
767 export_vertices_to_off(output, num_vertices,
false, point_projection);
769 export_simplices_to_off(*p_complex, output, num_simplices, p_simpl_to_color_in_red, p_simpl_to_color_in_green,
770 p_simpl_to_color_in_blue);
772 export_simplices_to_off(output, num_simplices, color_inconsistencies, p_simpl_to_color_in_red,
773 p_simpl_to_color_in_green, p_simpl_to_color_in_blue);
776 #ifdef GUDHI_TC_EXPORT_NORMALS 781 << num_vertices <<
" " << num_simplices <<
" " 789 void refresh_tangential_complex() {
790 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 791 std::cerr << yellow <<
"\nRefreshing TC... " << white;
794 #ifdef GUDHI_TC_PROFILING 799 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
800 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*
this));
802 #endif // GUDHI_USE_TBB 804 for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
807 #endif // GUDHI_USE_TBB 809 #ifdef GUDHI_TC_PROFILING 811 std::cerr << yellow <<
"done in " << t.num_seconds() <<
" seconds.\n" << white;
812 #elif defined(DEBUG_TRACES) 813 std::cerr << yellow <<
"done.\n" << white;
818 template <
typename Po
int_indices_range>
819 void refresh_tangential_complex(Point_indices_range
const &perturbed_points_indices) {
820 #if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING) 821 std::cerr << yellow <<
"\nRefreshing TC... " << white;
824 #ifdef GUDHI_TC_PROFILING 829 Points_ds updated_pts_ds(m_points, perturbed_points_indices);
833 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
834 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()),
835 Refresh_tangent_triangulation(*
this, updated_pts_ds));
837 #endif // GUDHI_USE_TBB 839 for (std::size_t i = 0; i < m_points.size(); ++i) refresh_tangent_triangulation(i, updated_pts_ds);
842 #endif // GUDHI_USE_TBB 844 #ifdef GUDHI_TC_PROFILING 846 std::cerr << yellow <<
"done in " << t.num_seconds() <<
" seconds.\n" << white;
847 #elif defined(DEBUG_TRACES) 848 std::cerr << yellow <<
"done.\n" << white;
852 void export_inconsistent_stars_to_OFF_files(std::string
const &filename_base)
const {
854 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
856 Simplicial_complex sc;
858 bool is_inconsistent =
false;
859 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
860 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
861 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
863 if (is_infinite(*it_inc_simplex))
continue;
865 Simplex c = *it_inc_simplex;
871 if (!is_inconsistent && !is_simplex_consistent(c)) is_inconsistent =
true;
874 if (is_inconsistent) {
876 std::stringstream output_filename;
877 output_filename << filename_base <<
"_" << idx <<
".off";
878 std::ofstream off_stream(output_filename.str().c_str());
879 export_to_off(sc, off_stream);
884 class Compare_distance_to_ref_point {
886 Compare_distance_to_ref_point(Point
const &ref, K
const &k) : m_ref(ref), m_k(k) {}
888 bool operator()(Point
const &p1, Point
const &p2) {
889 typename K::Squared_distance_d sqdist = m_k.squared_distance_d_object();
890 return sqdist(p1, m_ref) < sqdist(p2, m_ref);
900 class Compute_tangent_triangulation {
908 Compute_tangent_triangulation(
const Compute_tangent_triangulation &ctt) : m_tc(ctt.m_tc) {}
911 void operator()(
const tbb::blocked_range<size_t> &r)
const {
912 for (
size_t i = r.begin(); i != r.end(); ++i) m_tc.compute_tangent_triangulation(i);
917 class Refresh_tangent_triangulation {
919 Points_ds
const &m_updated_pts_ds;
923 Refresh_tangent_triangulation(
Tangential_complex &tc, Points_ds
const &updated_pts_ds)
924 : m_tc(tc), m_updated_pts_ds(updated_pts_ds) {}
927 Refresh_tangent_triangulation(
const Refresh_tangent_triangulation &ctt)
928 : m_tc(ctt.m_tc), m_updated_pts_ds(ctt.m_updated_pts_ds) {}
931 void operator()(
const tbb::blocked_range<size_t> &r)
const {
932 for (
size_t i = r.begin(); i != r.end(); ++i) m_tc.refresh_tangent_triangulation(i, m_updated_pts_ds);
935 #endif // GUDHI_USE_TBB 937 bool is_infinite(Simplex
const &s)
const {
return *s.rbegin() == (std::numeric_limits<std::size_t>::max)(); }
942 Tr_vertex_handle compute_star(std::size_t i,
const Point ¢er_pt,
const Tangent_space_basis &tsb,
943 Triangulation &triangulation,
bool verbose =
false) {
944 int tangent_space_dim = tsb.dimension();
945 const Tr_traits &local_tr_traits = triangulation.geom_traits();
948 typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
951 typename Tr_traits::Compute_weight_d point_weight = local_tr_traits.compute_weight_d_object();
952 typename Tr_traits::Power_center_d power_center = local_tr_traits.power_center_d_object();
961 if (i == tsb.origin()) {
963 proj_wp = local_tr_traits.construct_weighted_point_d_object()(
964 local_tr_traits.construct_point_d_object()(tangent_space_dim, CGAL::ORIGIN), m_weights[i]);
966 const Weighted_point &wp = compute_perturbed_weighted_point(i);
967 proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
970 Tr_vertex_handle center_vertex = triangulation.insert(proj_wp);
971 center_vertex->data() = i;
972 if (verbose) std::cerr <<
"* Inserted point #" << i <<
"\n";
974 #ifdef GUDHI_TC_VERY_VERBOSE 975 std::size_t num_attempts_to_insert_points = 1;
976 std::size_t num_inserted_points = 1;
980 INS_range ins_range = m_points_ds.incremental_nearest_neighbors(center_pt);
986 boost::optional<FT> squared_star_sphere_radius_plus_margin = boost::make_optional(
false, FT());
992 for (
auto nn_it = ins_range.begin(); nn_it != ins_range.end(); ++nn_it) {
993 std::size_t neighbor_point_idx = nn_it->first;
996 if (neighbor_point_idx != i) {
1001 compute_perturbed_weighted_point(neighbor_point_idx, neighbor_pt, neighbor_weight);
1003 if (squared_star_sphere_radius_plus_margin &&
1004 k_sqdist(center_pt, neighbor_pt) > *squared_star_sphere_radius_plus_margin)
1007 Tr_point proj_pt = project_point_and_compute_weight(neighbor_pt, neighbor_weight, tsb, local_tr_traits);
1009 #ifdef GUDHI_TC_VERY_VERBOSE 1010 ++num_attempts_to_insert_points;
1013 Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex);
1015 if (vh != Tr_vertex_handle() && vh->data() == (std::numeric_limits<std::size_t>::max)()) {
1016 #ifdef GUDHI_TC_VERY_VERBOSE 1017 ++num_inserted_points;
1019 if (verbose) std::cerr <<
"* Inserted point #" << neighbor_point_idx <<
"\n";
1021 vh->data() = neighbor_point_idx;
1024 if (triangulation.current_dimension() >= tangent_space_dim) {
1025 squared_star_sphere_radius_plus_margin = boost::none;
1027 std::vector<Tr_full_cell_handle> incident_cells;
1028 triangulation.incident_full_cells(center_vertex, std::back_inserter(incident_cells));
1029 for (
typename std::vector<Tr_full_cell_handle>::iterator cit = incident_cells.begin();
1030 cit != incident_cells.end(); ++cit) {
1031 Tr_full_cell_handle cell = *cit;
1032 if (triangulation.is_infinite(cell)) {
1033 squared_star_sphere_radius_plus_margin = boost::none;
1039 power_center(boost::make_transform_iterator(cell->vertices_begin(),
1040 vertex_handle_to_point<Tr_point, Tr_vertex_handle>),
1041 boost::make_transform_iterator(cell->vertices_end(),
1042 vertex_handle_to_point<Tr_point, Tr_vertex_handle>));
1044 FT sq_power_sphere_diam = 4 * point_weight(c);
1046 if (!squared_star_sphere_radius_plus_margin ||
1047 sq_power_sphere_diam > *squared_star_sphere_radius_plus_margin) {
1048 squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
1055 if (squared_star_sphere_radius_plus_margin) {
1057 squared_star_sphere_radius_plus_margin =
1058 CGAL::square(std::sqrt(*squared_star_sphere_radius_plus_margin) + 2 * m_last_max_perturb);
1061 m_squared_star_spheres_radii_incl_margin[i] = *squared_star_sphere_radius_plus_margin;
1063 m_squared_star_spheres_radii_incl_margin[i] = FT(-1);
1070 return center_vertex;
1073 void refresh_tangent_triangulation(std::size_t i, Points_ds
const &updated_pts_ds,
bool verbose =
false) {
1074 if (verbose) std::cerr <<
"** Refreshing tangent tri #" << i <<
" **\n";
1076 if (m_squared_star_spheres_radii_incl_margin[i] == FT(-1))
return compute_tangent_triangulation(i, verbose);
1078 Point center_point = compute_perturbed_point(i);
1080 std::size_t closest_pt_index = updated_pts_ds.
k_nearest_neighbors(center_point, 1,
false).begin()->first;
1082 typename K::Construct_weighted_point_d k_constr_wp = m_k.construct_weighted_point_d_object();
1083 typename K::Power_distance_d k_power_dist = m_k.power_distance_d_object();
1086 Weighted_point star_sphere = k_constr_wp(compute_perturbed_point(i), m_squared_star_spheres_radii_incl_margin[i]);
1087 Weighted_point closest_updated_point = compute_perturbed_weighted_point(closest_pt_index);
1090 if (k_power_dist(star_sphere, closest_updated_point) <= FT(0)) compute_tangent_triangulation(i, verbose);
1093 void compute_tangent_triangulation(std::size_t i,
bool verbose =
false) {
1094 if (verbose) std::cerr <<
"** Computing tangent tri #" << i <<
" **\n";
1099 const Point center_pt = compute_perturbed_point(i);
1100 Tangent_space_basis &tsb = m_tangent_spaces[i];
1103 if (!m_are_tangent_spaces_computed[i]) {
1104 #ifdef GUDHI_TC_EXPORT_NORMALS 1105 tsb = compute_tangent_space(center_pt, i,
true , &m_orth_spaces[i]);
1107 tsb = compute_tangent_space(center_pt, i);
1111 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1114 int tangent_space_dim = tangent_basis_dim(i);
1115 Triangulation &local_tr = m_triangulations[i].construct_triangulation(tangent_space_dim);
1117 m_triangulations[i].center_vertex() = compute_star(i, center_pt, tsb, local_tr, verbose);
1119 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1121 std::cerr <<
" - triangulation construction: " << t.num_seconds() <<
" s.\n";
1125 #ifdef GUDHI_TC_VERY_VERBOSE 1126 std::cerr <<
"Inserted " << num_inserted_points <<
" points / " << num_attempts_to_insert_points
1127 <<
" attemps to compute the star\n";
1132 #if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE) 1134 std::cerr <<
" - update_star: " << t.num_seconds() <<
" s.\n";
1140 void update_star(std::size_t i) {
1141 Star &star = m_stars[i];
1143 Triangulation &local_tr = m_triangulations[i].tr();
1144 Tr_vertex_handle center_vertex = m_triangulations[i].center_vertex();
1145 int cur_dim_plus_1 = local_tr.current_dimension() + 1;
1147 std::vector<Tr_full_cell_handle> incident_cells;
1148 local_tr.incident_full_cells(center_vertex, std::back_inserter(incident_cells));
1150 typename std::vector<Tr_full_cell_handle>::const_iterator it_c = incident_cells.begin();
1151 typename std::vector<Tr_full_cell_handle>::const_iterator it_c_end = incident_cells.end();
1153 for (; it_c != it_c_end; ++it_c) {
1155 Incident_simplex incident_simplex;
1156 for (
int j = 0; j < cur_dim_plus_1; ++j) {
1157 std::size_t index = (*it_c)->vertex(j)->data();
1158 if (index != i) incident_simplex.insert(index);
1160 GUDHI_CHECK(incident_simplex.size() == cur_dim_plus_1 - 1,
1161 std::logic_error(
"update_star: wrong size of incident simplex"));
1162 star.push_back(incident_simplex);
1168 Tangent_space_basis compute_tangent_space(
const Point &p,
const std::size_t i,
bool normalize_basis =
true,
1169 Orthogonal_space_basis *p_orth_space_basis = NULL) {
1170 unsigned int num_pts_for_pca =
1171 (std::min)(static_cast<unsigned int>(std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1172 static_cast<unsigned int>(m_points.size()));
1175 typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object();
1176 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1178 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 1179 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca,
false);
1180 const Points &points_for_pca = m_points_for_tse;
1182 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca,
false);
1183 const Points &points_for_pca = m_points;
1187 Eigen::MatrixXd mat_points(num_pts_for_pca, m_ambient_dim);
1188 auto nn_it = kns_range.begin();
1189 for (
unsigned int j = 0; j < num_pts_for_pca && nn_it != kns_range.end(); ++j, ++nn_it) {
1190 for (
int i = 0; i < m_ambient_dim; ++i) {
1191 mat_points(j, i) = CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1194 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1195 Eigen::MatrixXd cov = centered.adjoint() * centered;
1196 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1198 Tangent_space_basis tsb(i);
1202 for (
int j = m_ambient_dim - 1; j >= m_ambient_dim - m_intrinsic_dim; --j) {
1203 if (normalize_basis) {
1204 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1205 eig.eigenvectors().col(j).data() + m_ambient_dim);
1206 tsb.push_back(normalize_vector(v, m_k));
1208 tsb.push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1209 eig.eigenvectors().col(j).data() + m_ambient_dim));
1213 if (p_orth_space_basis) {
1214 p_orth_space_basis->set_origin(i);
1215 for (
int j = m_ambient_dim - m_intrinsic_dim - 1; j >= 0; --j) {
1216 if (normalize_basis) {
1217 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1218 eig.eigenvectors().col(j).data() + m_ambient_dim);
1219 p_orth_space_basis->push_back(normalize_vector(v, m_k));
1221 p_orth_space_basis->push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1222 eig.eigenvectors().col(j).data() + m_ambient_dim));
1227 m_are_tangent_spaces_computed[i] =
true;
1237 Tangent_space_basis compute_tangent_space(
const Simplex &s,
bool normalize_basis =
true) {
1238 unsigned int num_pts_for_pca =
1239 (std::min)(static_cast<unsigned int>(std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1240 static_cast<unsigned int>(m_points.size()));
1243 typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object();
1244 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1245 typename K::Squared_length_d sqlen = m_k.squared_length_d_object();
1246 typename K::Scaled_vector_d scaled_vec = m_k.scaled_vector_d_object();
1247 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1248 typename K::Difference_of_vectors_d diff_vec = m_k.difference_of_vectors_d_object();
1251 Eigen::MatrixXd mat_points(s.size() * num_pts_for_pca, m_ambient_dim);
1252 unsigned int current_row = 0;
1254 for (Simplex::const_iterator it_index = s.begin(); it_index != s.end(); ++it_index) {
1255 const Point &p = m_points[*it_index];
1257 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 1258 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca,
false);
1259 const Points &points_for_pca = m_points_for_tse;
1261 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca,
false);
1262 const Points &points_for_pca = m_points;
1265 auto nn_it = kns_range.begin();
1266 for (; current_row < num_pts_for_pca && nn_it != kns_range.end(); ++current_row, ++nn_it) {
1267 for (
int i = 0; i < m_ambient_dim; ++i) {
1268 mat_points(current_row, i) = CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1272 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1273 Eigen::MatrixXd cov = centered.adjoint() * centered;
1274 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1276 Tangent_space_basis tsb;
1280 for (
int j = m_ambient_dim - 1; j >= m_ambient_dim - m_intrinsic_dim; --j) {
1281 if (normalize_basis) {
1282 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1283 eig.eigenvectors().col(j).data() + m_ambient_dim);
1284 tsb.push_back(normalize_vector(v, m_k));
1286 tsb.push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1287 eig.eigenvectors().col(j).data() + m_ambient_dim));
1296 int tangent_basis_dim(std::size_t i)
const {
return m_tangent_spaces[i].dimension(); }
1298 Point compute_perturbed_point(std::size_t pt_idx)
const {
1299 #ifdef GUDHI_TC_PERTURB_POSITION 1300 return m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]);
1302 return m_points[pt_idx];
1306 void compute_perturbed_weighted_point(std::size_t pt_idx, Point &p, FT &w)
const {
1307 #ifdef GUDHI_TC_PERTURB_POSITION 1308 p = m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]);
1310 p = m_points[pt_idx];
1312 w = m_weights[pt_idx];
1315 Weighted_point compute_perturbed_weighted_point(std::size_t pt_idx)
const {
1316 typename K::Construct_weighted_point_d k_constr_wp = m_k.construct_weighted_point_d_object();
1318 Weighted_point wp = k_constr_wp(
1319 #ifdef GUDHI_TC_PERTURB_POSITION
1320 m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]),
1329 Point unproject_point(
const Tr_point &p,
const Tangent_space_basis &tsb,
const Tr_traits &tr_traits)
const {
1330 typename K::Translated_point_d k_transl = m_k.translated_point_d_object();
1331 typename K::Scaled_vector_d k_scaled_vec = m_k.scaled_vector_d_object();
1332 typename Tr_traits::Compute_coordinate_d coord = tr_traits.compute_coordinate_d_object();
1334 Point global_point = compute_perturbed_point(tsb.origin());
1335 for (
int i = 0; i < m_intrinsic_dim; ++i) global_point = k_transl(global_point, k_scaled_vec(tsb[i], coord(p, i)));
1337 return global_point;
1342 Tr_bare_point project_point(
const Point &p,
const Tangent_space_basis &tsb,
const Tr_traits &tr_traits)
const {
1343 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1344 typename K::Difference_of_points_d diff_points = m_k.difference_of_points_d_object();
1346 Vector v = diff_points(p, compute_perturbed_point(tsb.origin()));
1348 std::vector<FT> coords;
1350 coords.reserve(tsb.dimension());
1351 for (std::size_t i = 0; i < m_intrinsic_dim; ++i) {
1353 FT coord = scalar_pdct(v, tsb[i]);
1354 coords.push_back(coord);
1357 return tr_traits.construct_point_d_object()(
static_cast<int>(coords.size()), coords.begin(), coords.end());
1364 Tr_point project_point_and_compute_weight(
const Weighted_point &wp,
const Tangent_space_basis &tsb,
1365 const Tr_traits &tr_traits)
const {
1366 typename K::Point_drop_weight_d k_drop_w = m_k.point_drop_weight_d_object();
1367 typename K::Compute_weight_d k_point_weight = m_k.compute_weight_d_object();
1368 return project_point_and_compute_weight(k_drop_w(wp), k_point_weight(wp), tsb, tr_traits);
1372 Tr_point project_point_and_compute_weight(
const Point &p,
const FT w,
const Tangent_space_basis &tsb,
1373 const Tr_traits &tr_traits)
const {
1374 const int point_dim = m_k.point_dimension_d_object()(p);
1376 typename K::Construct_point_d constr_pt = m_k.construct_point_d_object();
1377 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1378 typename K::Difference_of_points_d diff_points = m_k.difference_of_points_d_object();
1379 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1380 typename K::Construct_cartesian_const_iterator_d ccci = m_k.construct_cartesian_const_iterator_d_object();
1382 Point origin = compute_perturbed_point(tsb.origin());
1383 Vector v = diff_points(p, origin);
1386 bool same_dim = (point_dim == tsb.dimension());
1388 std::vector<FT> coords;
1390 std::vector<FT> p_proj(ccci(origin), ccci(origin, 0));
1391 coords.reserve(tsb.dimension());
1392 for (
int i = 0; i < tsb.dimension(); ++i) {
1394 FT c = scalar_pdct(v, tsb[i]);
1395 coords.push_back(c);
1399 for (
int j = 0; j < point_dim; ++j) p_proj[j] += c * coord(tsb[i], j);
1404 FT sq_dist_to_proj_pt = 0;
1406 Point projected_pt = constr_pt(point_dim, p_proj.begin(), p_proj.end());
1407 sq_dist_to_proj_pt = m_k.squared_distance_d_object()(p, projected_pt);
1410 return tr_traits.construct_weighted_point_d_object()(
1411 tr_traits.construct_point_d_object()(
static_cast<int>(coords.size()), coords.begin(), coords.end()),
1412 w - sq_dist_to_proj_pt);
1417 template <
typename Indexed_po
int_range>
1418 std::vector<Tr_point> project_points_and_compute_weights(
const Indexed_point_range &point_indices,
1419 const Tangent_space_basis &tsb,
1420 const Tr_traits &tr_traits)
const {
1421 std::vector<Tr_point> ret;
1422 for (
typename Indexed_point_range::const_iterator it = point_indices.begin(), it_end = point_indices.end();
1423 it != it_end; ++it) {
1424 ret.push_back(project_point_and_compute_weight(compute_perturbed_weighted_point(*it), tsb, tr_traits));
1431 bool is_simplex_consistent(Tr_full_cell_handle fch,
int cur_dim)
const {
1433 for (
int i = 0; i < cur_dim + 1; ++i) {
1434 std::size_t data = fch->vertex(i)->data();
1437 return is_simplex_consistent(c);
1443 bool is_simplex_consistent(Simplex
const &simplex)
const {
1445 Simplex::const_iterator it_point_idx = simplex.begin();
1448 for (; it_point_idx != simplex.end(); ++it_point_idx) {
1449 std::size_t point_idx = *it_point_idx;
1451 if (point_idx == (std::numeric_limits<std::size_t>::max)())
continue;
1453 Star
const &star = m_stars[point_idx];
1456 Incident_simplex is_to_find = simplex;
1457 is_to_find.erase(point_idx);
1460 if (std::find(star.begin(), star.end(), is_to_find) == star.end())
return false;
1472 template <
typename OutputIterator>
1473 bool is_simplex_consistent(std::size_t center_point,
1474 Incident_simplex
const &s,
1475 OutputIterator points_whose_star_does_not_contain_s,
1476 bool check_also_in_non_maximal_faces =
false)
const {
1477 Simplex full_simplex = s;
1478 full_simplex.insert(center_point);
1481 Incident_simplex::const_iterator it_point_idx = s.begin();
1484 for (; it_point_idx != s.end(); ++it_point_idx) {
1485 std::size_t point_idx = *it_point_idx;
1487 if (point_idx == (std::numeric_limits<std::size_t>::max)())
continue;
1489 Star
const &star = m_stars[point_idx];
1492 Incident_simplex is_to_find = full_simplex;
1493 is_to_find.erase(point_idx);
1495 if (check_also_in_non_maximal_faces) {
1499 for (Star::const_iterator is = star.begin(), is_end = star.end(); !found && is != is_end; ++is) {
1500 if (std::includes(is->begin(), is->end(), is_to_find.begin(), is_to_find.end())) found =
true;
1503 if (!found) *points_whose_star_does_not_contain_s++ = point_idx;
1506 if (std::find(star.begin(), star.end(), is_to_find) == star.end())
1507 *points_whose_star_does_not_contain_s++ = point_idx;
1517 bool is_simplex_in_star(std::size_t p, Incident_simplex
const &s,
bool check_also_in_non_maximal_faces =
true)
const {
1518 Star
const &star = m_stars[p];
1520 if (check_also_in_non_maximal_faces) {
1524 for (Star::const_iterator is = star.begin(), is_end = star.end(); !found && is != is_end; ++is) {
1525 if (std::includes(is->begin(), is->end(), s.begin(), s.end())) found =
true;
1530 return !(std::find(star.begin(), star.end(), s) == star.end());
1534 #ifdef GUDHI_USE_TBB 1536 class Try_to_solve_inconsistencies_in_a_local_triangulation {
1538 double m_max_perturb;
1539 tbb::combinable<std::size_t> &m_num_inconsistencies;
1540 tbb::combinable<std::vector<std::size_t> > &m_updated_points;
1544 Try_to_solve_inconsistencies_in_a_local_triangulation(
Tangential_complex &tc,
double max_perturb,
1545 tbb::combinable<std::size_t> &num_inconsistencies,
1546 tbb::combinable<std::vector<std::size_t> > &updated_points)
1548 m_max_perturb(max_perturb),
1549 m_num_inconsistencies(num_inconsistencies),
1550 m_updated_points(updated_points) {}
1553 Try_to_solve_inconsistencies_in_a_local_triangulation(
1554 const Try_to_solve_inconsistencies_in_a_local_triangulation &tsilt)
1556 m_max_perturb(tsilt.m_max_perturb),
1557 m_num_inconsistencies(tsilt.m_num_inconsistencies),
1558 m_updated_points(tsilt.m_updated_points) {}
1561 void operator()(
const tbb::blocked_range<size_t> &r)
const {
1562 for (
size_t i = r.begin(); i != r.end(); ++i) {
1563 m_num_inconsistencies.local() += m_tc.try_to_solve_inconsistencies_in_a_local_triangulation(
1564 i, m_max_perturb, std::back_inserter(m_updated_points.local()));
1568 #endif // GUDHI_USE_TBB 1570 void perturb(std::size_t point_idx,
double max_perturb) {
1571 const Tr_traits &local_tr_traits = m_triangulations[point_idx].tr().geom_traits();
1572 typename Tr_traits::Compute_coordinate_d coord = local_tr_traits.compute_coordinate_d_object();
1573 typename K::Translated_point_d k_transl = m_k.translated_point_d_object();
1574 typename K::Construct_vector_d k_constr_vec = m_k.construct_vector_d_object();
1575 typename K::Scaled_vector_d k_scaled_vec = m_k.scaled_vector_d_object();
1577 CGAL::Random_points_in_ball_d<Tr_bare_point> tr_point_in_ball_generator(
1578 m_intrinsic_dim, m_random_generator.get_double(0., max_perturb));
1580 Tr_point local_random_transl =
1581 local_tr_traits.construct_weighted_point_d_object()(*tr_point_in_ball_generator++, 0);
1582 Translation_for_perturb global_transl = k_constr_vec(m_ambient_dim);
1583 const Tangent_space_basis &tsb = m_tangent_spaces[point_idx];
1584 for (
int i = 0; i < m_intrinsic_dim; ++i) {
1585 global_transl = k_transl(global_transl, k_scaled_vec(tsb[i], coord(local_random_transl, i)));
1588 #if defined(GUDHI_USE_TBB) 1589 m_p_perturb_mutexes[point_idx].lock();
1590 m_translations[point_idx] = global_transl;
1591 m_p_perturb_mutexes[point_idx].unlock();
1594 m_translations[point_idx] = global_transl;
1599 template <
typename OutputIt>
1600 bool try_to_solve_inconsistencies_in_a_local_triangulation(
1601 std::size_t tr_index,
double max_perturb, OutputIt perturbed_pts_indices = CGAL::Emptyset_iterator()) {
1602 bool is_inconsistent =
false;
1604 Star
const &star = m_stars[tr_index];
1607 Star::const_iterator it_inc_simplex = star.begin();
1608 Star::const_iterator it_inc_simplex_end = star.end();
1609 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1610 const Incident_simplex &incident_simplex = *it_inc_simplex;
1613 if (is_infinite(incident_simplex))
continue;
1615 Simplex c = incident_simplex;
1619 if (!is_simplex_consistent(c)) {
1620 is_inconsistent =
true;
1622 std::size_t idx = tr_index;
1624 perturb(tr_index, max_perturb);
1625 *perturbed_pts_indices++ = idx;
1632 return is_inconsistent;
1637 std::ostream &export_point_set(std::ostream &os,
bool use_perturbed_points =
false,
1638 const char *coord_separator =
" ")
const {
1639 if (use_perturbed_points) {
1640 std::vector<Point> perturbed_points;
1641 perturbed_points.reserve(m_points.size());
1642 for (std::size_t i = 0; i < m_points.size(); ++i) perturbed_points.push_back(compute_perturbed_point(i));
1644 return export_point_set(m_k, perturbed_points, os, coord_separator);
1646 return export_point_set(m_k, m_points, os, coord_separator);
1650 template <
typename ProjectionFunctor = CGAL::Identity<Po
int> >
1651 std::ostream &export_vertices_to_off(std::ostream &os, std::size_t &num_vertices,
bool use_perturbed_points =
false,
1652 ProjectionFunctor
const &point_projection = ProjectionFunctor())
const {
1653 if (m_points.empty()) {
1661 const int N = (m_intrinsic_dim == 1 ? 2 : 1);
1664 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1666 #ifdef GUDHI_TC_EXPORT_ALL_COORDS_IN_OFF 1667 int num_coords = m_ambient_dim;
1669 int num_coords = (std::min)(m_ambient_dim, 3);
1672 #ifdef GUDHI_TC_EXPORT_NORMALS 1673 OS_container::const_iterator it_os = m_orth_spaces.begin();
1675 typename Points::const_iterator it_p = m_points.begin();
1676 typename Points::const_iterator it_p_end = m_points.end();
1678 for (std::size_t i = 0; it_p != it_p_end; ++it_p, ++i) {
1679 Point p = point_projection(use_perturbed_points ? compute_perturbed_point(i) : *it_p);
1680 for (
int ii = 0; ii < N; ++ii) {
1682 for (; j < num_coords; ++j) os << CGAL::to_double(coord(p, j)) <<
" ";
1683 if (j == 2) os <<
"0";
1685 #ifdef GUDHI_TC_EXPORT_NORMALS 1686 for (j = 0; j < num_coords; ++j) os <<
" " << CGAL::to_double(coord(*it_os->begin(), j));
1690 #ifdef GUDHI_TC_EXPORT_NORMALS 1695 num_vertices = N * m_points.size();
1699 std::ostream &export_simplices_to_off(std::ostream &os, std::size_t &num_OFF_simplices,
1700 bool color_inconsistencies =
false,
1701 Simplex_set
const *p_simpl_to_color_in_red = NULL,
1702 Simplex_set
const *p_simpl_to_color_in_green = NULL,
1703 Simplex_set
const *p_simpl_to_color_in_blue = NULL)
const {
1706 num_OFF_simplices = 0;
1707 std::size_t num_maximal_simplices = 0;
1708 std::size_t num_inconsistent_maximal_simplices = 0;
1709 std::size_t num_inconsistent_stars = 0;
1710 typename Tr_container::const_iterator it_tr = m_triangulations.begin();
1711 typename Tr_container::const_iterator it_tr_end = m_triangulations.end();
1713 for (std::size_t idx = 0; it_tr != it_tr_end; ++it_tr, ++idx) {
1714 bool is_star_inconsistent =
false;
1716 Triangulation
const &tr = it_tr->tr();
1718 if (tr.current_dimension() < m_intrinsic_dim)
continue;
1721 std::stringstream color;
1723 color << 128 <<
" " << 128 <<
" " << 128;
1726 typedef std::vector<std::pair<Simplex, int> > Star_using_triangles;
1727 Star_using_triangles star_using_triangles;
1730 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
1731 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
1732 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1733 Simplex c = *it_inc_simplex;
1735 std::size_t num_vertices = c.size();
1736 ++num_maximal_simplices;
1738 int color_simplex = -1;
1739 if (color_inconsistencies && !is_simplex_consistent(c)) {
1740 ++num_inconsistent_maximal_simplices;
1742 is_star_inconsistent =
true;
1744 if (p_simpl_to_color_in_red && std::find(p_simpl_to_color_in_red->begin(), p_simpl_to_color_in_red->end(),
1745 c) != p_simpl_to_color_in_red->end()) {
1747 }
else if (p_simpl_to_color_in_green &&
1748 std::find(p_simpl_to_color_in_green->begin(), p_simpl_to_color_in_green->end(), c) !=
1749 p_simpl_to_color_in_green->end()) {
1751 }
else if (p_simpl_to_color_in_blue &&
1752 std::find(p_simpl_to_color_in_blue->begin(), p_simpl_to_color_in_blue->end(), c) !=
1753 p_simpl_to_color_in_blue->end()) {
1762 if (m_intrinsic_dim == 1) {
1764 Simplex::iterator it = c.begin();
1765 for (; it != c.end(); ++it) tmp_c.insert(*it * 2);
1766 if (num_vertices == 2) tmp_c.insert(*tmp_c.rbegin() + 1);
1771 if (num_vertices <= 3) {
1772 star_using_triangles.push_back(std::make_pair(c, color_simplex));
1775 std::vector<bool> booleans(num_vertices,
false);
1776 std::fill(booleans.begin() + num_vertices - 3, booleans.end(),
true);
1779 Simplex::iterator it = c.begin();
1780 for (
int i = 0; it != c.end(); ++i, ++it) {
1781 if (booleans[i]) triangle.insert(*it);
1783 star_using_triangles.push_back(std::make_pair(triangle, color_simplex));
1784 }
while (std::next_permutation(booleans.begin(), booleans.end()));
1789 Star_using_triangles::const_iterator it_simplex = star_using_triangles.begin();
1790 Star_using_triangles::const_iterator it_simplex_end = star_using_triangles.end();
1791 for (; it_simplex != it_simplex_end; ++it_simplex) {
1792 const Simplex &c = it_simplex->first;
1795 if (is_infinite(c))
continue;
1797 int color_simplex = it_simplex->second;
1799 std::stringstream sstr_c;
1801 Simplex::const_iterator it_point_idx = c.begin();
1802 for (; it_point_idx != c.end(); ++it_point_idx) {
1803 sstr_c << *it_point_idx <<
" ";
1806 os << 3 <<
" " << sstr_c.str();
1807 if (color_inconsistencies || p_simpl_to_color_in_red || p_simpl_to_color_in_green || p_simpl_to_color_in_blue) {
1808 switch (color_simplex) {
1822 os <<
" " << color.str();
1826 ++num_OFF_simplices;
1829 if (is_star_inconsistent) ++num_inconsistent_stars;
1833 std::cerr <<
"\n==========================================================\n" 1834 <<
"Export from list of stars to OFF:\n" 1835 <<
" * Number of vertices: " << m_points.size() <<
"\n" 1836 <<
" * Total number of maximal simplices: " << num_maximal_simplices <<
"\n";
1837 if (color_inconsistencies) {
1838 std::cerr <<
" * Number of inconsistent stars: " << num_inconsistent_stars <<
" (" 1839 << (m_points.size() > 0 ? 100. * num_inconsistent_stars / m_points.size() : 0.) <<
"%)\n" 1840 <<
" * Number of inconsistent maximal simplices: " << num_inconsistent_maximal_simplices <<
" (" 1841 << (num_maximal_simplices > 0 ? 100. * num_inconsistent_maximal_simplices / num_maximal_simplices : 0.)
1844 std::cerr <<
"==========================================================\n";
1851 std::ostream &export_simplices_to_off(
const Simplicial_complex &complex, std::ostream &os,
1852 std::size_t &num_OFF_simplices,
1853 Simplex_set
const *p_simpl_to_color_in_red = NULL,
1854 Simplex_set
const *p_simpl_to_color_in_green = NULL,
1855 Simplex_set
const *p_simpl_to_color_in_blue = NULL)
const {
1856 typedef Simplicial_complex::Simplex Simplex;
1857 typedef Simplicial_complex::Simplex_set Simplex_set;
1861 num_OFF_simplices = 0;
1862 std::size_t num_maximal_simplices = 0;
1864 typename Simplex_set::const_iterator it_s = complex.simplex_range().begin();
1865 typename Simplex_set::const_iterator it_s_end = complex.simplex_range().end();
1867 for (; it_s != it_s_end; ++it_s) {
1869 ++num_maximal_simplices;
1871 int color_simplex = -1;
1872 if (p_simpl_to_color_in_red && std::find(p_simpl_to_color_in_red->begin(), p_simpl_to_color_in_red->end(), c) !=
1873 p_simpl_to_color_in_red->end()) {
1875 }
else if (p_simpl_to_color_in_green &&
1876 std::find(p_simpl_to_color_in_green->begin(), p_simpl_to_color_in_green->end(), c) !=
1877 p_simpl_to_color_in_green->end()) {
1879 }
else if (p_simpl_to_color_in_blue &&
1880 std::find(p_simpl_to_color_in_blue->begin(), p_simpl_to_color_in_blue->end(), c) !=
1881 p_simpl_to_color_in_blue->end()) {
1886 typedef std::vector<Simplex> Triangles;
1887 Triangles triangles;
1889 int num_vertices =
static_cast<int>(c.size());
1891 if (num_vertices < m_intrinsic_dim + 1)
continue;
1897 if (m_intrinsic_dim == 1) {
1899 Simplex::iterator it = c.begin();
1900 for (; it != c.end(); ++it) tmp_c.insert(*it * 2);
1901 if (num_vertices == 2) tmp_c.insert(*tmp_c.rbegin() + 1);
1906 if (num_vertices <= 3) {
1907 triangles.push_back(c);
1910 std::vector<bool> booleans(num_vertices,
false);
1911 std::fill(booleans.begin() + num_vertices - 3, booleans.end(),
true);
1914 Simplex::iterator it = c.begin();
1915 for (
int i = 0; it != c.end(); ++i, ++it) {
1916 if (booleans[i]) triangle.insert(*it);
1918 triangles.push_back(triangle);
1919 }
while (std::next_permutation(booleans.begin(), booleans.end()));
1923 Triangles::const_iterator it_tri = triangles.begin();
1924 Triangles::const_iterator it_tri_end = triangles.end();
1925 for (; it_tri != it_tri_end; ++it_tri) {
1927 if (is_infinite(*it_tri))
continue;
1930 Simplex::const_iterator it_point_idx = it_tri->begin();
1931 for (; it_point_idx != it_tri->end(); ++it_point_idx) {
1932 os << *it_point_idx <<
" ";
1935 if (p_simpl_to_color_in_red || p_simpl_to_color_in_green || p_simpl_to_color_in_blue) {
1936 switch (color_simplex) {
1950 os <<
" 128 128 128";
1955 ++num_OFF_simplices;
1961 std::cerr <<
"\n==========================================================\n" 1962 <<
"Export from complex to OFF:\n" 1963 <<
" * Number of vertices: " << m_points.size() <<
"\n" 1964 <<
" * Total number of maximal simplices: " << num_maximal_simplices <<
"\n" 1965 <<
"==========================================================\n";
1973 const int m_intrinsic_dim;
1974 const int m_ambient_dim;
1978 #ifdef GUDHI_TC_PERTURB_POSITION 1979 Translations_for_perturb m_translations;
1980 #if defined(GUDHI_USE_TBB) 1981 Mutex_for_perturb *m_p_perturb_mutexes;
1985 Points_ds m_points_ds;
1986 double m_last_max_perturb;
1987 std::vector<bool> m_are_tangent_spaces_computed;
1988 TS_container m_tangent_spaces;
1989 #ifdef GUDHI_TC_EXPORT_NORMALS 1990 OS_container m_orth_spaces;
1992 Tr_container m_triangulations;
1994 Stars_container m_stars;
1995 std::vector<FT> m_squared_star_spheres_radii_incl_margin;
1997 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM 1998 Points m_points_for_tse;
1999 Points_ds m_points_ds_for_tse;
2002 mutable CGAL::Random m_random_generator;
2008 #endif // TANGENTIAL_COMPLEX_H_
std::size_t num_inconsistent_simplices
Number of inconsistent simplices.
Definition: Tangential_complex.h:538
Type returned by Tangential_complex::fix_inconsistencies_using_perturbation.
Definition: Tangential_complex.h:372
Point get_point(std::size_t vertex) const
Returns the point corresponding to the vertex given as parameter.
Definition: Tangential_complex.h:290
Tangential complex data structure.
Definition: Tangential_complex.h:124
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:382
Definition: SimplicialComplexForAlpha.h:26
~Tangential_complex()
Destructor.
Definition: Tangential_complex.h:271
std::size_t num_simplices
Total number of simplices in stars (including duplicates that appear in several stars) ...
Definition: Tangential_complex.h:536
std::size_t best_num_inconsistent_stars
best number of inconsistent stars during the process
Definition: Tangential_complex.h:380
std::size_t num_inconsistent_stars
Number of stars containing at least one inconsistent simplex.
Definition: Tangential_complex.h:540
int create_complex(Simplex_tree_ &tree, bool export_inconsistent_simplices=true) const
Exports the complex into a Simplex_tree.
Definition: Tangential_complex.h:605
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:390
std::size_t initial_num_inconsistent_stars
initial number of inconsistent stars
Definition: Tangential_complex.h:378
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:239
Num_inconsistencies number_of_inconsistent_simplices(bool verbose=false) const
Definition: Tangential_complex.h:546
bool success
true if all inconsistencies could be removed, false if the time limit has been reached before ...
Definition: Tangential_complex.h:374
int intrinsic_dimension() const
Returns the intrinsic dimension of the manifold.
Definition: Tangential_complex.h:278
Definition: Intro_spatial_searching.h:28
void compute_tangential_complex()
Computes the tangential complex.
Definition: Tangential_complex.h:325
int ambient_dimension() const
Returns the ambient dimension.
Definition: Tangential_complex.h:281
Type returned by Tangential_complex::number_of_inconsistent_simplices.
Definition: Tangential_complex.h:534
std::size_t number_of_vertices() const
Returns the number of vertices.
Definition: Tangential_complex.h:301
unsigned int num_steps
number of steps performed
Definition: Tangential_complex.h:376
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:297