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