Use Epeck_d for examples. Add an Epick_d (fast version) example and test it

5 files changed, 94 insertions, 12 deletions

diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
index d411e90a..6b5f1a74 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
@@ -2,8 +2,6 @@
 // to construct a simplex_tree from alpha complex
 #include <gudhi/Simplex_tree.h>
 
-#include <CGAL/Epick_d.h>
-
 #include <iostream>
 #include <string>
 
@@ -23,8 +21,7 @@ int main(int argc, char **argv) {
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from an OFF file
   // ----------------------------------------------------------------------------
-  using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
-  Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name);
+  Gudhi::alpha_complex::Alpha_complex<> alpha_complex_from_file(off_file_name);
 
   std::streambuf* streambufffer;
   std::ofstream ouput_file_stream;
diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
index 981aa470..6526ca3a 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
@@ -2,12 +2,13 @@
 // to construct a simplex_tree from alpha complex
 #include <gudhi/Simplex_tree.h>
 
-#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
 
 #include <iostream>
 #include <vector>
 
-using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >;
+// Explicit dimension 2 Epeck_d kernel
+using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<2> >;
 using Point = Kernel::Point_d;
 using Vector_of_points = std::vector<Point>;
 
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index b069b443..b0337934 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -5,9 +5,12 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
   target_link_libraries(Alpha_complex_example_from_points  ${CGAL_LIBRARY})
   add_executable ( Alpha_complex_example_from_off Alpha_complex_from_off.cpp )
   target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY})
+  add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp )
+  target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY})
   if (TBB_FOUND)
     target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES})
     target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES})
+    target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES})
   endif()
 
   add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>)
@@ -16,7 +19,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
       "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt")
   add_test(NAME Alpha_complex_example_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off>
       "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt")
-  if (DIFF_PATH)
+
+  add_test(NAME Alpha_complex_example_fast_from_off_60 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off>
+      "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt")
+  add_test(NAME Alpha_complex_example_fast_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off>
+      "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt")
+
+if (DIFF_PATH)
     # Do not forget to copy test results files in current binary dir
     file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
     file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
@@ -25,6 +34,11 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
         ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
     add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH}
         ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
+
+    add_test(Alpha_complex_example_fast_from_off_60_diff_files ${DIFF_PATH}
+        ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
+    add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH}
+        ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
   endif()
 
   add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp )
diff --git a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
new file mode 100644
index 00000000..d6a39a76
--- /dev/null
+++ b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
@@ -0,0 +1,65 @@
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <CGAL/Epick_d.h>
+
+#include <iostream>
+#include <string>
+
+void usage(int nbArgs, char * const progName) {
+  std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
+  std::cerr << "Usage: " << progName << " filename.off alpha_square_max_value [ouput_file.txt]\n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0\n";
+  exit(-1);  // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
+
+  std::string off_file_name {argv[1]};
+  double alpha_square_max_value {atof(argv[2])};
+
+  // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d)
+  // (i.e. when the points are on a grid)
+  using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from an OFF file
+  // ----------------------------------------------------------------------------
+  Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_name);
+
+  std::streambuf* streambufffer;
+  std::ofstream ouput_file_stream;
+
+  if (argc == 4) {
+    ouput_file_stream.open(std::string(argv[3]));
+    streambufffer = ouput_file_stream.rdbuf();
+  } else {
+    streambufffer = std::cout.rdbuf();
+  }
+
+  Gudhi::Simplex_tree<> simplex;
+  if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
+    std::ostream output_stream(streambufffer);
+
+    // ----------------------------------------------------------------------------
+    // Display information about the alpha complex
+    // ----------------------------------------------------------------------------
+    output_stream << "Alpha complex is of dimension " << simplex.dimension() <<
+        " - " << simplex.num_simplices() << " simplices - " <<
+        simplex.num_vertices() << " vertices." << std::endl;
+
+    output_stream << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" <<
+        std::endl;
+    for (auto f_simplex : simplex.filtration_simplex_range()) {
+      output_stream << "   ( ";
+      for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
+        output_stream << vertex << " ";
+      }
+      output_stream << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
+      output_stream << std::endl;
+    }
+  }
+  ouput_file_stream.close();
+  return 0;
+}
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 27a4ba04..4e0e5159 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -24,7 +24,6 @@
 #include <CGAL/Epick_d.h>  // For FAST version
 #include <CGAL/Spatial_sort_traits_adapter_d.h>
 #include <CGAL/property_map.h>  // for CGAL::Identity_property_map
-#include <CGAL/NT_converter.h>
 #include <CGAL/version.h>  // for CGAL_VERSION_NR
 
 #include <Eigen/src/Core/util/Macros.h>  // for EIGEN_VERSION_AT_LEAST
@@ -249,6 +248,8 @@ class Alpha_complex {
    * @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
    * little point using anything else since it does not save time.
+   * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank"
+   * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>.
    * 
    * @return true if creation succeeds, false otherwise.
    * 
@@ -260,7 +261,8 @@ class Alpha_complex {
   template <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()) {
+                      Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
+                      bool exact = false) {
     // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
     typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
     typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle;
@@ -336,9 +338,12 @@ class Alpha_complex {
             if (f_simplex_dim > 0) {
               // squared_radius function initialization
               Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
-              CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
-
-              alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end()));
+              
+              if (exact) {
+                alpha_complex_filtration = CGAL::to_double(CGAL::exact(squared_radius(pointVector.begin(), pointVector.end())));
+              } else {
+                alpha_complex_filtration = CGAL::to_double(squared_radius(pointVector.begin(), pointVector.end()));
+              }
             }
             complex.assign_filtration(f_simplex, alpha_complex_filtration);
 #ifdef DEBUG_TRACES