Merge branch 'master' into toplex_map

3 files changed, 92 insertions, 28 deletions

diff --git a/src/Tangential_complex/doc/Intro_tangential_complex.h b/src/Tangential_complex/doc/Intro_tangential_complex.h
index f4fc8ac7..501f4a8b 100644
--- a/src/Tangential_complex/doc/Intro_tangential_complex.h
+++ b/src/Tangential_complex/doc/Intro_tangential_complex.h
@@ -35,9 +35,11 @@ namespace tangential_complex {
 
 \section tangentialdefinition Definition
 
-A Tangential Delaunay complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
+A Tangential Delaunay complex is a
+<a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
 designed to reconstruct a \f$k\f$-dimensional smooth manifold embedded in \f$d\f$-dimensional Euclidean space. 
-The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of "small" intrinsic dimension.
+The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of
+"small" intrinsic dimension.
 The running time depends only linearly on the extrinsic dimension \f$ d \f$
 and exponentially on the intrinsic dimension \f$ k \f$.
 
@@ -46,17 +48,19 @@ An extensive description of the Tangential complex can be found in \cite tangent
 \subsection whatisthetc What is a Tangential Complex?
 
 Let us start with the description of the Tangential complex of a simple example, with \f$ k=1 \f$ and \f$ d=2 \f$.
-The input data is 4 points \f$ P \f$ located on a curve embedded in 2D.
+The point set \f$ \mathscr P \f$ is located on a closed curve embedded in 2D.
+Only 4 points will be displayed (more are required for PCA) to simplify the figures.
 \image html "tc_example_01.png" "The input"
-For each point \f$ p \f$, estimate its tangent subspace \f$ T_p \f$ (e.g. using PCA).
+For each point \f$ P \f$, estimate its tangent subspace \f$ T_P \f$ using PCA.
 \image html "tc_example_02.png" "The estimated normals"
-Let us add the Voronoi diagram of the points in orange. For each point \f$ p \f$, construct its star in the Delaunay triangulation of \f$ P \f$ restricted to \f$ T_p \f$.
+Let us add the Voronoi diagram of the points in orange. For each point \f$ P \f$, construct its star in the Delaunay
+triangulation of \f$ \mathscr P \f$ restricted to \f$ T_P \f$.
 \image html "tc_example_03.png" "The Voronoi diagram"
 The Tangential Delaunay complex is the union of those stars.
 
 In practice, neither the ambient Voronoi diagram nor the ambient Delaunay triangulation is computed.
-Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only need the star of each point.
-More details can be found in \cite tangentialcomplex2014.
+Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only
+need the star of each point. More details can be found in \cite tangentialcomplex2014.
 
 \subsection inconsistencies Inconsistencies
 
@@ -65,7 +69,7 @@ An inconsistency occurs when a simplex is not in the star of all its vertices.
 
 Let us take the same example.
 \image html "tc_example_07_before.png" "Before"
-Let us slightly move the tangent subspace \f$ T_q \f$
+Let us slightly move the tangent subspace \f$ T_Q \f$
 \image html "tc_example_07_after.png" "After"
 Now, the star of \f$ Q \f$ contains \f$ QP \f$, but the star of \f$ P \f$ does not contain \f$ QP \f$. We have an inconsistency.
 \image html "tc_example_08.png" "After"
diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h
index d1c846cf..4a78127c 100644
--- a/src/Tangential_complex/include/gudhi/Tangential_complex.h
+++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h
@@ -30,6 +30,7 @@
 #include <gudhi/console_color.h>
 #include <gudhi/Clock.h>
 #include <gudhi/Simplex_tree.h>
+#include <gudhi/Debug_utils.h>
 
 #include <CGAL/Default.h>
 #include <CGAL/Dimension.h>
@@ -321,7 +322,11 @@ class Tangential_complex {
     for (std::size_t i = 0; i < m_points.size(); ++i) m_are_tangent_spaces_computed[i] = true;
   }
 
-  /// Computes the tangential complex.
+  /** \brief Computes the tangential complex.
+   *  \exception std::invalid_argument In debug mode, if the computed star dimension is too low. Try to set a bigger
+   *  maximal edge length value with `Tangential_complex::set_max_squared_edge_length` if
+   *  this happens.
+   */
   void compute_tangential_complex() {
 #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS
     std::cerr << red << "WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. "
@@ -356,7 +361,7 @@ class Tangential_complex {
       tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*this));
     } else {
 #endif  // GUDHI_USE_TBB
-      // Sequential
+        // Sequential
       for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
 #ifdef GUDHI_USE_TBB
     }
