Square root filtration values for alpha and delaunay-cech complex

src/Alpha_complex/include/gudhi/Alpha_complex.h

@@ -6,6 +6,7 @@
  *
  *    Modification(s):
  *      - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3
+ *      - 2024/10 Vincent Rouvreau: Add square root filtration values interface
  *      - YYYY/MM Author: Description of the modification
  */
 
@@ -65,14 +66,14 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val
 /**
  * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h
  * \brief Alpha complex data structure.
- * 
+ *
  * \ingroup alpha_complex
- * 
+ *
  * \details
- * The data structure is constructing a CGAL Delaunay triangulation (for more information on CGAL Delaunay 
+ * The data structure is constructing a CGAL Delaunay triangulation (for more information on CGAL Delaunay
  * triangulation, please refer to the corresponding chapter in page http://doc.cgal.org/latest/Triangulation/) from a
  * range of points or from an OFF file (cf. Points_off_reader).
- * 
+ *
  * Please refer to \ref alpha_complex for examples.
  *
  * The complex is a template class requiring an <a target="_blank"
@@ -84,7 +85,7 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val
  * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>
  * < <a target="_blank" href="http://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Dynamic__dimension__tag.html">
  * CGAL::Dynamic_dimension_tag </a> >
- * 
+ *
  * \remark
  * - When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex.
  * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the
@@ -161,10 +162,10 @@ class Alpha_complex {
 
  public:
   /** \brief Alpha_complex constructor from an OFF file name.
-   * 
-   * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize 
+   *
+   * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize
    * the Alpha_complex.
-   * 
+   *
    * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
    *
    * @param[in] off_file_name OFF file [path and] name.
@@ -183,9 +184,9 @@ class Alpha_complex {
    *
    * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points,
    * weighted hidden point, ...).
-   * 
+   *
    * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d.
-   * 
+   *
    * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
    * Kernel::Point_d or Kernel::Weighted_point_d.
    */
@@ -198,11 +199,11 @@ class Alpha_complex {
    *
    * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points,
    * weighted hidden point, ...).
-   * 
+   *
    * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d.
-   * 
+   *
    * @param[in] weights Range of points weights. Weights must be in Kernel::FT.
-   * 
+   *
    * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
    * Kernel::Point_d.
    */
