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