@@ -629,12 +634,12 @@ class Tangential_complex {
         // Don't export infinite cells
         if (!export_infinite_simplices && is_infinite(c)) continue;
 
-        if (!export_inconsistent_simplices && !is_simplex_consistent(c)) continue;
-
         if (static_cast<int>(c.size()) > max_dim) max_dim = static_cast<int>(c.size());
         // Add the missing center vertex
         c.insert(idx);
 
+        if (!export_inconsistent_simplices && !is_simplex_consistent(c)) continue;
+
         // Try to insert the simplex
         bool inserted = tree.insert_simplex_and_subfaces(c).second;
 
@@ -689,6 +694,10 @@ class Tangential_complex {
         // Don't export infinite cells
         if (!export_infinite_simplices && is_infinite(c)) continue;
 
+        if (static_cast<int>(c.size()) > max_dim) max_dim = static_cast<int>(c.size());
+        // Add the missing center vertex
+        c.insert(idx);
+
         if (!export_inconsistent_simplices && !is_simplex_consistent(c)) continue;
 
         // Unusual simplex dim?
@@ -701,10 +710,6 @@ class Tangential_complex {
           check_lower_and_higher_dim_simplices = 1;
         }
 
-        if (static_cast<int>(c.size()) > max_dim) max_dim = static_cast<int>(c.size());
-        // Add the missing center vertex
-        c.insert(idx);
-
         // Try to insert the simplex
         bool added = complex.add_simplex(c, check_lower_and_higher_dim_simplices == 1);
 
@@ -800,7 +805,7 @@ class Tangential_complex {
       tbb::parallel_for(tbb::blocked_range<size_t>(0, m_points.size()), Compute_tangent_triangulation(*this));
     } else {
 #endif  // GUDHI_USE_TBB
-      // Sequential
+        // Sequential
       for (std::size_t i = 0; i < m_points.size(); ++i) compute_tangent_triangulation(i);
 #ifdef GUDHI_USE_TBB
     }
@@ -835,7 +840,7 @@ class Tangential_complex {
                         Refresh_tangent_triangulation(*this, updated_pts_ds));
     } else {
 #endif  // GUDHI_USE_TBB
-      // Sequential
+        // Sequential
       for (std::size_t i = 0; i < m_points.size(); ++i) refresh_tangent_triangulation(i, updated_pts_ds);
 #ifdef GUDHI_USE_TBB
     }
@@ -983,10 +988,9 @@ class Tangential_complex {
     // of the sphere "star sphere" centered at "center_vertex"
     // and which contains all the
     // circumspheres of the star of "center_vertex"
-    boost::optional<FT> squared_star_sphere_radius_plus_margin = boost::make_optional(false, FT());
-    // This is the strange way boost is recommending to get rid of "may be used uninitialized in this function".
-    // Former code was :
-    // boost::optional<FT> squared_star_sphere_radius_plus_margin;
+    // If th the m_max_squared_edge_length is set the maximal radius of the "star sphere"
+    // is at most square root of m_max_squared_edge_length
+    boost::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length;
 
     // Insert points until we find a point which is outside "star sphere"
     for (auto nn_it = ins_range.begin(); nn_it != ins_range.end(); ++nn_it) {
@@ -999,10 +1003,16 @@ class Tangential_complex {
         Point neighbor_pt;
         FT neighbor_weight;
         compute_perturbed_weighted_point(neighbor_point_idx, neighbor_pt, neighbor_weight);
-
+        GUDHI_CHECK(!m_max_squared_edge_length ||
+                    squared_star_sphere_radius_plus_margin.value() <= m_max_squared_edge_length.value(),
+                    std::invalid_argument("Tangential_complex::compute_star - set a bigger value with set_max_squared_edge_length."));
         if (squared_star_sphere_radius_plus_margin &&
-            k_sqdist(center_pt, neighbor_pt) > *squared_star_sphere_radius_plus_margin)
+            k_sqdist(center_pt, neighbor_pt) > squared_star_sphere_radius_plus_margin.value()) {
+          GUDHI_CHECK(triangulation.current_dimension() >= tangent_space_dim,
+                      std::invalid_argument("Tangential_complex::compute_star - Dimension of the star is only " + \
+                                            std::to_string(triangulation.current_dimension())));
           break;
+        }
 
         Tr_point proj_pt = project_point_and_compute_weight(neighbor_pt, neighbor_weight, tsb, local_tr_traits);
 
@@ -1044,7 +1054,7 @@ class Tangential_complex {
                 FT sq_power_sphere_diam = 4 * point_weight(c);
 
                 if (!squared_star_sphere_radius_plus_margin ||
-                    sq_power_sphere_diam > *squared_star_sphere_radius_plus_margin) {
+                    sq_power_sphere_diam > squared_star_sphere_radius_plus_margin.value()) {
                   squared_star_sphere_radius_plus_margin = sq_power_sphere_diam;
                 }
               }
@@ -1055,12 +1065,22 @@ class Tangential_complex {
             if (squared_star_sphere_radius_plus_margin) {
               // "2*m_last_max_perturb" because both points can be perturbed
               squared_star_sphere_radius_plus_margin =
-                  CGAL::square(std::sqrt(*squared_star_sphere_radius_plus_margin) + 2 * m_last_max_perturb);
+                  CGAL::square(std::sqrt(squared_star_sphere_radius_plus_margin.value()) + 2 * m_last_max_perturb);
+
+              // Reduce the square radius to  m_max_squared_edge_length if necessary
+              if (m_max_squared_edge_length && squared_star_sphere_radius_plus_margin.value() > m_max_squared_edge_length.value()) {
+                squared_star_sphere_radius_plus_margin = m_max_squared_edge_length.value();
+              }
 
               // Save it in `m_squared_star_spheres_radii_incl_margin`
-              m_squared_star_spheres_radii_incl_margin[i] = *squared_star_sphere_radius_plus_margin;
+              m_squared_star_spheres_radii_incl_margin[i] = squared_star_sphere_radius_plus_margin.value();
             } else {
-              m_squared_star_spheres_radii_incl_margin[i] = FT(-1);
+              if (m_max_squared_edge_length) {
+                squared_star_sphere_radius_plus_margin = m_max_squared_edge_length.value();
+                m_squared_star_spheres_radii_incl_margin[i] = m_max_squared_edge_length.value();
+              } else {
+                m_squared_star_spheres_radii_incl_margin[i] = FT(-1);
+              }
             }
           }
         }
@@ -1968,6 +1988,15 @@ class Tangential_complex {
     return os;
   }
 
+  /** \brief Sets the maximal possible squared edge length for the edges in the triangulations.
+   *
+   * @param[in] max_squared_edge_length Maximal possible squared edge length.
+   *
+   * If the maximal edge length value is too low `Tangential_complex::compute_tangential_complex` will throw an
+   * exception in debug mode.
+   */
+  void set_max_squared_edge_length(FT max_squared_edge_length) { m_max_squared_edge_length = max_squared_edge_length; }
+
  private:
   const K m_k;
   const int m_intrinsic_dim;
@@ -1993,6 +2022,7 @@ class Tangential_complex {
   // and their center vertex
   Stars_container m_stars;
   std::vector<FT> m_squared_star_spheres_radii_incl_margin;
+  boost::optional<FT> m_max_squared_edge_length;
 
 #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
   Points m_points_for_tse;
diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp
index 4e2d4f65..103b8b30 100644
--- a/src/Tangential_complex/test/test_tangential_complex.cpp
+++ b/src/Tangential_complex/test/test_tangential_complex.cpp
@@ -126,3 +126,33 @@ BOOST_AUTO_TEST_CASE(test_mini_tangential) {
   BOOST_CHECK(stree.num_vertices() == 4);
   BOOST_CHECK(stree.num_simplices() == 6);
 }
+
+#ifdef GUDHI_DEBUG
+BOOST_AUTO_TEST_CASE(test_basic_example_throw) {
+  typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
+  typedef Kernel::FT FT;
+  typedef Kernel::Point_d Point;
+  typedef Kernel::Vector_d Vector;
+  typedef tc::Tangential_complex<Kernel, CGAL::Dynamic_dimension_tag,CGAL::Parallel_tag> TC;
+
+  const int INTRINSIC_DIM = 2;
+  const int AMBIENT_DIM = 3;
+  const int NUM_POINTS = 1000;
+
+  Kernel k;
+
+  // Generate points on a 2-sphere
+  CGAL::Random_points_on_sphere_d<Point> generator(AMBIENT_DIM, 3.);
+  std::vector<Point> points;
+  points.reserve(NUM_POINTS);
+  for (int i = 0; i < NUM_POINTS; ++i)
+    points.push_back(*generator++);
+
+  // Compute the TC
+  TC tc(points, INTRINSIC_DIM, k);
+  tc.set_max_squared_edge_length(0.01);
+  std::cout << "test_basic_example_throw - set_max_squared_edge_length(0.01) to make GUDHI_CHECK fail" << std::endl;
+  BOOST_CHECK_THROW(tc.compute_tangential_complex(), std::invalid_argument);
+
+}
+#endif