@@ -354,10 +355,12 @@ class Alpha_complex {
    * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value
    * is not set.
    *
+   * \tparam output_squared_values If `true` (default value), it assigns to each simplex a filtration value equal to
+   * \f$ \alpha^2 \f$, or to \f$ \alpha \f$ when `output_squared_values` is `false`. cf. \ref createcomplexalgorithm
    * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept.
-   * 
+   *
    * @param[in] complex SimplicialComplexForAlpha to be created.
-   * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very
+   * @param[in] max_alpha_square maximum for \f$ \alpha^2 \f$ value. Default value is +\f$\infty\f$, and there is very
    * little point using anything else since it does not save time. Useless if `default_filtration_value` is set to
    * `true`.
    * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank"
@@ -366,14 +369,18 @@ class Alpha_complex {
    * (will be set to `NaN`).
    * Default value is `false` (which means compute the filtration values).
    *
-   * @return true if creation succeeds, false otherwise.
-   * 
+   * @return `true` if creation succeeds, `false` otherwise.
+   *
    * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0.
    * @pre The simplicial complex must be empty (no vertices)
-   * 
+   *
    * Initialization can be launched once.
+   *
+   * \note Weighted Alpha complex can have negative filtration values. If `output_squared_values` is set to `false`,
+   * filtration values will be `NaN` in this case.
    */
-  template <typename SimplicialComplexForAlpha,
+  template <bool output_squared_values = true,
+            typename SimplicialComplexForAlpha,
             typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value>
   bool create_complex(SimplicialComplexForAlpha& complex,
                       Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
@@ -389,6 +396,9 @@ class Alpha_complex {
     using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle;
     using Vector_vertex = std::vector<Vertex_handle>;
 
+    // For users to be able to define their own sqrt function on their desired Filtration_value type
+    using std::sqrt;
+
     if (triangulation_ == nullptr) {
       std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
       return false;  // ----- >>
@@ -438,7 +448,6 @@ class Alpha_complex {
       }
     }
     // --------------------------------------------------------------------------------------------
-
     if (!default_filtration_value) {
       CGAL::NT_converter<FT, Filtration_value> cgal_converter;
       // --------------------------------------------------------------------------------------------
@@ -458,6 +467,9 @@ class Alpha_complex {
                 if(exact) CGAL::exact(sqrad);
 #endif
                 alpha_complex_filtration = cgal_converter(sqrad);
+                if constexpr (!output_squared_values) {
+                  alpha_complex_filtration = sqrt(alpha_complex_filtration);
+                }
               }
               complex.assign_filtration(f_simplex, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
@@ -473,14 +485,18 @@ class Alpha_complex {
         cache_.clear();
       }
       // --------------------------------------------------------------------------------------------
-  
+
       // --------------------------------------------------------------------------------------------
       if (!exact)
         // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
         // Only in not exact version, cf. https://github.com/GUDHI/gudhi-devel/issues/57
         complex.make_filtration_non_decreasing();
       // Remove all simplices that have a filtration value greater than max_alpha_square
-      complex.prune_above_filtration(max_alpha_square);
+      if constexpr (!output_squared_values) {
+        complex.prune_above_filtration(sqrt(max_alpha_square));
+      } else {
+        complex.prune_above_filtration(max_alpha_square);
+      }
       // --------------------------------------------------------------------------------------------
     }
     return true;

src/Alpha_complex/test/Delaunay_complex_unit_test.h

@@ -24,7 +24,7 @@ using Simplex_handle = Simplex_tree::Simplex_handle;
 
 template<class CGAL_kernel>
 void compare_delaunay_complex_simplices() {
-  std::cout << "*****************************************************************************************************";
+  std::clog << "*****************************************************************************************************\n";
   using Point = typename CGAL_kernel::Point_d;
   std::vector<Point> points;
   // 50 points on a 4-sphere
@@ -40,10 +40,25 @@ void compare_delaunay_complex_simplices() {
   Simplex_tree stree_from_delaunay_complex;
   BOOST_CHECK(alpha_complex.create_complex(stree_from_delaunay_complex, 0., false, true));
 
-  // Check all the simplices from alpha complex are in the Delaunay complex
+  std::clog << "Check all the simplices from alpha complex are in the Delaunay complex\n";
   for (auto f_simplex : stree_from_alpha_complex.complex_simplex_range()) {
     Simplex_handle sh = stree_from_delaunay_complex.find(stree_from_alpha_complex.simplex_vertex_range(f_simplex));
     BOOST_CHECK(std::isnan(stree_from_delaunay_complex.filtration(sh)));
     BOOST_CHECK(sh != stree_from_delaunay_complex.null_simplex());
   }
+
+  std::clog << "Alpha complex with square root filtration values\n";
+  Simplex_tree stree_from_alpha_sqrt;
+  // set output_squared_values to false
+  BOOST_CHECK(alpha_complex.template create_complex<false>(stree_from_alpha_sqrt));
+
+  std::clog << "Check simplices from alpha complex filtration values when output_squared_values is true\n";
+  // Check all the simplices from alpha complex are in the Delaunay complex
+  for (auto f_simplex : stree_from_alpha_complex.complex_simplex_range()) {
+    Simplex_handle sh = stree_from_alpha_sqrt.find(stree_from_alpha_complex.simplex_vertex_range(f_simplex));
+    BOOST_CHECK(sh != stree_from_alpha_sqrt.null_simplex());
+    GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree_from_alpha_sqrt.filtration(sh),
+                                    std::sqrt(stree_from_alpha_complex.filtration(f_simplex)));
+  }
+
 }

src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp

@@ -124,4 +124,48 @@ BOOST_AUTO_TEST_CASE(Weighted_alpha_complex_3d_comparison) {
     }
     ++dD_itr;
   }
-}
+}
+
+BOOST_AUTO_TEST_CASE(Is_weighted_alpha_complex_nan) {
+  using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<3> >;
+  using Bare_point = Kernel::Point_d;
+  using Weighted_point = Kernel::Weighted_point_d;
+  using Vector_of_points = std::vector<Weighted_point>;
+
+  Vector_of_points points;
+  points.emplace_back(Bare_point(1, -1, -1), 4.);
+  points.emplace_back(Bare_point(-1, 1, -1), 4.);
+  points.emplace_back(Bare_point(-1, -1, 1), 4.);
+  points.emplace_back(Bare_point(1, 1, 1),   4.);
+  points.emplace_back(Bare_point(2, 2, 2),   1.);
+
+  Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_weighted_points(points);
+
+  std::clog << "Weighted alpha complex\n";
+  Gudhi::Simplex_tree<> stree;
+  if (alpha_complex_from_weighted_points.create_complex(stree)) {
+    for (auto f_simplex : stree.filtration_simplex_range()) {
+      std::clog << "   ( ";
+      for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
+        std::clog << vertex << " ";
+      }
+      std::clog << ") -> " << "[" << stree.filtration(f_simplex) << "]\n";
+
+      BOOST_CHECK(!std::isnan(stree.filtration(f_simplex)));
+    }
+  }
+  std::clog << "Weighted alpha complex with output_squared_values\n";
+  Gudhi::Simplex_tree<> stree_sqrt;
+  // set output_squared_values to false means that negative squared values will be NaN
+  if (alpha_complex_from_weighted_points.create_complex<false>(stree_sqrt)) {
+    for (auto f_simplex : stree_sqrt.filtration_simplex_range()) {
+        std::clog << "   ( ";
+        for (auto vertex : stree_sqrt.simplex_vertex_range(f_simplex)) {
+          std::clog << vertex << " ";
+        }
+        std::clog << ") -> " << "[" << stree_sqrt.filtration(f_simplex) << "]\n";
+
+        BOOST_CHECK(std::isnan(stree_sqrt.filtration(f_simplex)));
+    }
+  }
+}

src/Alpha_complex/utilities/CMakeLists.txt

@@ -7,6 +7,13 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options
       "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f")
   add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
       "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e")
