Tangential_complex.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Clement Jamin
4 *
5 * Copyright (C) 2016 Inria
6 *
7 * Modification(s):
8 * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3
9 * - YYYY/MM Author: Description of the modification
10 */
11
12#ifndef TANGENTIAL_COMPLEX_H_
13#define TANGENTIAL_COMPLEX_H_
14
15#include <gudhi/Tangential_complex/config.h>
16#include <gudhi/Tangential_complex/Simplicial_complex.h>
17#include <gudhi/Tangential_complex/utilities.h>
18#include <gudhi/Kd_tree_search.h>
19#include <gudhi/console_color.h>
20#include <gudhi/Clock.h>
21#include <gudhi/Simplex_tree.h>
22#include <gudhi/Debug_utils.h>
23
24#include <CGAL/Default.h>
25#include <CGAL/Dimension.h>
26#include <CGAL/function_objects.h> // for CGAL::Identity
27#include <CGAL/Epick_d.h>
28#include <CGAL/Regular_triangulation_traits_adapter.h>
29#include <CGAL/Regular_triangulation.h>
30#include <CGAL/Delaunay_triangulation.h>
31#include <CGAL/Combination_enumerator.h>
32#include <CGAL/point_generators_d.h>
33#include <CGAL/version.h> // for CGAL_VERSION_NR
34
35#include <Eigen/Core>
36#include <Eigen/Eigen>
37#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST
38
39#include <boost/optional.hpp>
40#include <boost/iterator/transform_iterator.hpp>
41#include <boost/range/adaptor/transformed.hpp>
42#include <boost/range/counting_range.hpp>
43#include <boost/math/special_functions/factorials.hpp>
44#include <boost/container/flat_set.hpp>
45
46#include <tuple>
47#include <vector>
48#include <set>
49#include <utility>
50#include <sstream>
51#include <iostream>
52#include <limits>
53#include <algorithm>
54#include <functional>
55#include <iterator>
56#include <cmath> // for std::sqrt
57#include <string>
58#include <cstddef> // for std::size_t
59
60#ifdef GUDHI_USE_TBB
61#include <tbb/parallel_for.h>
62#include <tbb/combinable.h>
63#include <mutex>
64#endif
65
66// #define GUDHI_TC_EXPORT_NORMALS // Only for 3D surfaces (k=2, d=3)
67
68// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
69#if CGAL_VERSION_NR < 1041101000
70# error Tangential_complex is only available for CGAL >= 4.11
71#endif
72
73#if !EIGEN_VERSION_AT_LEAST(3,1,0)
74# error Tangential_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
75#endif
76
77namespace sps = Gudhi::spatial_searching;
78
79namespace Gudhi {
80
81namespace tangential_complex {
82
83using namespace internal;
84
85class Vertex_data {
86 public:
87 Vertex_data(std::size_t data = (std::numeric_limits<std::size_t>::max)()) : m_data(data) {}
88
89 operator std::size_t() { return m_data; }
90
91 operator std::size_t() const { return m_data; }
92
93 private:
94 std::size_t m_data;
95};
96
122template <typename Kernel_, // ambiant kernel
123 typename DimensionTag, // intrinsic dimension
124 typename Concurrency_tag = CGAL::Parallel_tag, typename Triangulation_ = CGAL::Default>
126 typedef Kernel_ K;
127 typedef typename K::FT FT;
128 typedef typename K::Point_d Point;
129 typedef typename K::Weighted_point_d Weighted_point;
130 typedef typename K::Vector_d Vector;
131
132 typedef typename CGAL::Default::Get<
133 Triangulation_,
134 CGAL::Regular_triangulation<
135 CGAL::Epick_d<DimensionTag>,
136 CGAL::Triangulation_data_structure<
137 typename CGAL::Epick_d<DimensionTag>::Dimension,
138 CGAL::Triangulation_vertex<CGAL::Regular_triangulation_traits_adapter<CGAL::Epick_d<DimensionTag> >,
139 Vertex_data>,
140 CGAL::Triangulation_full_cell<
141 CGAL::Regular_triangulation_traits_adapter<CGAL::Epick_d<DimensionTag> > > > > >::type Triangulation;
142 typedef typename Triangulation::Geom_traits Tr_traits;
143 typedef typename Triangulation::Weighted_point Tr_point;
144 typedef typename Tr_traits::Base::Point_d Tr_bare_point;
145 typedef typename Triangulation::Vertex_handle Tr_vertex_handle;
146 typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle;
147 typedef typename Tr_traits::Vector_d Tr_vector;
148
149#if defined(GUDHI_USE_TBB)
150 typedef std::mutex Mutex_for_perturb;
151 typedef Vector Translation_for_perturb;
152 typedef std::vector<Atomic_wrapper<FT> > Weights;
153#else
154 typedef Vector Translation_for_perturb;
155 typedef std::vector<FT> Weights;
156#endif
157 typedef std::vector<Translation_for_perturb> Translations_for_perturb;
158
159 // Store a local triangulation and a handle to its center vertex
160
161 struct Tr_and_VH {
162 public:
163 Tr_and_VH() : m_tr(NULL) {}
164
165 Tr_and_VH(int dim) : m_tr(new Triangulation(dim)) {}
166
167 ~Tr_and_VH() { destroy_triangulation(); }
168
169 Triangulation &construct_triangulation(int dim) {
170 delete m_tr;
171 m_tr = new Triangulation(dim);
172 return tr();
173 }
174
175 void destroy_triangulation() {
176 delete m_tr;
177 m_tr = NULL;
178 }
179
180 Triangulation &tr() { return *m_tr; }
181
182 Triangulation const &tr() const { return *m_tr; }
183
184 Tr_vertex_handle const &center_vertex() const { return m_center_vertex; }
185
186 Tr_vertex_handle &center_vertex() { return m_center_vertex; }
187
188 private:
189 Triangulation *m_tr;
190 Tr_vertex_handle m_center_vertex;
191 };
192
193 public:
194 typedef Basis<K> Tangent_space_basis;
195 typedef Basis<K> Orthogonal_space_basis;
196 typedef std::vector<Tangent_space_basis> TS_container;
197 typedef std::vector<Orthogonal_space_basis> OS_container;
198
199 typedef std::vector<Point> Points;
200
201 typedef boost::container::flat_set<std::size_t> Simplex;
202 typedef std::set<Simplex> Simplex_set;
203
204 private:
206 typedef typename Points_ds::KNS_range KNS_range;
207 typedef typename Points_ds::INS_range INS_range;
208
209 typedef std::vector<Tr_and_VH> Tr_container;
210 typedef std::vector<Vector> Vectors;
211
212 // An Incident_simplex is the list of the vertex indices
213 // except the center vertex
214 typedef boost::container::flat_set<std::size_t> Incident_simplex;
215 typedef std::vector<Incident_simplex> Star;
216 typedef std::vector<Star> Stars_container;
217
218 // For transform_iterator
219
220 static const Tr_point &vertex_handle_to_point(Tr_vertex_handle vh) { return vh->point(); }
221
222 template <typename P, typename VH>
223 static const P &vertex_handle_to_point(VH vh) {
224 return vh->point();
225 }
226
227 public:
228 typedef internal::Simplicial_complex Simplicial_complex;
229
239 template <typename Point_range>
241#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
242 InputIterator first_for_tse, InputIterator last_for_tse,
243#endif
244 const K &k = K())
245 : m_k(k),
246 m_intrinsic_dim(intrinsic_dimension),
247 m_ambient_dim(points.empty() ? 0 : k.point_dimension_d_object()(*points.begin())),
248 m_points(points.begin(), points.end()),
249 m_weights(m_points.size(), FT(0))
250#if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION)
251 ,
252 m_p_perturb_mutexes(NULL)
253#endif
254 ,
255 m_points_ds(m_points),
256 m_last_max_perturb(0.),
257 m_are_tangent_spaces_computed(m_points.size(), false),
258 m_tangent_spaces(m_points.size(), Tangent_space_basis())
259#ifdef GUDHI_TC_EXPORT_NORMALS
260 ,
261 m_orth_spaces(m_points.size(), Orthogonal_space_basis())
262#endif
263#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
264 ,
265 m_points_for_tse(first_for_tse, last_for_tse),
266 m_points_ds_for_tse(m_points_for_tse)
267#endif
268 {
269 }
270
273#if defined(GUDHI_USE_TBB) && defined(GUDHI_TC_PERTURB_POSITION)
274 delete[] m_p_perturb_mutexes;
275#endif
276 }
277
279 int intrinsic_dimension() const { return m_intrinsic_dim; }
280
282 int ambient_dimension() const { return m_ambient_dim; }
283
284 Points const &points() const { return m_points; }
285
291 Point get_point(std::size_t vertex) const { return m_points[vertex]; }
292
298 Point get_perturbed_point(std::size_t vertex) const { return compute_perturbed_point(vertex); }
299
301
302 std::size_t number_of_vertices() const { return m_points.size(); }
303
304 void set_weights(const Weights &weights) { m_weights = weights; }
305
306 void set_tangent_planes(const TS_container &tangent_spaces
307#ifdef GUDHI_TC_EXPORT_NORMALS
308 ,
309 const OS_container &orthogonal_spaces
310#endif
311 ) {
312#ifdef GUDHI_TC_EXPORT_NORMALS
313 GUDHI_CHECK(m_points.size() == tangent_spaces.size() && m_points.size() == orthogonal_spaces.size(),
314 std::logic_error("Wrong sizes"));
315#else
316 GUDHI_CHECK(m_points.size() == tangent_spaces.size(), std::logic_error("Wrong sizes"));
317#endif
318 m_tangent_spaces = tangent_spaces;
319#ifdef GUDHI_TC_EXPORT_NORMALS
320 m_orth_spaces = orthogonal_spaces;
321#endif
322 for (std::size_t i = 0; i < m_points.size(); ++i) m_are_tangent_spaces_computed[i] = true;
323 }
324
331#ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS
332 std::cerr << red << "WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. "
333 << "Computation might be slower than usual.\n"
334 << white;
335#endif
336
337#if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB)
338 Gudhi::Clock t;
339#endif
340
341 // We need to do that because we don't want the container to copy the
342 // already-computed triangulations (while resizing) since it would
343 // invalidate the vertex handles stored beside the triangulations
344 m_triangulations.resize(m_points.size());
345 m_stars.resize(m_points.size());
346 m_squared_star_spheres_radii_incl_margin.resize(m_points.size(), FT(-1));
347#ifdef GUDHI_TC_PERTURB_POSITION
348 if (m_points.empty())
349 m_translations.clear();
350 else
351 m_translations.resize(m_points.size(), m_k.construct_vector_d_object()(m_ambient_dim));
352#if defined(GUDHI_USE_TBB)
353 delete[] m_p_perturb_mutexes;
354 m_p_perturb_mutexes = new Mutex_for_perturb[m_points.size()];
355#endif
356#endif
357
358#ifdef GUDHI_USE_TBB
359 // Parallel
360 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
361 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*this));
362 } else {
363#endif // GUDHI_USE_TBB
364 // Sequential
365 for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
366#ifdef GUDHI_USE_TBB
367 }
368#endif // GUDHI_USE_TBB
369
370#if defined(GUDHI_TC_PROFILING) && defined(GUDHI_USE_TBB)
371 t.end();
372 std::cerr << "Tangential complex computed in " << t.num_seconds() << " seconds.\n";
373#endif
374 }
375
379 bool success = false;
381 unsigned int num_steps = 0;
388 };
389
395 Fix_inconsistencies_info fix_inconsistencies_using_perturbation(double max_perturb, double time_limit = -1.) {
397
398 if (time_limit == 0.) return info;
399
400 Gudhi::Clock t;
401
402#ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES
403 std::tuple<std::size_t, std::size_t, std::size_t> stats_before = number_of_inconsistent_simplices(false);
404
405 if (std::get<1>(stats_before) == 0) {
406#ifdef DEBUG_TRACES
407 std::cerr << "Nothing to fix.\n";
408#endif
409 info.success = false;
410 return info;
411 }
412#endif // GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES
413
414 m_last_max_perturb = max_perturb;
415
416 bool done = false;
417 info.best_num_inconsistent_stars = m_triangulations.size();
418 info.num_steps = 0;
419 while (!done) {
420#ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES
421 std::cerr << "\nBefore fix step:\n"
422 << " * Total number of simplices in stars (incl. duplicates): " << std::get<0>(stats_before) << "\n"
423 << " * Num inconsistent simplices in stars (incl. duplicates): " << red << std::get<1>(stats_before)
424 << white << " (" << 100. * std::get<1>(stats_before) / std::get<0>(stats_before) << "%)\n"
425 << " * Number of stars containing inconsistent simplices: " << red << std::get<2>(stats_before)
426 << white << " (" << 100. * std::get<2>(stats_before) / m_points.size() << "%)\n";
427#endif
428
429#if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING)
430 std::cerr << yellow << "\nAttempt to fix inconsistencies using perturbations - step #" << info.num_steps + 1
431 << "... " << white;
432#endif
433
434 std::size_t num_inconsistent_stars = 0;
435 std::vector<std::size_t> updated_points;
436
437#ifdef GUDHI_TC_PROFILING
438 Gudhi::Clock t_fix_step;
439#endif
440
441 // Parallel
442#if defined(GUDHI_USE_TBB)
443 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
444 tbb::combinable<std::size_t> num_inconsistencies;
445 tbb::combinable<std::vector<std::size_t> > tls_updated_points;
446 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_triangulations.size()),
447 Try_to_solve_inconsistencies_in_a_local_triangulation(*this, max_perturb, num_inconsistencies,
448 tls_updated_points));
449 num_inconsistent_stars = num_inconsistencies.combine(std::plus<std::size_t>());
450 updated_points =
451 tls_updated_points.combine([](std::vector<std::size_t> const &x, std::vector<std::size_t> const &y) {
452 std::vector<std::size_t> res;
453 res.reserve(x.size() + y.size());
454 res.insert(res.end(), x.begin(), x.end());
455 res.insert(res.end(), y.begin(), y.end());
456 return res;
457 });
458 } else {
459#endif // GUDHI_USE_TBB
460 // Sequential
461 for (std::size_t i = 0; i < m_triangulations.size(); ++i) {
462 num_inconsistent_stars +=
463 try_to_solve_inconsistencies_in_a_local_triangulation(i, max_perturb, std::back_inserter(updated_points));
464 }
465#if defined(GUDHI_USE_TBB)
466 }
467#endif // GUDHI_USE_TBB
468
469#ifdef GUDHI_TC_PROFILING
470 t_fix_step.end();
471#endif
472
473#if defined(GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES) || defined(DEBUG_TRACES)
474 std::cerr << "\nEncountered during fix:\n"
475 << " * Num stars containing inconsistent simplices: " << red << num_inconsistent_stars << white << " ("
476 << 100. * num_inconsistent_stars / m_points.size() << "%)\n";
477#endif
478
479#ifdef GUDHI_TC_PROFILING
480 std::cerr << yellow << "done in " << t_fix_step.num_seconds() << " seconds.\n" << white;
481#elif defined(DEBUG_TRACES)
482 std::cerr << yellow << "done.\n" << white;
483#endif
484
485 if (num_inconsistent_stars > 0) refresh_tangential_complex(updated_points);
486
487#ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS
488 // Confirm that all stars were actually refreshed
489 std::size_t num_inc_1 = std::get<1>(number_of_inconsistent_simplices(false));
490 refresh_tangential_complex();
491 std::size_t num_inc_2 = std::get<1>(number_of_inconsistent_simplices(false));
492 if (num_inc_1 != num_inc_2)
493 std::cerr << red << "REFRESHMENT CHECK: FAILED. (" << num_inc_1 << " vs " << num_inc_2 << ")\n" << white;
494 else
495 std::cerr << green << "REFRESHMENT CHECK: PASSED.\n" << white;
496#endif
497
498#ifdef GUDHI_TC_SHOW_DETAILED_STATS_FOR_INCONSISTENCIES
499 std::tuple<std::size_t, std::size_t, std::size_t> stats_after = number_of_inconsistent_simplices(false);
500
501 std::cerr << "\nAfter fix:\n"
502 << " * Total number of simplices in stars (incl. duplicates): " << std::get<0>(stats_after) << "\n"
503 << " * Num inconsistent simplices in stars (incl. duplicates): " << red << std::get<1>(stats_after)
504 << white << " (" << 100. * std::get<1>(stats_after) / std::get<0>(stats_after) << "%)\n"
505 << " * Number of stars containing inconsistent simplices: " << red << std::get<2>(stats_after) << white
506 << " (" << 100. * std::get<2>(stats_after) / m_points.size() << "%)\n";
507
508 stats_before = stats_after;
509#endif
510
511 if (info.num_steps == 0) info.initial_num_inconsistent_stars = num_inconsistent_stars;
512
513 if (num_inconsistent_stars < info.best_num_inconsistent_stars)
514 info.best_num_inconsistent_stars = num_inconsistent_stars;
515
516 info.final_num_inconsistent_stars = num_inconsistent_stars;
517
518 done = (num_inconsistent_stars == 0);
519 if (!done) {
520 ++info.num_steps;
521 if (time_limit > 0. && t.num_seconds() > time_limit) {
522#ifdef DEBUG_TRACES
523 std::cerr << red << "Time limit reached.\n" << white;
524#endif
525 info.success = false;
526 return info;
527 }
528 }
529 }
530
531#ifdef DEBUG_TRACES
532 std::cerr << green << "Fixed!\n" << white;
533#endif
534 info.success = true;
535 return info;
536 }
537
541 std::size_t num_simplices = 0;
545 std::size_t num_inconsistent_stars = 0;
546 };
547
550
552#ifdef DEBUG_TRACES
553 bool verbose = true
554#else
555 bool verbose = false
556#endif
557 ) const {
559
560 // For each triangulation
561 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
562 bool is_star_inconsistent = false;
563
564 // For each cell
565 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
566 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
567 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
568 // Don't check infinite cells
569 if (is_infinite(*it_inc_simplex)) continue;
570
571 Simplex c = *it_inc_simplex;
572 c.insert(idx); // Add the missing index
573
574 if (!is_simplex_consistent(c)) {
576 is_star_inconsistent = true;
577 }
578
579 ++stats.num_simplices;
580 }
581 stats.num_inconsistent_stars += is_star_inconsistent;
582 }
583
584 if (verbose) {
585 std::cerr << "\n==========================================================\n"
586 << "Inconsistencies:\n"
587 << " * Total number of simplices in stars (incl. duplicates): " << stats.num_simplices << "\n"
588 << " * Number of inconsistent simplices in stars (incl. duplicates): "
589 << stats.num_inconsistent_simplices << " ("
590 << 100. * stats.num_inconsistent_simplices / stats.num_simplices << "%)\n"
591 << " * Number of stars containing inconsistent simplices: " << stats.num_inconsistent_stars << " ("
592 << 100. * stats.num_inconsistent_stars / m_points.size() << "%)\n"
593 << "==========================================================\n";
594 }
595
596 return stats;
597 }
598
609 template <typename Simplex_tree_>
610 int create_complex(Simplex_tree_ &tree,
611 bool export_inconsistent_simplices = true
613 ,
614 bool export_infinite_simplices = false, Simplex_set *p_inconsistent_simplices = NULL
616 ) const {
617#if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING)
618 std::cerr << yellow << "\nExporting the TC as a Simplex_tree... " << white;
619#endif
620#ifdef GUDHI_TC_PROFILING
621 Gudhi::Clock t;
622#endif
623
624 int max_dim = -1;
625
626 // For each triangulation
627 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
628 // For each cell of the star
629 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
630 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
631 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
632 Simplex c = *it_inc_simplex;
633
634 // Don't export infinite cells
635 if (!export_infinite_simplices && is_infinite(c)) continue;
636
637 if (static_cast<int>(c.size()) > max_dim) max_dim = static_cast<int>(c.size());
638 // Add the missing center vertex
639 c.insert(idx);
640
641 if (!export_inconsistent_simplices && !is_simplex_consistent(c)) continue;
642
643 // Try to insert the simplex
644 bool inserted = tree.insert_simplex_and_subfaces(c).second;
645
646 // Inconsistent?
647 if (p_inconsistent_simplices && inserted && !is_simplex_consistent(c)) {
648 p_inconsistent_simplices->insert(c);
649 }
650 }
651 }
652
653#ifdef GUDHI_TC_PROFILING
654 t.end();
655 std::cerr << yellow << "done in " << t.num_seconds() << " seconds.\n" << white;
656#elif defined(DEBUG_TRACES)
657 std::cerr << yellow << "done.\n" << white;
658#endif
659
660 return max_dim;
661 }
662
663 // First clears the complex then exports the TC into it
664 // Returns the max dimension of the simplices
665 // check_lower_and_higher_dim_simplices : 0 (false), 1 (true), 2 (auto)
666 // If the check is enabled, the function:
667 // - won't insert the simplex if it is already in a higher dim simplex
668 // - will erase any lower-dim simplices that are faces of the new simplex
669 // "auto" (= 2) will enable the check as a soon as it encounters a
670 // simplex whose dimension is different from the previous ones.
671 // N.B.: The check is quite expensive.
672
673 int create_complex(Simplicial_complex &complex, bool export_inconsistent_simplices = true,
674 bool export_infinite_simplices = false, int check_lower_and_higher_dim_simplices = 2,
675 Simplex_set *p_inconsistent_simplices = NULL) const {
676#if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING)
677 std::cerr << yellow << "\nExporting the TC as a Simplicial_complex... " << white;
678#endif
679#ifdef GUDHI_TC_PROFILING
680 Gudhi::Clock t;
681#endif
682
683 int max_dim = -1;
684 complex.clear();
685
686 // For each triangulation
687 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
688 // For each cell of the star
689 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
690 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
691 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
692 Simplex c = *it_inc_simplex;
693
694 // Don't export infinite cells
695 if (!export_infinite_simplices && is_infinite(c)) continue;
696
697 if (static_cast<int>(c.size()) > max_dim) max_dim = static_cast<int>(c.size());
698 // Add the missing center vertex
699 c.insert(idx);
700
701 if (!export_inconsistent_simplices && !is_simplex_consistent(c)) continue;
702
703 // Unusual simplex dim?
704 if (check_lower_and_higher_dim_simplices == 2 && max_dim != -1 && static_cast<int>(c.size()) != max_dim) {
705 // Let's activate the check
706 std::cerr << red
707 << "Info: check_lower_and_higher_dim_simplices ACTIVATED. "
708 "Export might be take some time...\n"
709 << white;
710 check_lower_and_higher_dim_simplices = 1;
711 }
712
713 // Try to insert the simplex
714 bool added = complex.add_simplex(c, check_lower_and_higher_dim_simplices == 1);
715
716 // Inconsistent?
717 if (p_inconsistent_simplices && added && !is_simplex_consistent(c)) {
718 p_inconsistent_simplices->insert(c);
719 }
720 }
721 }
722
723#ifdef GUDHI_TC_PROFILING
724 t.end();
725 std::cerr << yellow << "done in " << t.num_seconds() << " seconds.\n" << white;
726#elif defined(DEBUG_TRACES)
727 std::cerr << yellow << "done.\n" << white;
728#endif
729
730 return max_dim;
731 }
732
733 template <typename ProjectionFunctor = CGAL::Identity<Point> >
734 std::ostream &export_to_off(const Simplicial_complex &complex, std::ostream &os,
735 Simplex_set const *p_simpl_to_color_in_red = NULL,
736 Simplex_set const *p_simpl_to_color_in_green = NULL,
737 Simplex_set const *p_simpl_to_color_in_blue = NULL,
738 ProjectionFunctor const &point_projection = ProjectionFunctor()) const {
739 return export_to_off(os, false, p_simpl_to_color_in_red, p_simpl_to_color_in_green, p_simpl_to_color_in_blue,
740 &complex, point_projection);
741 }
742
743 template <typename ProjectionFunctor = CGAL::Identity<Point> >
744 std::ostream &export_to_off(std::ostream &os, bool color_inconsistencies = false,
745 Simplex_set const *p_simpl_to_color_in_red = NULL,
746 Simplex_set const *p_simpl_to_color_in_green = NULL,
747 Simplex_set const *p_simpl_to_color_in_blue = NULL,
748 const Simplicial_complex *p_complex = NULL,
749 ProjectionFunctor const &point_projection = ProjectionFunctor()) const {
750 if (m_points.empty()) return os;
751
752 if (m_ambient_dim < 2) {
753 std::cerr << "Error: export_to_off => ambient dimension should be >= 2.\n";
754 os << "Error: export_to_off => ambient dimension should be >= 2.\n";
755 return os;
756 }
757 if (m_ambient_dim > 3) {
758 std::cerr << "Warning: export_to_off => ambient dimension should be "
759 "<= 3. Only the first 3 coordinates will be exported.\n";
760 }
761
762 if (m_intrinsic_dim < 1 || m_intrinsic_dim > 3) {
763 std::cerr << "Error: export_to_off => intrinsic dimension should be "
764 "between 1 and 3.\n";
765 os << "Error: export_to_off => intrinsic dimension should be "
766 "between 1 and 3.\n";
767 return os;
768 }
769
770 std::stringstream output;
771 std::size_t num_simplices, num_vertices;
772 export_vertices_to_off(output, num_vertices, false, point_projection);
773 if (p_complex) {
774 export_simplices_to_off(*p_complex, output, num_simplices, p_simpl_to_color_in_red, p_simpl_to_color_in_green,
775 p_simpl_to_color_in_blue);
776 } else {
777 export_simplices_to_off(output, num_simplices, color_inconsistencies, p_simpl_to_color_in_red,
778 p_simpl_to_color_in_green, p_simpl_to_color_in_blue);
779 }
780
781#ifdef GUDHI_TC_EXPORT_NORMALS
782 os << "N";
783#endif
784
785 os << "OFF \n"
786 << num_vertices << " " << num_simplices << " "
787 << "0 \n"
788 << output.str();
789
790 return os;
791 }
792
793 private:
794 void refresh_tangential_complex() {
795#if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING)
796 std::cerr << yellow << "\nRefreshing TC... " << white;
797#endif
798
799#ifdef GUDHI_TC_PROFILING
800 Gudhi::Clock t;
801#endif
802#ifdef GUDHI_USE_TBB
803 // Parallel
804 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
805 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*this));
806 } else {
807#endif // GUDHI_USE_TBB
808 // Sequential
809 for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
810#ifdef GUDHI_USE_TBB
811 }
812#endif // GUDHI_USE_TBB
813
814#ifdef GUDHI_TC_PROFILING
815 t.end();
816 std::cerr << yellow << "done in " << t.num_seconds() << " seconds.\n" << white;
817#elif defined(DEBUG_TRACES)
818 std::cerr << yellow << "done.\n" << white;
819#endif
820 }
821
822 // If the list of perturbed points is provided, it is much faster
823 template <typename Point_indices_range>
824 void refresh_tangential_complex(Point_indices_range const &perturbed_points_indices) {
825#if defined(DEBUG_TRACES) || defined(GUDHI_TC_PROFILING)
826 std::cerr << yellow << "\nRefreshing TC... " << white;
827#endif
828
829#ifdef GUDHI_TC_PROFILING
830 Gudhi::Clock t;
831#endif
832
833 // ANN tree containing only the perturbed points
834 Points_ds updated_pts_ds(m_points, perturbed_points_indices);
835
836#ifdef GUDHI_USE_TBB
837 // Parallel
838 if (boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value) {
839 tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()),
840 Refresh_tangent_triangulation(*this, updated_pts_ds));
841 } else {
842#endif // GUDHI_USE_TBB
843 // Sequential
844 for (std::size_t i = 0; i < m_points.size(); ++i) refresh_tangent_triangulation(i, updated_pts_ds);
845#ifdef GUDHI_USE_TBB
846 }
847#endif // GUDHI_USE_TBB
848
849#ifdef GUDHI_TC_PROFILING
850 t.end();
851 std::cerr << yellow << "done in " << t.num_seconds() << " seconds.\n" << white;
852#elif defined(DEBUG_TRACES)
853 std::cerr << yellow << "done.\n" << white;
854#endif
855 }
856
857 void export_inconsistent_stars_to_OFF_files(std::string const &filename_base) const {
858 // For each triangulation
859 for (std::size_t idx = 0; idx < m_points.size(); ++idx) {
860 // We build a SC along the way in case it's inconsistent
861 Simplicial_complex sc;
862 // For each cell
863 bool is_inconsistent = false;
864 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
865 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
866 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
867 // Skip infinite cells
868 if (is_infinite(*it_inc_simplex)) continue;
869
870 Simplex c = *it_inc_simplex;
871 c.insert(idx); // Add the missing index
872
873 sc.add_simplex(c);
874
875 // If we do not already know this star is inconsistent, test it
876 if (!is_inconsistent && !is_simplex_consistent(c)) is_inconsistent = true;
877 }
878
879 if (is_inconsistent) {
880 // Export star to OFF file
881 std::stringstream output_filename;
882 output_filename << filename_base << "_" << idx << ".off";
883 std::ofstream off_stream(output_filename.str().c_str());
884 export_to_off(sc, off_stream);
885 }
886 }
887 }
888
889 class Compare_distance_to_ref_point {
890 public:
891 Compare_distance_to_ref_point(Point const &ref, K const &k) : m_ref(ref), m_k(k) {}
892
893 bool operator()(Point const &p1, Point const &p2) {
894 typename K::Squared_distance_d sqdist = m_k.squared_distance_d_object();
895 return sqdist(p1, m_ref) < sqdist(p2, m_ref);
896 }
897
898 private:
899 Point const &m_ref;
900 K const &m_k;
901 };
902
903#ifdef GUDHI_USE_TBB
904 // Functor for compute_tangential_complex function
905 class Compute_tangent_triangulation {
906 Tangential_complex &m_tc;
907
908 public:
909 // Constructor
910 Compute_tangent_triangulation(Tangential_complex &tc) : m_tc(tc) {}
911
912 // Constructor
913 Compute_tangent_triangulation(const Compute_tangent_triangulation &ctt) : m_tc(ctt.m_tc) {}
914
915 // operator()
916 void operator()(const tbb::blocked_range<size_t> &r) const {
917 for (size_t i = r.begin(); i != r.end(); ++i) m_tc.compute_tangent_triangulation(i);
918 }
919 };
920
921 // Functor for refresh_tangential_complex function
922 class Refresh_tangent_triangulation {
923 Tangential_complex &m_tc;
924 Points_ds const &m_updated_pts_ds;
925
926 public:
927 // Constructor
928 Refresh_tangent_triangulation(Tangential_complex &tc, Points_ds const &updated_pts_ds)
929 : m_tc(tc), m_updated_pts_ds(updated_pts_ds) {}
930
931 // Constructor
932 Refresh_tangent_triangulation(const Refresh_tangent_triangulation &ctt)
933 : m_tc(ctt.m_tc), m_updated_pts_ds(ctt.m_updated_pts_ds) {}
934
935 // operator()
936 void operator()(const tbb::blocked_range<size_t> &r) const {
937 for (size_t i = r.begin(); i != r.end(); ++i) m_tc.refresh_tangent_triangulation(i, m_updated_pts_ds);
938 }
939 };
940#endif // GUDHI_USE_TBB
941
942 bool is_infinite(Simplex const &s) const { return *s.rbegin() == (std::numeric_limits<std::size_t>::max)(); }
943
944 // Output: "triangulation" is a Regular Triangulation containing at least the
945 // star of "center_pt"
946 // Returns the handle of the center vertex
947 Tr_vertex_handle compute_star(std::size_t i, const Point &center_pt, const Tangent_space_basis &tsb,
948 Triangulation &triangulation, bool verbose = false) {
949 int tangent_space_dim = tsb.dimension();
950 const Tr_traits &local_tr_traits = triangulation.geom_traits();
951
952 // Kernel functor & objects
953 typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
954
955 // Triangulation's traits functor & objects
956 typename Tr_traits::Compute_weight_d point_weight = local_tr_traits.compute_weight_d_object();
957#if CGAL_VERSION_NR < 1050200000
958 typename Tr_traits::Power_center_d power_center = local_tr_traits.power_center_d_object();
959#else
960 typename Tr_traits::Construct_power_sphere_d power_center = local_tr_traits.construct_power_sphere_d_object();
961#endif
962
963 //***************************************************
964 // Build a minimal triangulation in the tangent space
965 // (we only need the star of p)
966 //***************************************************
967
968 // Insert p
969 Tr_point proj_wp;
970 if (i == tsb.origin()) {
971 // Insert {(0, 0, 0...), m_weights[i]}
972 proj_wp = local_tr_traits.construct_weighted_point_d_object()(
973 local_tr_traits.construct_point_d_object()(tangent_space_dim, CGAL::ORIGIN), m_weights[i]);
974 } else {
975 const Weighted_point &wp = compute_perturbed_weighted_point(i);
976 proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
977 }
978
979 Tr_vertex_handle center_vertex = triangulation.insert(proj_wp);
980 center_vertex->data() = i;
981 if (verbose) std::cerr << "* Inserted point #" << i << "\n";
982
983#ifdef GUDHI_TC_VERY_VERBOSE
984 std::size_t num_attempts_to_insert_points = 1;
985 std::size_t num_inserted_points = 1;
986#endif
987 // const int NUM_NEIGHBORS = 150;
988 // KNS_range ins_range = m_points_ds.k_nearest_neighbors(center_pt, NUM_NEIGHBORS);
989 INS_range ins_range = m_points_ds.incremental_nearest_neighbors(center_pt);
990
991 // While building the local triangulation, we keep the radius
992 // of the sphere "star sphere" centered at "center_vertex"
993 // and which contains all the
994 // circumspheres of the star of "center_vertex"
995 // If th the m_max_squared_edge_length is set the maximal radius of the "star sphere"
996 // is at most square root of m_max_squared_edge_length
997 boost::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length;
998
999 // Insert points until we find a point which is outside "star sphere"
1000 for (auto nn_it = ins_range.begin(); nn_it != ins_range.end(); ++nn_it) {
1001 std::size_t neighbor_point_idx = nn_it->first;
1002
1003 // ith point = p, which is already inserted
1004 if (neighbor_point_idx != i) {
1005 // No need to lock the Mutex_for_perturb here since this will not be
1006 // called while other threads are perturbing the positions
1007 Point neighbor_pt;
1008 FT neighbor_weight;
1009 compute_perturbed_weighted_point(neighbor_point_idx, neighbor_pt, neighbor_weight);
1010 GUDHI_CHECK(!m_max_squared_edge_length ||
1011 squared_star_sphere_radius_plus_margin.value() <= m_max_squared_edge_length.value(),
1012 std::invalid_argument("Tangential_complex::compute_star - set a bigger value with set_max_squared_edge_length."));
1013 if (squared_star_sphere_radius_plus_margin &&
1014 k_sqdist(center_pt, neighbor_pt) > squared_star_sphere_radius_plus_margin.value()) {
1015 GUDHI_CHECK(triangulation.current_dimension() >= tangent_space_dim,
1016 std::invalid_argument("Tangential_complex::compute_star - Dimension of the star is only " + \
1017 std::to_string(triangulation.current_dimension())));
1018 break;
1019 }
1020
1021 Tr_point proj_pt = project_point_and_compute_weight(neighbor_pt, neighbor_weight, tsb, local_tr_traits);
1022
1023#ifdef GUDHI_TC_VERY_VERBOSE
1024 ++num_attempts_to_insert_points;
1025#endif
1026
1027 Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex);
1028 // Tr_vertex_handle vh = triangulation.insert(proj_pt);
1029 if (vh != Tr_vertex_handle() && vh->data() == (std::numeric_limits<std::size_t>::max)()) {
1030#ifdef GUDHI_TC_VERY_VERBOSE
1031 ++num_inserted_points;
1032#endif
1033 if (verbose) std::cerr << "* Inserted point #" << neighbor_point_idx << "\n";
1034
1035 vh->data() = neighbor_point_idx;
1036
1037 // Let's recompute squared_star_sphere_radius_plus_margin
1038 if (triangulation.current_dimension() >= tangent_space_dim) {
1039 squared_star_sphere_radius_plus_margin = boost::none;
1040 // Get the incident cells and look for the biggest circumsphere
1041 std::vector<Tr_full_cell_handle> incident_cells;
1042 triangulation.incident_full_cells(center_vertex, std::back_inserter(incident_cells));
1043 for (typename std::vector<Tr_full_cell_handle>::iterator cit = incident_cells.begin();
1044 cit != incident_cells.end(); ++cit) {
1045 Tr_full_cell_handle cell = *cit;
1046 if (triangulation.is_infinite(cell)) {
1047 squared_star_sphere_radius_plus_margin = boost::none;
1048 break;
1049 } else {
1050 // Note that this uses the perturbed point since it uses
1051 // the points of the local triangulation
1052 Tr_point c =
1053 power_center(boost::make_transform_iterator(cell->vertices_begin(),
1054 vertex_handle_to_point<Tr_point, Tr_vertex_handle>),
1055 boost::make_transform_iterator(cell->vertices_end(),
1056 vertex_handle_to_point<Tr_point, Tr_vertex_handle>));
1057
1058 FT sq_power_sphere_diam = 4 * point_weight(c);
1059
1060 if (!squared_star_sphere_radius_plus_margin ||
1061 sq_power_sphere_diam > squared_star_sphere_radius_plus_margin.value()) {
1062 squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
1063 }
1064 }
1065 }
1066
1067 // Let's add the margin, now
1068 // The value depends on whether we perturb weight or position
1069 if (squared_star_sphere_radius_plus_margin) {
1070 // "2*m_last_max_perturb" because both points can be perturbed
1071 squared_star_sphere_radius_plus_margin =
1072 CGAL::square(std::sqrt(squared_star_sphere_radius_plus_margin.value()) + 2 * m_last_max_perturb);
1073
1074 // Reduce the square radius to m_max_squared_edge_length if necessary
1075 if (m_max_squared_edge_length && squared_star_sphere_radius_plus_margin.value() > m_max_squared_edge_length.value()) {
1076 squared_star_sphere_radius_plus_margin = m_max_squared_edge_length.value();
1077 }
1078
1079 // Save it in `m_squared_star_spheres_radii_incl_margin`
1080 m_squared_star_spheres_radii_incl_margin[i] = squared_star_sphere_radius_plus_margin.value();
1081 } else {
1082 if (m_max_squared_edge_length) {
1083 squared_star_sphere_radius_plus_margin = m_max_squared_edge_length.value();
1084 m_squared_star_spheres_radii_incl_margin[i] = m_max_squared_edge_length.value();
1085 } else {
1086 m_squared_star_spheres_radii_incl_margin[i] = FT(-1);
1087 }
1088 }
1089 }
1090 }
1091 }
1092 }
1093
1094 return center_vertex;
1095 }
1096
1097 void refresh_tangent_triangulation(std::size_t i, Points_ds const &updated_pts_ds, bool verbose = false) {
1098 if (verbose) std::cerr << "** Refreshing tangent tri #" << i << " **\n";
1099
1100 if (m_squared_star_spheres_radii_incl_margin[i] == FT(-1)) return compute_tangent_triangulation(i, verbose);
1101
1102 Point center_point = compute_perturbed_point(i);
1103 // Among updated point, what is the closer from our center point?
1104 std::size_t closest_pt_index = updated_pts_ds.k_nearest_neighbors(center_point, 1, false).begin()->first;
1105
1106 typename K::Construct_weighted_point_d k_constr_wp = m_k.construct_weighted_point_d_object();
1107#if CGAL_VERSION_NR < 1050200000
1108 typename K::Power_distance_d k_power_dist = m_k.power_distance_d_object();
1109#else
1110 typename K::Compute_power_product_d k_power_dist = m_k.compute_power_product_d_object();
1111#endif
1112
1113 // Construct a weighted point equivalent to the star sphere
1114 Weighted_point star_sphere = k_constr_wp(compute_perturbed_point(i), m_squared_star_spheres_radii_incl_margin[i]);
1115 Weighted_point closest_updated_point = compute_perturbed_weighted_point(closest_pt_index);
1116
1117 // Is the "closest point" inside our star sphere?
1118 if (k_power_dist(star_sphere, closest_updated_point) <= FT(0)) compute_tangent_triangulation(i, verbose);
1119 }
1120
1121 void compute_tangent_triangulation(std::size_t i, bool verbose = false) {
1122 if (verbose) std::cerr << "** Computing tangent tri #" << i << " **\n";
1123 // std::cerr << "***********************************************\n";
1124
1125 // No need to lock the mutex here since this will not be called while
1126 // other threads are perturbing the positions
1127 const Point center_pt = compute_perturbed_point(i);
1128 Tangent_space_basis &tsb = m_tangent_spaces[i];
1129
1130 // Estimate the tangent space
1131 if (!m_are_tangent_spaces_computed[i]) {
1132#ifdef GUDHI_TC_EXPORT_NORMALS
1133 tsb = compute_tangent_space(center_pt, i, true /*normalize*/, &m_orth_spaces[i]);
1134#else
1135 tsb = compute_tangent_space(center_pt, i);
1136#endif
1137 }
1138
1139#if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE)
1140 Gudhi::Clock t;
1141#endif
1142 int tangent_space_dim = tangent_basis_dim(i);
1143 Triangulation &local_tr = m_triangulations[i].construct_triangulation(tangent_space_dim);
1144
1145 m_triangulations[i].center_vertex() = compute_star(i, center_pt, tsb, local_tr, verbose);
1146
1147#if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE)
1148 t.end();
1149 std::cerr << " - triangulation construction: " << t.num_seconds() << " s.\n";
1150 t.reset();
1151#endif
1152
1153#ifdef GUDHI_TC_VERY_VERBOSE
1154 std::cerr << "Inserted " << num_inserted_points << " points / " << num_attempts_to_insert_points
1155 << " attempts to compute the star\n";
1156#endif
1157
1158 update_star(i);
1159
1160#if defined(GUDHI_TC_PROFILING) && defined(GUDHI_TC_VERY_VERBOSE)
1161 t.end();
1162 std::cerr << " - update_star: " << t.num_seconds() << " s.\n";
1163#endif
1164 }
1165
1166 // Updates m_stars[i] directly from m_triangulations[i]
1167
1168 void update_star(std::size_t i) {
1169 Star &star = m_stars[i];
1170 star.clear();
1171 Triangulation &local_tr = m_triangulations[i].tr();
1172 Tr_vertex_handle center_vertex = m_triangulations[i].center_vertex();
1173 int cur_dim_plus_1 = local_tr.current_dimension() + 1;
1174
1175 std::vector<Tr_full_cell_handle> incident_cells;
1176 local_tr.incident_full_cells(center_vertex, std::back_inserter(incident_cells));
1177
1178 typename std::vector<Tr_full_cell_handle>::const_iterator it_c = incident_cells.begin();
1179 typename std::vector<Tr_full_cell_handle>::const_iterator it_c_end = incident_cells.end();
1180 // For each cell
1181 for (; it_c != it_c_end; ++it_c) {
1182 // Will contain all indices except center_vertex
1183 Incident_simplex incident_simplex;
1184 for (int j = 0; j < cur_dim_plus_1; ++j) {
1185 std::size_t index = (*it_c)->vertex(j)->data();
1186 if (index != i) incident_simplex.insert(index);
1187 }
1188 GUDHI_CHECK(incident_simplex.size() == cur_dim_plus_1 - 1,
1189 std::logic_error("update_star: wrong size of incident simplex"));
1190 star.push_back(incident_simplex);
1191 }
1192 }
1193
1194 // Estimates tangent subspaces using PCA
1195
1196 Tangent_space_basis compute_tangent_space(const Point &p, const std::size_t i, bool normalize_basis = true,
1197 Orthogonal_space_basis *p_orth_space_basis = NULL) {
1198 unsigned int num_pts_for_pca =
1199 (std::min)(static_cast<unsigned int>(std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1200 static_cast<unsigned int>(m_points.size()));
1201
1202 // Kernel functors
1203 typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object();
1204 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1205
1206#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
1207 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca, false);
1208 const Points &points_for_pca = m_points_for_tse;
1209#else
1210 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca, false);
1211 const Points &points_for_pca = m_points;
1212#endif
1213
1214 // One row = one point
1215 Eigen::MatrixXd mat_points(num_pts_for_pca, m_ambient_dim);
1216 auto nn_it = kns_range.begin();
1217 for (unsigned int j = 0; j < num_pts_for_pca && nn_it != kns_range.end(); ++j, ++nn_it) {
1218 for (int i = 0; i < m_ambient_dim; ++i) {
1219 mat_points(j, i) = CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1220 }
1221 }
1222 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1223 Eigen::MatrixXd cov = centered.adjoint() * centered;
1224 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1225
1226 Tangent_space_basis tsb(i); // p = compute_perturbed_point(i) here
1227
1228 // The eigenvectors are sorted in increasing order of their corresponding
1229 // eigenvalues
1230 for (int j = m_ambient_dim - 1; j >= m_ambient_dim - m_intrinsic_dim; --j) {
1231 if (normalize_basis) {
1232 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1233 eig.eigenvectors().col(j).data() + m_ambient_dim);
1234 tsb.push_back(normalize_vector(v, m_k));
1235 } else {
1236 tsb.push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1237 eig.eigenvectors().col(j).data() + m_ambient_dim));
1238 }
1239 }
1240
1241 if (p_orth_space_basis) {
1242 p_orth_space_basis->set_origin(i);
1243 for (int j = m_ambient_dim - m_intrinsic_dim - 1; j >= 0; --j) {
1244 if (normalize_basis) {
1245 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1246 eig.eigenvectors().col(j).data() + m_ambient_dim);
1247 p_orth_space_basis->push_back(normalize_vector(v, m_k));
1248 } else {
1249 p_orth_space_basis->push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1250 eig.eigenvectors().col(j).data() + m_ambient_dim));
1251 }
1252 }
1253 }
1254
1255 m_are_tangent_spaces_computed[i] = true;
1256
1257 return tsb;
1258 }
1259
1260 // Compute the space tangent to a simplex (p1, p2, ... pn)
1261 // TODO(CJ): Improve this?
1262 // Basically, it takes all the neighbor points to p1, p2... pn and runs PCA
1263 // on it. Note that most points are duplicated.
1264
1265 Tangent_space_basis compute_tangent_space(const Simplex &s, bool normalize_basis = true) {
1266 unsigned int num_pts_for_pca =
1267 (std::min)(static_cast<unsigned int>(std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)),
1268 static_cast<unsigned int>(m_points.size()));
1269
1270 // Kernel functors
1271 typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object();
1272 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1273 typename K::Squared_length_d sqlen = m_k.squared_length_d_object();
1274 typename K::Scaled_vector_d scaled_vec = m_k.scaled_vector_d_object();
1275 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1276 typename K::Difference_of_vectors_d diff_vec = m_k.difference_of_vectors_d_object();
1277
1278 // One row = one point
1279 Eigen::MatrixXd mat_points(s.size() * num_pts_for_pca, m_ambient_dim);
1280 unsigned int current_row = 0;
1281
1282 for (Simplex::const_iterator it_index = s.begin(); it_index != s.end(); ++it_index) {
1283 const Point &p = m_points[*it_index];
1284
1285#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
1286 KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca, false);
1287 const Points &points_for_pca = m_points_for_tse;
1288#else
1289 KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca, false);
1290 const Points &points_for_pca = m_points;
1291#endif
1292
1293 auto nn_it = kns_range.begin();
1294 for (; current_row < num_pts_for_pca && nn_it != kns_range.end(); ++current_row, ++nn_it) {
1295 for (int i = 0; i < m_ambient_dim; ++i) {
1296 mat_points(current_row, i) = CGAL::to_double(coord(points_for_pca[nn_it->first], i));
1297 }
1298 }
1299 }
1300 Eigen::MatrixXd centered = mat_points.rowwise() - mat_points.colwise().mean();
1301 Eigen::MatrixXd cov = centered.adjoint() * centered;
1302 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
1303
1304 Tangent_space_basis tsb;
1305
1306 // The eigenvectors are sorted in increasing order of their corresponding
1307 // eigenvalues
1308 for (int j = m_ambient_dim - 1; j >= m_ambient_dim - m_intrinsic_dim; --j) {
1309 if (normalize_basis) {
1310 Vector v = constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1311 eig.eigenvectors().col(j).data() + m_ambient_dim);
1312 tsb.push_back(normalize_vector(v, m_k));
1313 } else {
1314 tsb.push_back(constr_vec(m_ambient_dim, eig.eigenvectors().col(j).data(),
1315 eig.eigenvectors().col(j).data() + m_ambient_dim));
1316 }
1317 }
1318
1319 return tsb;
1320 }
1321
1322 // Returns the dimension of the ith local triangulation
1323
1324 int tangent_basis_dim(std::size_t i) const { return m_tangent_spaces[i].dimension(); }
1325
1326 Point compute_perturbed_point(std::size_t pt_idx) const {
1327#ifdef GUDHI_TC_PERTURB_POSITION
1328 return m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]);
1329#else
1330 return m_points[pt_idx];
1331#endif
1332 }
1333
1334 void compute_perturbed_weighted_point(std::size_t pt_idx, Point &p, FT &w) const {
1335#ifdef GUDHI_TC_PERTURB_POSITION
1336 p = m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]);
1337#else
1338 p = m_points[pt_idx];
1339#endif
1340 w = m_weights[pt_idx];
1341 }
1342
1343 Weighted_point compute_perturbed_weighted_point(std::size_t pt_idx) const {
1344 typename K::Construct_weighted_point_d k_constr_wp = m_k.construct_weighted_point_d_object();
1345
1346 Weighted_point wp = k_constr_wp(
1347#ifdef GUDHI_TC_PERTURB_POSITION
1348 m_k.translated_point_d_object()(m_points[pt_idx], m_translations[pt_idx]),
1349#else
1350 m_points[pt_idx],
1351#endif
1352 m_weights[pt_idx]);
1353
1354 return wp;
1355 }
1356
1357 Point unproject_point(const Tr_point &p, const Tangent_space_basis &tsb, const Tr_traits &tr_traits) const {
1358 typename K::Translated_point_d k_transl = m_k.translated_point_d_object();
1359 typename K::Scaled_vector_d k_scaled_vec = m_k.scaled_vector_d_object();
1360 typename Tr_traits::Compute_coordinate_d coord = tr_traits.compute_coordinate_d_object();
1361
1362 Point global_point = compute_perturbed_point(tsb.origin());
1363 for (int i = 0; i < m_intrinsic_dim; ++i) global_point = k_transl(global_point, k_scaled_vec(tsb[i], coord(p, i)));
1364
1365 return global_point;
1366 }
1367
1368 // Project the point in the tangent space
1369 // Resulting point coords are expressed in tsb's space
1370 Tr_bare_point project_point(const Point &p, const Tangent_space_basis &tsb, const Tr_traits &tr_traits) const {
1371 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1372 typename K::Difference_of_points_d diff_points = m_k.difference_of_points_d_object();
1373
1374 Vector v = diff_points(p, compute_perturbed_point(tsb.origin()));
1375
1376 std::vector<FT> coords;
1377 // Ambiant-space coords of the projected point
1378 coords.reserve(tsb.dimension());
1379 for (std::size_t i = 0; i < m_intrinsic_dim; ++i) {
1380 // Local coords are given by the scalar product with the vectors of tsb
1381 FT coord = scalar_pdct(v, tsb[i]);
1382 coords.push_back(coord);
1383 }
1384
1385 return tr_traits.construct_point_d_object()(static_cast<int>(coords.size()), coords.begin(), coords.end());
1386 }
1387
1388 // Project the point in the tangent space
1389 // The weight will be the squared distance between p and the projection of p
1390 // Resulting point coords are expressed in tsb's space
1391
1392 Tr_point project_point_and_compute_weight(const Weighted_point &wp, const Tangent_space_basis &tsb,
1393 const Tr_traits &tr_traits) const {
1394 typename K::Point_drop_weight_d k_drop_w = m_k.point_drop_weight_d_object();
1395 typename K::Compute_weight_d k_point_weight = m_k.compute_weight_d_object();
1396 return project_point_and_compute_weight(k_drop_w(wp), k_point_weight(wp), tsb, tr_traits);
1397 }
1398
1399 // Same as above, with slightly different parameters
1400 Tr_point project_point_and_compute_weight(const Point &p, const FT w, const Tangent_space_basis &tsb,
1401 const Tr_traits &tr_traits) const {
1402 const int point_dim = m_k.point_dimension_d_object()(p);
1403
1404 typename K::Construct_point_d constr_pt = m_k.construct_point_d_object();
1405 typename K::Scalar_product_d scalar_pdct = m_k.scalar_product_d_object();
1406 typename K::Difference_of_points_d diff_points = m_k.difference_of_points_d_object();
1407 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1408 typename K::Construct_cartesian_const_iterator_d ccci = m_k.construct_cartesian_const_iterator_d_object();
1409
1410 Point origin = compute_perturbed_point(tsb.origin());
1411 Vector v = diff_points(p, origin);
1412
1413 // Same dimension? Then the weight is 0
1414 bool same_dim = (point_dim == tsb.dimension());
1415
1416 std::vector<FT> coords;
1417 // Ambiant-space coords of the projected point
1418 std::vector<FT> p_proj(ccci(origin), ccci(origin, 0));
1419 coords.reserve(tsb.dimension());
1420 for (int i = 0; i < tsb.dimension(); ++i) {
1421 // Local coords are given by the scalar product with the vectors of tsb
1422 FT c = scalar_pdct(v, tsb[i]);
1423 coords.push_back(c);
1424
1425 // p_proj += c * tsb[i]
1426 if (!same_dim) {
1427 for (int j = 0; j < point_dim; ++j) p_proj[j] += c * coord(tsb[i], j);
1428 }
1429 }
1430
1431 // Same dimension? Then the weight is 0
1432 FT sq_dist_to_proj_pt = 0;
1433 if (!same_dim) {
1434 Point projected_pt = constr_pt(point_dim, p_proj.begin(), p_proj.end());
1435 sq_dist_to_proj_pt = m_k.squared_distance_d_object()(p, projected_pt);
1436 }
1437
1438 return tr_traits.construct_weighted_point_d_object()(
1439 tr_traits.construct_point_d_object()(static_cast<int>(coords.size()), coords.begin(), coords.end()),
1440 w - sq_dist_to_proj_pt);
1441 }
1442
1443 // Project all the points in the tangent space
1444
1445 template <typename Indexed_point_range>
1446 std::vector<Tr_point> project_points_and_compute_weights(const Indexed_point_range &point_indices,
1447 const Tangent_space_basis &tsb,
1448 const Tr_traits &tr_traits) const {
1449 std::vector<Tr_point> ret;
1450 for (typename Indexed_point_range::const_iterator it = point_indices.begin(), it_end = point_indices.end();
1451 it != it_end; ++it) {
1452 ret.push_back(project_point_and_compute_weight(compute_perturbed_weighted_point(*it), tsb, tr_traits));
1453 }
1454 return ret;
1455 }
1456
1457 // A simplex here is a local tri's full cell handle
1458
1459 bool is_simplex_consistent(Tr_full_cell_handle fch, int cur_dim) const {
1460 Simplex c;
1461 for (int i = 0; i < cur_dim + 1; ++i) {
1462 std::size_t data = fch->vertex(i)->data();
1463 c.insert(data);
1464 }
1465 return is_simplex_consistent(c);
1466 }
1467
1468 // A simplex here is a list of point indices
1469 // TODO(CJ): improve it like the other "is_simplex_consistent" below
1470
1471 bool is_simplex_consistent(Simplex const &simplex) const {
1472 // Check if the simplex is in the stars of all its vertices
1473 Simplex::const_iterator it_point_idx = simplex.begin();
1474 // For each point p of the simplex, we parse the incidents cells of p
1475 // and we check if "simplex" is among them
1476 for (; it_point_idx != simplex.end(); ++it_point_idx) {
1477 std::size_t point_idx = *it_point_idx;
1478 // Don't check infinite simplices
1479 if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue;
1480
1481 Star const &star = m_stars[point_idx];
1482
1483 // What we're looking for is "simplex" \ point_idx
1484 Incident_simplex is_to_find = simplex;
1485 is_to_find.erase(point_idx);
1486
1487 // For each cell
1488 if (std::find(star.begin(), star.end(), is_to_find) == star.end()) return false;
1489 }
1490
1491 return true;
1492 }
1493
1494 // A simplex here is a list of point indices
1495 // "s" contains all the points of the simplex except "center_point"
1496 // This function returns the points whose star doesn't contain the simplex
1497 // N.B.: the function assumes that the simplex is contained in
1498 // star(center_point)
1499
1500 template <typename OutputIterator> // value_type = std::size_t
1501 bool is_simplex_consistent(std::size_t center_point,
1502 Incident_simplex const &s, // without "center_point"
1503 OutputIterator points_whose_star_does_not_contain_s,
1504 bool check_also_in_non_maximal_faces = false) const {
1505 Simplex full_simplex = s;
1506 full_simplex.insert(center_point);
1507
1508 // Check if the simplex is in the stars of all its vertices
1509 Incident_simplex::const_iterator it_point_idx = s.begin();
1510 // For each point p of the simplex, we parse the incidents cells of p
1511 // and we check if "simplex" is among them
1512 for (; it_point_idx != s.end(); ++it_point_idx) {
1513 std::size_t point_idx = *it_point_idx;
1514 // Don't check infinite simplices
1515 if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue;
1516
1517 Star const &star = m_stars[point_idx];
1518
1519 // What we're looking for is full_simplex \ point_idx
1520 Incident_simplex is_to_find = full_simplex;
1521 is_to_find.erase(point_idx);
1522
1523 if (check_also_in_non_maximal_faces) {
1524 // For each simplex "is" of the star, check if ic_to_simplex is
1525 // included in "is"
1526 bool found = false;
1527 for (Star::const_iterator is = star.begin(), is_end = star.end(); !found && is != is_end; ++is) {
1528 if (std::includes(is->begin(), is->end(), is_to_find.begin(), is_to_find.end())) found = true;
1529 }
1530
1531 if (!found) *points_whose_star_does_not_contain_s++ = point_idx;
1532 } else {
1533 // Does the star contain is_to_find?
1534 if (std::find(star.begin(), star.end(), is_to_find) == star.end())
1535 *points_whose_star_does_not_contain_s++ = point_idx;
1536 }
1537 }
1538
1539 return true;
1540 }
1541
1542 // A simplex here is a list of point indices
1543 // It looks for s in star(p).
1544 // "s" contains all the points of the simplex except p.
1545 bool is_simplex_in_star(std::size_t p, Incident_simplex const &s, bool check_also_in_non_maximal_faces = true) const {
1546 Star const &star = m_stars[p];
1547
1548 if (check_also_in_non_maximal_faces) {
1549 // For each simplex "is" of the star, check if ic_to_simplex is
1550 // included in "is"
1551 bool found = false;
1552 for (Star::const_iterator is = star.begin(), is_end = star.end(); !found && is != is_end; ++is) {
1553 if (std::includes(is->begin(), is->end(), s.begin(), s.end())) found = true;
1554 }
1555
1556 return found;
1557 } else {
1558 return !(std::find(star.begin(), star.end(), s) == star.end());
1559 }
1560 }
1561
1562#ifdef GUDHI_USE_TBB
1563 // Functor for try_to_solve_inconsistencies_in_a_local_triangulation function
1564 class Try_to_solve_inconsistencies_in_a_local_triangulation {
1565 Tangential_complex &m_tc;
1566 double m_max_perturb;
1567 tbb::combinable<std::size_t> &m_num_inconsistencies;
1568 tbb::combinable<std::vector<std::size_t> > &m_updated_points;
1569
1570 public:
1571 // Constructor
1572 Try_to_solve_inconsistencies_in_a_local_triangulation(Tangential_complex &tc, double max_perturb,
1573 tbb::combinable<std::size_t> &num_inconsistencies,
1574 tbb::combinable<std::vector<std::size_t> > &updated_points)
1575 : m_tc(tc),
1576 m_max_perturb(max_perturb),
1577 m_num_inconsistencies(num_inconsistencies),
1578 m_updated_points(updated_points) {}
1579
1580 // Constructor
1581 Try_to_solve_inconsistencies_in_a_local_triangulation(
1582 const Try_to_solve_inconsistencies_in_a_local_triangulation &tsilt)
1583 : m_tc(tsilt.m_tc),
1584 m_max_perturb(tsilt.m_max_perturb),
1585 m_num_inconsistencies(tsilt.m_num_inconsistencies),
1586 m_updated_points(tsilt.m_updated_points) {}
1587
1588 // operator()
1589 void operator()(const tbb::blocked_range<size_t> &r) const {
1590 for (size_t i = r.begin(); i != r.end(); ++i) {
1591 m_num_inconsistencies.local() += m_tc.try_to_solve_inconsistencies_in_a_local_triangulation(
1592 i, m_max_perturb, std::back_inserter(m_updated_points.local()));
1593 }
1594 }
1595 };
1596#endif // GUDHI_USE_TBB
1597
1598 void perturb(std::size_t point_idx, double max_perturb) {
1599 const Tr_traits &local_tr_traits = m_triangulations[point_idx].tr().geom_traits();
1600 typename Tr_traits::Compute_coordinate_d coord = local_tr_traits.compute_coordinate_d_object();
1601 typename K::Translated_point_d k_transl = m_k.translated_point_d_object();
1602 typename K::Construct_vector_d k_constr_vec = m_k.construct_vector_d_object();
1603 typename K::Scaled_vector_d k_scaled_vec = m_k.scaled_vector_d_object();
1604
1605 CGAL::Random_points_in_ball_d<Tr_bare_point> tr_point_in_ball_generator(
1606 m_intrinsic_dim, m_random_generator.get_double(0., max_perturb));
1607
1608 Tr_point local_random_transl =
1609 local_tr_traits.construct_weighted_point_d_object()(*tr_point_in_ball_generator++, 0);
1610 Translation_for_perturb global_transl = k_constr_vec(m_ambient_dim);
1611 const Tangent_space_basis &tsb = m_tangent_spaces[point_idx];
1612 for (int i = 0; i < m_intrinsic_dim; ++i) {
1613 global_transl = k_transl(global_transl, k_scaled_vec(tsb[i], coord(local_random_transl, i)));
1614 }
1615 // Parallel
1616#if defined(GUDHI_USE_TBB)
1617 m_p_perturb_mutexes[point_idx].lock();
1618 m_translations[point_idx] = global_transl;
1619 m_p_perturb_mutexes[point_idx].unlock();
1620 // Sequential
1621#else
1622 m_translations[point_idx] = global_transl;
1623#endif
1624 }
1625
1626 // Return true if inconsistencies were found
1627 template <typename OutputIt>
1628 bool try_to_solve_inconsistencies_in_a_local_triangulation(
1629 std::size_t tr_index, double max_perturb, OutputIt perturbed_pts_indices = CGAL::Emptyset_iterator()) {
1630 bool is_inconsistent = false;
1631
1632 Star const &star = m_stars[tr_index];
1633
1634 // For each incident simplex
1635 Star::const_iterator it_inc_simplex = star.begin();
1636 Star::const_iterator it_inc_simplex_end = star.end();
1637 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1638 const Incident_simplex &incident_simplex = *it_inc_simplex;
1639
1640 // Don't check infinite cells
1641 if (is_infinite(incident_simplex)) continue;
1642
1643 Simplex c = incident_simplex;
1644 c.insert(tr_index); // Add the missing index
1645
1646 // Perturb the center point
1647 if (!is_simplex_consistent(c)) {
1648 is_inconsistent = true;
1649
1650 std::size_t idx = tr_index;
1651
1652 perturb(tr_index, max_perturb);
1653 *perturbed_pts_indices++ = idx;
1654
1655 // We will try the other cells next time
1656 break;
1657 }
1658 }
1659
1660 return is_inconsistent;
1661 }
1662
1663 // 1st line: number of points
1664 // Then one point per line
1665 std::ostream &export_point_set(std::ostream &os, bool use_perturbed_points = false,
1666 const char *coord_separator = " ") const {
1667 if (use_perturbed_points) {
1668 std::vector<Point> perturbed_points;
1669 perturbed_points.reserve(m_points.size());
1670 for (std::size_t i = 0; i < m_points.size(); ++i) perturbed_points.push_back(compute_perturbed_point(i));
1671
1672 return export_point_set(m_k, perturbed_points, os, coord_separator);
1673 } else {
1674 return export_point_set(m_k, m_points, os, coord_separator);
1675 }
1676 }
1677
1678 template <typename ProjectionFunctor = CGAL::Identity<Point> >
1679 std::ostream &export_vertices_to_off(std::ostream &os, std::size_t &num_vertices, bool use_perturbed_points = false,
1680 ProjectionFunctor const &point_projection = ProjectionFunctor()) const {
1681 if (m_points.empty()) {
1682 num_vertices = 0;
1683 return os;
1684 }
1685
1686 // If m_intrinsic_dim = 1, we output each point two times
1687 // to be able to export each segment as a flat triangle with 3 different
1688 // indices (otherwise, Meshlab detects degenerated simplices)
1689 const int N = (m_intrinsic_dim == 1 ? 2 : 1);
1690
1691 // Kernel functors
1692 typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object();
1693
1694#ifdef GUDHI_TC_EXPORT_ALL_COORDS_IN_OFF
1695 int num_coords = m_ambient_dim;
1696#else
1697 int num_coords = (std::min)(m_ambient_dim, 3);
1698#endif
1699
1700#ifdef GUDHI_TC_EXPORT_NORMALS
1701 OS_container::const_iterator it_os = m_orth_spaces.begin();
1702#endif
1703 typename Points::const_iterator it_p = m_points.begin();
1704 typename Points::const_iterator it_p_end = m_points.end();
1705 // For each point p
1706 for (std::size_t i = 0; it_p != it_p_end; ++it_p, ++i) {
1707 Point p = point_projection(use_perturbed_points ? compute_perturbed_point(i) : *it_p);
1708 for (int ii = 0; ii < N; ++ii) {
1709 int j = 0;
1710 for (; j < num_coords; ++j) os << CGAL::to_double(coord(p, j)) << " ";
1711 if (j == 2) os << "0";
1712
1713#ifdef GUDHI_TC_EXPORT_NORMALS
1714 for (j = 0; j < num_coords; ++j) os << " " << CGAL::to_double(coord(*it_os->begin(), j));
1715#endif
1716 os << "\n";
1717 }
1718#ifdef GUDHI_TC_EXPORT_NORMALS
1719 ++it_os;
1720#endif
1721 }
1722
1723 num_vertices = N * m_points.size();
1724 return os;
1725 }
1726
1727 std::ostream &export_simplices_to_off(std::ostream &os, std::size_t &num_OFF_simplices,
1728 bool color_inconsistencies = false,
1729 Simplex_set const *p_simpl_to_color_in_red = NULL,
1730 Simplex_set const *p_simpl_to_color_in_green = NULL,
1731 Simplex_set const *p_simpl_to_color_in_blue = NULL) const {
1732 // If m_intrinsic_dim = 1, each point is output two times
1733 // (see export_vertices_to_off)
1734 num_OFF_simplices = 0;
1735 std::size_t num_maximal_simplices = 0;
1736 std::size_t num_inconsistent_maximal_simplices = 0;
1737 std::size_t num_inconsistent_stars = 0;
1738 typename Tr_container::const_iterator it_tr = m_triangulations.begin();
1739 typename Tr_container::const_iterator it_tr_end = m_triangulations.end();
1740 // For each triangulation
1741 for (std::size_t idx = 0; it_tr != it_tr_end; ++it_tr, ++idx) {
1742 bool is_star_inconsistent = false;
1743
1744 Triangulation const &tr = it_tr->tr();
1745
1746 if (tr.current_dimension() < m_intrinsic_dim) continue;
1747
1748 // Color for this star
1749 std::stringstream color;
1750 // color << rand()%256 << " " << 100+rand()%156 << " " << 100+rand()%156;
1751 color << 128 << " " << 128 << " " << 128;
1752
1753 // Gather the triangles here, with an int telling its color
1754 typedef std::vector<std::pair<Simplex, int> > Star_using_triangles;
1755 Star_using_triangles star_using_triangles;
1756
1757 // For each cell of the star
1758 Star::const_iterator it_inc_simplex = m_stars[idx].begin();
1759 Star::const_iterator it_inc_simplex_end = m_stars[idx].end();
1760 for (; it_inc_simplex != it_inc_simplex_end; ++it_inc_simplex) {
1761 Simplex c = *it_inc_simplex;
1762 c.insert(idx);
1763 std::size_t num_vertices = c.size();
1764 ++num_maximal_simplices;
1765
1766 int color_simplex = -1; // -1=no color, 0=yellow, 1=red, 2=green, 3=blue
1767 if (color_inconsistencies && !is_simplex_consistent(c)) {
1768 ++num_inconsistent_maximal_simplices;
1769 color_simplex = 0;
1770 is_star_inconsistent = true;
1771 } else {
1772 if (p_simpl_to_color_in_red && std::find(p_simpl_to_color_in_red->begin(), p_simpl_to_color_in_red->end(),
1773 c) != p_simpl_to_color_in_red->end()) {
1774 color_simplex = 1;
1775 } else if (p_simpl_to_color_in_green &&
1776 std::find(p_simpl_to_color_in_green->begin(), p_simpl_to_color_in_green->end(), c) !=
1777 p_simpl_to_color_in_green->end()) {
1778 color_simplex = 2;
1779 } else if (p_simpl_to_color_in_blue &&
1780 std::find(p_simpl_to_color_in_blue->begin(), p_simpl_to_color_in_blue->end(), c) !=
1781 p_simpl_to_color_in_blue->end()) {
1782 color_simplex = 3;
1783 }
1784 }
1785
1786 // If m_intrinsic_dim = 1, each point is output two times,
1787 // so we need to multiply each index by 2
1788 // And if only 2 vertices, add a third one (each vertex is duplicated in
1789 // the file when m_intrinsic dim = 2)
1790 if (m_intrinsic_dim == 1) {
1791 Simplex tmp_c;
1792 Simplex::iterator it = c.begin();
1793 for (; it != c.end(); ++it) tmp_c.insert(*it * 2);
1794 if (num_vertices == 2) tmp_c.insert(*tmp_c.rbegin() + 1);
1795
1796 c = tmp_c;
1797 }
1798
1799 if (num_vertices <= 3) {
1800 star_using_triangles.push_back(std::make_pair(c, color_simplex));
1801 } else {
1802 // num_vertices >= 4: decompose the simplex in triangles
1803 std::vector<bool> booleans(num_vertices, false);
1804 std::fill(booleans.begin() + num_vertices - 3, booleans.end(), true);
1805 do {
1806 Simplex triangle;
1807 Simplex::iterator it = c.begin();
1808 for (int i = 0; it != c.end(); ++i, ++it) {
1809 if (booleans[i]) triangle.insert(*it);
1810 }
1811 star_using_triangles.push_back(std::make_pair(triangle, color_simplex));
1812 } while (std::next_permutation(booleans.begin(), booleans.end()));
1813 }
1814 }
1815
1816 // For each cell
1817 Star_using_triangles::const_iterator it_simplex = star_using_triangles.begin();
1818 Star_using_triangles::const_iterator it_simplex_end = star_using_triangles.end();
1819 for (; it_simplex != it_simplex_end; ++it_simplex) {
1820 const Simplex &c = it_simplex->first;
1821
1822 // Don't export infinite cells
1823 if (is_infinite(c)) continue;
1824
1825 int color_simplex = it_simplex->second;
1826
1827 std::stringstream sstr_c;
1828
1829 Simplex::const_iterator it_point_idx = c.begin();
1830 for (; it_point_idx != c.end(); ++it_point_idx) {
1831 sstr_c << *it_point_idx << " ";
1832 }
1833
1834 os << 3 << " " << sstr_c.str();
1835 if (color_inconsistencies || p_simpl_to_color_in_red || p_simpl_to_color_in_green || p_simpl_to_color_in_blue) {
1836 switch (color_simplex) {
1837 case 0:
1838 os << " 255 255 0";
1839 break;
1840 case 1:
1841 os << " 255 0 0";
1842 break;
1843 case 2:
1844 os << " 0 255 0";
1845 break;
1846 case 3:
1847 os << " 0 0 255";
1848 break;
1849 default:
1850 os << " " << color.str();
1851 break;
1852 }
1853 }
1854 ++num_OFF_simplices;
1855 os << "\n";
1856 }
1857 if (is_star_inconsistent) ++num_inconsistent_stars;
1858 }
1859
1860#ifdef DEBUG_TRACES
1861 std::cerr << "\n==========================================================\n"
1862 << "Export from list of stars to OFF:\n"
1863 << " * Number of vertices: " << m_points.size() << "\n"
1864 << " * Total number of maximal simplices: " << num_maximal_simplices << "\n";
1865 if (color_inconsistencies) {
1866 std::cerr << " * Number of inconsistent stars: " << num_inconsistent_stars << " ("
1867 << (m_points.size() > 0 ? 100. * num_inconsistent_stars / m_points.size() : 0.) << "%)\n"
1868 << " * Number of inconsistent maximal simplices: " << num_inconsistent_maximal_simplices << " ("
1869 << (num_maximal_simplices > 0 ? 100. * num_inconsistent_maximal_simplices / num_maximal_simplices : 0.)
1870 << "%)\n";
1871 }
1872 std::cerr << "==========================================================\n";
1873#endif
1874
1875 return os;
1876 }
1877
1878 public:
1879 std::ostream &export_simplices_to_off(const Simplicial_complex &complex, std::ostream &os,
1880 std::size_t &num_OFF_simplices,
1881 Simplex_set const *p_simpl_to_color_in_red = NULL,
1882 Simplex_set const *p_simpl_to_color_in_green = NULL,
1883 Simplex_set const *p_simpl_to_color_in_blue = NULL) const {
1884 typedef Simplicial_complex::Simplex Simplex;
1885 typedef Simplicial_complex::Simplex_set Simplex_set;
1886
1887 // If m_intrinsic_dim = 1, each point is output two times
1888 // (see export_vertices_to_off)
1889 num_OFF_simplices = 0;
1890 std::size_t num_maximal_simplices = 0;
1891
1892 typename Simplex_set::const_iterator it_s = complex.simplex_range().begin();
1893 typename Simplex_set::const_iterator it_s_end = complex.simplex_range().end();
1894 // For each simplex
1895 for (; it_s != it_s_end; ++it_s) {
1896 Simplex c = *it_s;
1897 ++num_maximal_simplices;
1898
1899 int color_simplex = -1; // -1=no color, 0=yellow, 1=red, 2=green, 3=blue
1900 if (p_simpl_to_color_in_red && std::find(p_simpl_to_color_in_red->begin(), p_simpl_to_color_in_red->end(), c) !=
1901 p_simpl_to_color_in_red->end()) {
1902 color_simplex = 1;
1903 } else if (p_simpl_to_color_in_green &&
1904 std::find(p_simpl_to_color_in_green->begin(), p_simpl_to_color_in_green->end(), c) !=
1905 p_simpl_to_color_in_green->end()) {
1906 color_simplex = 2;
1907 } else if (p_simpl_to_color_in_blue &&
1908 std::find(p_simpl_to_color_in_blue->begin(), p_simpl_to_color_in_blue->end(), c) !=
1909 p_simpl_to_color_in_blue->end()) {
1910 color_simplex = 3;
1911 }
1912
1913 // Gather the triangles here
1914 typedef std::vector<Simplex> Triangles;
1915 Triangles triangles;
1916
1917 int num_vertices = static_cast<int>(c.size());
1918 // Do not export smaller dimension simplices
1919 if (num_vertices < m_intrinsic_dim + 1) continue;
1920
1921 // If m_intrinsic_dim = 1, each point is output two times,
1922 // so we need to multiply each index by 2
1923 // And if only 2 vertices, add a third one (each vertex is duplicated in
1924 // the file when m_intrinsic dim = 2)
1925 if (m_intrinsic_dim == 1) {
1926 Simplex tmp_c;
1927 Simplex::iterator it = c.begin();
1928 for (; it != c.end(); ++it) tmp_c.insert(*it * 2);
1929 if (num_vertices == 2) tmp_c.insert(*tmp_c.rbegin() + 1);
1930
1931 c = tmp_c;
1932 }
1933
1934 if (num_vertices <= 3) {
1935 triangles.push_back(c);
1936 } else {
1937 // num_vertices >= 4: decompose the simplex in triangles
1938 std::vector<bool> booleans(num_vertices, false);
1939 std::fill(booleans.begin() + num_vertices - 3, booleans.end(), true);
1940 do {
1941 Simplex triangle;
1942 Simplex::iterator it = c.begin();
1943 for (int i = 0; it != c.end(); ++i, ++it) {
1944 if (booleans[i]) triangle.insert(*it);
1945 }
1946 triangles.push_back(triangle);
1947 } while (std::next_permutation(booleans.begin(), booleans.end()));
1948 }
1949
1950 // For each cell
1951 Triangles::const_iterator it_tri = triangles.begin();
1952 Triangles::const_iterator it_tri_end = triangles.end();
1953 for (; it_tri != it_tri_end; ++it_tri) {
1954 // Don't export infinite cells
1955 if (is_infinite(*it_tri)) continue;
1956
1957 os << 3 << " ";
1958 Simplex::const_iterator it_point_idx = it_tri->begin();
1959 for (; it_point_idx != it_tri->end(); ++it_point_idx) {
1960 os << *it_point_idx << " ";
1961 }
1962
1963 if (p_simpl_to_color_in_red || p_simpl_to_color_in_green || p_simpl_to_color_in_blue) {
1964 switch (color_simplex) {
1965 case 0:
1966 os << " 255 255 0";
1967 break;
1968 case 1:
1969 os << " 255 0 0";
1970 break;
1971 case 2:
1972 os << " 0 255 0";
1973 break;
1974 case 3:
1975 os << " 0 0 255";
1976 break;
1977 default:
1978 os << " 128 128 128";
1979 break;
1980 }
1981 }
1982
1983 ++num_OFF_simplices;
1984 os << "\n";
1985 }
1986 }
1987
1988#ifdef DEBUG_TRACES
1989 std::cerr << "\n==========================================================\n"
1990 << "Export from complex to OFF:\n"
1991 << " * Number of vertices: " << m_points.size() << "\n"
1992 << " * Total number of maximal simplices: " << num_maximal_simplices << "\n"
1993 << "==========================================================\n";
1994#endif
1995
1996 return os;
1997 }
1998
2006 void set_max_squared_edge_length(FT max_squared_edge_length) { m_max_squared_edge_length = max_squared_edge_length; }
2007
2008 private:
2009 const K m_k;
2010 const int m_intrinsic_dim;
2011 const int m_ambient_dim;
2012
2013 Points m_points;
2014 Weights m_weights;
2015#ifdef GUDHI_TC_PERTURB_POSITION
2016 Translations_for_perturb m_translations;
2017#if defined(GUDHI_USE_TBB)
2018 Mutex_for_perturb *m_p_perturb_mutexes;
2019#endif
2020#endif
2021
2022 Points_ds m_points_ds;
2023 double m_last_max_perturb;
2024 std::vector<bool> m_are_tangent_spaces_computed;
2025 TS_container m_tangent_spaces;
2026#ifdef GUDHI_TC_EXPORT_NORMALS
2027 OS_container m_orth_spaces;
2028#endif
2029 Tr_container m_triangulations; // Contains the triangulations
2030 // and their center vertex
2031 Stars_container m_stars;
2032 std::vector<FT> m_squared_star_spheres_radii_incl_margin;
2033 boost::optional<FT> m_max_squared_edge_length;
2034
2035#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
2036 Points m_points_for_tse;
2037 Points_ds m_points_ds_for_tse;
2038#endif
2039
2040 mutable CGAL::Random m_random_generator;
2041}; // /class Tangential_complex
2042
2043} // end namespace tangential_complex
2044} // end namespace Gudhi
2045
2046#endif // TANGENTIAL_COMPLEX_H_
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:104
INS_range incremental_nearest_neighbors(Point const &p, FT eps=FT(0)) const
Search incrementally for the nearest neighbors from a query point.
Definition: Kd_tree_search.h:271
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:245
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:112
Tangential complex data structure.
Definition: Tangential_complex.h:125
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:240
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:395
void compute_tangential_complex()
Computes the tangential complex.
Definition: Tangential_complex.h:330
int intrinsic_dimension() const
Returns the intrinsic dimension of the manifold.
Definition: Tangential_complex.h:279
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:298
void set_max_squared_edge_length(FT max_squared_edge_length)
Sets the maximal possible squared edge length for the edges in the triangulations.
Definition: Tangential_complex.h:2006
Point get_point(std::size_t vertex) const
Returns the point corresponding to the vertex given as parameter.
Definition: Tangential_complex.h:291
Num_inconsistencies number_of_inconsistent_simplices(bool verbose=false) const
Definition: Tangential_complex.h:551
int ambient_dimension() const
Returns the ambient dimension.
Definition: Tangential_complex.h:282
std::size_t number_of_vertices() const
Returns the number of vertices.
Definition: Tangential_complex.h:302
int create_complex(Simplex_tree_ &tree, bool export_inconsistent_simplices=true) const
Exports the complex into a Simplex_tree.
Definition: Tangential_complex.h:610
~Tangential_complex()
Destructor.
Definition: Tangential_complex.h:272
Type returned by Tangential_complex::fix_inconsistencies_using_perturbation.
Definition: Tangential_complex.h:377
std::size_t final_num_inconsistent_stars
final number of inconsistent stars
Definition: Tangential_complex.h:387
bool success
true if all inconsistencies could be removed, false if the time limit has been reached before
Definition: Tangential_complex.h:379
unsigned int num_steps
number of steps performed
Definition: Tangential_complex.h:381
std::size_t best_num_inconsistent_stars
best number of inconsistent stars during the process
Definition: Tangential_complex.h:385
std::size_t initial_num_inconsistent_stars
initial number of inconsistent stars
Definition: Tangential_complex.h:383
Type returned by Tangential_complex::number_of_inconsistent_simplices.
Definition: Tangential_complex.h:539
std::size_t num_inconsistent_simplices
Number of inconsistent simplices.
Definition: Tangential_complex.h:543
std::size_t num_inconsistent_stars
Number of stars containing at least one inconsistent simplex.
Definition: Tangential_complex.h:545
std::size_t num_simplices
Total number of simplices in stars (including duplicates that appear in several stars)
Definition: Tangential_complex.h:541