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