+  # Same with sqrt filtration values - "-s"
+  add_test(NAME Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "safe_sqrt.pers")
+  add_test(NAME Alpha_complex_utilities_fast_sqrt_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "fast_sqrt.pers" "-f")
+  add_test(NAME Alpha_complex_utilities_exact_sqrt_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-s" "-p" "2" "-m" "0.45" "-o" "exact_sqrt.pers" "-e")
   if (DIFF_PATH)
     add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH}
              "exact.pers" "safe.pers")
@@ -17,6 +24,17 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen AND TARGET Boost::program_options
              "fast.pers" "safe.pers")
     set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex PROPERTIES DEPENDS
              "Alpha_complex_utilities_fast_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence")
+
+    # Same with sqrt filtration values - "-s"
+    add_test(Alpha_complex_utilities_diff_exact_sqrt_alpha_complex ${DIFF_PATH}
+             "exact_sqrt.pers" "safe_sqrt.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_exact_sqrt_alpha_complex PROPERTIES DEPENDS
+             "Alpha_complex_utilities_exact_sqrt_alpha_complex_persistence;Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence")
+
+    add_test(Alpha_complex_utilities_diff_fast_sqrt_alpha_complex ${DIFF_PATH}
+             "fast_sqrt.pers" "safe_sqrt.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_fast_sqrt_alpha_complex PROPERTIES DEPENDS
+             "Alpha_complex_utilities_fast_sqrt_alpha_complex_persistence;Alpha_complex_utilities_safe_sqrt_alpha_complex_persistence")
   endif()
 
   install(TARGETS alpha_complex_persistence DESTINATION bin)
@@ -62,6 +80,6 @@ if (TARGET CGAL::CGAL AND TARGET Boost::program_options)
       "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights"
       "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt"
       "-p" "2" "-m" "0" "-e")
-  
+
   install(TARGETS alpha_complex_3d_persistence DESTINATION bin)
 endif()