Use unique_ptr to store cgal alpha_shapes_3. Required for not objects_ to be freed.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/alpha_complex_3d_module_vincent@3626 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 6dc9bbe25757a1a7bef2178ef63876c20d3a38b1

5 files changed, 317 insertions, 327 deletions

diff --git a/src/Alpha_complex/example/traits_test.cpp b/src/Alpha_complex/example/traits_test.cpp
index 99f1543c..3be62ad3 100644
--- a/src/Alpha_complex/example/traits_test.cpp
+++ b/src/Alpha_complex/example/traits_test.cpp
@@ -1,4 +1,5 @@
 #include <gudhi/Alpha_complex_3d.h>
+#include <gudhi/Simplex_tree.h>
 
 #include <iostream>
 #include <string>
@@ -15,31 +16,174 @@ void usage(int nbArgs, char * const progName) {
 int main(int argc, char **argv) {
   //if ((argc != 1) && (argc != 2)) usage(argc, (argv[0] - 1));
 
+  std::cout << "Alpha complex 3d" << std::endl;
   using Alpha_shapes_3d = Gudhi::alpha_complex::Alpha_shapes_3d;
   std::vector<Alpha_shapes_3d::Point_3> points;
-  points.push_back(Alpha_shapes_3d::Point_3(1., 2., 3.));
-  points.push_back(Alpha_shapes_3d::Point_3(6., 5., 4.));
+  points.push_back(Alpha_shapes_3d::Point_3(0.0, 0.0, 0.0));
+  points.push_back(Alpha_shapes_3d::Point_3(0.0, 0.0, 0.2));
+  points.push_back(Alpha_shapes_3d::Point_3(0.2, 0.0, 0.2));
+  points.push_back(Alpha_shapes_3d::Point_3(0.6, 0.6, 0.0));
+  points.push_back(Alpha_shapes_3d::Point_3(0.8, 0.8, 0.2));
+  points.push_back(Alpha_shapes_3d::Point_3(0.2, 0.8, 0.6));
 
   Gudhi::alpha_complex::Alpha_complex_3d<Alpha_shapes_3d> alpha_complex(points);
 
+  Gudhi::Simplex_tree<> stree;
+  alpha_complex.create_complex(stree);
+
+  std::cout << "Weighted alpha complex 3d" << std::endl;
   using Weighted_alpha_shapes_3d = Gudhi::alpha_complex::Weighted_alpha_shapes_3d;
   std::vector<Weighted_alpha_shapes_3d::Point_3> w_points;
-  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(1., 2., 3.));
-  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(6., 5., 4.));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.0, 0.0, 0.0));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.0, 0.0, 0.2));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.2, 0.0, 0.2));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.6, 0.6, 0.0));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.8, 0.8, 0.2));
+  w_points.push_back(Weighted_alpha_shapes_3d::Point_3(0.2, 0.8, 0.6));
+
+  std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001};
 
-  std::vector<double> weights = {0.01, 0.005};
+  Gudhi::alpha_complex::Alpha_complex_3d<Weighted_alpha_shapes_3d> weighted_alpha_complex(w_points, weights);
 
-  Gudhi::alpha_complex::Alpha_complex_3d<Weighted_alpha_shapes_3d> weighted_alpha_complex(points, weights);
+  Gudhi::Simplex_tree<> w_stree;
+  weighted_alpha_complex.create_complex(w_stree);
 
+  std::cout << "Periodic alpha complex 3d" << std::endl;
   using Periodic_alpha_shapes_3d = Gudhi::alpha_complex::Periodic_alpha_shapes_3d;
   std::vector<Periodic_alpha_shapes_3d::Point_3> p_points;
-  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(1., 2., 3.));
-  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(6., 5., 4.));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.0, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.0, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.0, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.0, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.2, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.2, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.2, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.2, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.2, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.4, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.4, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.4, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.4, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.4, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.6, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.6, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.6, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.6, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.6, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.8, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.8, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.8, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.8, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.0, 0.8, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.0, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.0, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.0, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.0, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.2, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.2, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.2, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.2, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.2, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.4, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.4, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.4, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.4, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.4, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.6, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.6, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.6, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.6, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.6, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.8, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.8, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.8, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.8, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.2, 0.8, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.0, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.0, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.0, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.0, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.2, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.2, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.2, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.2, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.2, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.4, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.4, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.4, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.4, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.4, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.6, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.6, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.6, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.6, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.6, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.8, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.8, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.8, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.8, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.4, 0.8, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.1, 0.0));  //
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.2, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.2, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.2, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.2, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.2, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.4, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.4, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.4, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.4, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.4, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.6, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.6, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.6, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.6, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.6, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.8, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.8, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.8, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.8, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.6, 0.8, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.0, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.0, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.0, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.0, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.0, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.2, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.2, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.2, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.2, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.2, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.4, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.4, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.4, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.4, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.4, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.6, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.6, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.6, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.6, 0.6));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.6, 0.8));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.8, 0.0));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.8, 0.2));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.8, 0.4));
+  p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.8, 0.6));
+  //p_points.push_back(Periodic_alpha_shapes_3d::Point_3(0.8, 0.8, 0.8));
 
   Gudhi::alpha_complex::Alpha_complex_3d<Periodic_alpha_shapes_3d> periodic_alpha_complex(points,
                                                                                           0., 0., 0.,
                                                                                           1., 1., 1.);
 
+  std::cout << "Weighted periodic alpha complex 3d" << std::endl;
   using Weighted_periodic_alpha_shapes_3d = Gudhi::alpha_complex::Weighted_periodic_alpha_shapes_3d;
   std::vector<Weighted_periodic_alpha_shapes_3d::Point_3> wp_points;
   wp_points.push_back(Weighted_periodic_alpha_shapes_3d::Point_3(0.1, 0.2, 0.3));
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index 79f8593d..c8bc9c8e 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -37,6 +37,12 @@
 #include <CGAL/tuple.h>
 #include <CGAL/iterator.h>
 
+#include <iostream>
+#include <vector>
+#include <map>
+#include <stdexcept>
+#include <cstddef>
+#include <memory>  // for std::unique_ptr
 
 namespace Gudhi {
 
@@ -44,19 +50,17 @@ namespace alpha_complex {
 
 template<typename AlphaComplex3dOptions>
 class Alpha_complex_3d {
-private:
   using Alpha_shape_3 = typename AlphaComplex3dOptions::Alpha_shape_3;
-  // filtration with alpha values needed type definition
   using Alpha_value_type = typename Alpha_shape_3::FT;
-  using Object = CGAL::Object;
-  using Dispatch = CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
-      CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
-          std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+  using Dispatch =
+  CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, Alpha_value_type>,
+  CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<CGAL::Object> >,
+      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
+
   using Cell_handle = typename Alpha_shape_3::Cell_handle;
   using Facet = typename Alpha_shape_3::Facet;
-  using Edge_3 = typename Alpha_shape_3::Edge;
+  using Edge = typename Alpha_shape_3::Edge;
   using Alpha_vertex_handle = typename Alpha_shape_3::Vertex_handle;
-
 #if BOOST_VERSION >= 105400
   using Vertex_list = boost::container::static_vector<Alpha_vertex_handle, 4>;
 #else
@@ -64,6 +68,9 @@ private:
 #endif
 
 public:
+  using Point_3 = typename AlphaComplex3dOptions::Point_3;
+
+public:
   /** \brief Alpha_complex constructor from a list of points.
   *
   * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
@@ -75,17 +82,19 @@ public:
   */
   template<typename InputPointRange >
   Alpha_complex_3d(const InputPointRange& points) {
-    static_assert(!AlphaComplex3dOptions::weighted, "This constructor is not available for weighted versions of Alpha_complex_3d");
-    static_assert(!AlphaComplex3dOptions::periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");
-    std::cout << points[0] << std::endl;
-    Alpha_shape_3 alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL);
+    static_assert(!AlphaComplex3dOptions::weighted,
+                  "This constructor is not available for weighted versions of Alpha_complex_3d");
+    static_assert(!AlphaComplex3dOptions::periodic,
+                  "This constructor is not available for periodic versions of Alpha_complex_3d");
 
-    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                    std::back_inserter(the_alpha_values));
+    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(points.begin(), points.end(), 0,
+                                                        Alpha_shape_3::GENERAL));
+    Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, Alpha_value_type>(std::back_inserter(objects_),
+                                                                                std::back_inserter(alpha_values_));
 
-    alpha_shape_3.filtration_with_alpha_values(disp);
+    alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
 #ifdef DEBUG_TRACES
-    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+    std::cout << "filtration_with_alpha_values returns : " << objects_.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
 
   }
@@ -106,7 +115,7 @@ public:
     static_assert(!AlphaComplex3dOptions::periodic,
                   "This constructor is not available for periodic versions of Alpha_complex_3d");
     GUDHI_CHECK((weights.size() == points.size()),
-                std::invalid_argument("Alpha_complex_3d constructor with weights requires points number to be equal with points number"));
+                std::invalid_argument("Points number in range different from weights range number"));
 
     using Weighted_point_3 = typename AlphaComplex3dOptions::Weighted_point_3;
     std::vector<Weighted_point_3> weighted_points_3;
@@ -118,14 +127,17 @@ public:
       index++;
     }
 
-    Alpha_shape_3 alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL);
+    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(std::begin(weighted_points_3),
+                                                                          std::end(weighted_points_3),
+                                                                          0,
+                                                                          Alpha_shape_3::GENERAL));
 
-    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                    std::back_inserter(the_alpha_values));
+    Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, Alpha_value_type>(std::back_inserter(objects_),
+                                                                    std::back_inserter(alpha_values_));
 
-    alpha_shape_3.filtration_with_alpha_values(disp);
+    alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
 #ifdef DEBUG_TRACES
-    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+    std::cout << "filtration_with_alpha_values returns : " << objects_.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
   }
 
@@ -140,14 +152,16 @@ public:
 */
   template<typename InputPointRange>
   Alpha_complex_3d(const InputPointRange& points,
-                   Alpha_value_type x_min, Alpha_value_type y_min, Alpha_value_type z_min,
-                   Alpha_value_type x_max, Alpha_value_type y_max, Alpha_value_type z_max) {
+                   double x_min, double y_min, double z_min,
+                   double x_max, double y_max, double z_max) {
     static_assert(!AlphaComplex3dOptions::weighted,
                   "This constructor is not available for weighted versions of Alpha_complex_3d");
     static_assert(AlphaComplex3dOptions::periodic,
                   "This constructor is not available for non-periodic versions of Alpha_complex_3d");
     // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
-    GUDHI_CHECK((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min),
+    GUDHI_CHECK((x_max - x_min == y_max - y_min) ||
+                (x_max - x_min == z_max - z_min) ||
+                (z_max - z_min == y_max - y_min),
                 std::invalid_argument("The size of the cuboid in every directions is not the same."));
 
     using Periodic_delaunay_triangulation_3 = typename AlphaComplex3dOptions::Periodic_delaunay_triangulation_3;
@@ -157,23 +171,27 @@ public:
     // Heuristic for inserting large point sets (if pts is reasonably large)
     pdt.insert(std::begin(points), std::end(points), true);
     // As pdt won't be modified anymore switch to 1-sheeted cover if possible
-    GUDHI_CHECK(pdt.is_triangulation_in_1_sheet(),
-                std::invalid_argument("Uable to construct a triangulation within a single periodic domain."));
+    // GUDHI_CHECK(pdt.is_triangulation_in_1_sheet(),
+    //            std::invalid_argument("Uable to construct a triangulation within a single periodic domain."));
+    if (pdt.is_triangulation_in_1_sheet()) {
+      pdt.convert_to_1_sheeted_covering();
+    } else {
+      std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl;
+      exit(-1);
+    }
+    std::cout << "Periodic Delaunay computed." << std::endl;
 
     // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
     // Maybe need to set it to GENERAL mode
-    Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL);
-
-    // filtration with alpha values from alpha shape
-    std::vector<Object> the_objects;
-    std::vector<Alpha_value_type> the_alpha_values;
+    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0,
+                                                                          Alpha_shape_3::GENERAL));
 
-    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                    std::back_inserter(the_alpha_values));
+    Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, Alpha_value_type>(std::back_inserter(objects_),
+                                                                    std::back_inserter(alpha_values_));
 
-    as.filtration_with_alpha_values(disp);
+    alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
 #ifdef DEBUG_TRACES
-    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+    std::cout << "filtration_with_alpha_values returns : " << objects_.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
 
   }
@@ -196,12 +214,17 @@ public:
     static_assert(AlphaComplex3dOptions::periodic,
                   "This constructor is not available for non-periodic versions of Alpha_complex_3d");
     GUDHI_CHECK((weights.size() == points.size()),
-                std::invalid_argument("Alpha_complex_3d constructor with weights requires points number to be equal with points number"));
+                std::invalid_argument("Points number in range different from weights range number"));
     // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
-    GUDHI_CHECK((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min),
+    GUDHI_CHECK((x_max - x_min == y_max - y_min) ||
+                (x_max - x_min == z_max - z_min) ||
+                (z_max - z_min == y_max - y_min),
                 std::invalid_argument("The size of the cuboid in every directions is not the same."));
 
+#ifdef GUDHI_DEBUG
+    // Defined in GUDHI_DEBUG to avoid unused variable warning
     double maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
+#endif
 
     using Weighted_point_3 = typename AlphaComplex3dOptions::Weighted_point_3;
     std::vector<Weighted_point_3> weighted_points_3;
@@ -225,25 +248,25 @@ public:
     pdt.insert(std::begin(weighted_points_3), std::end(weighted_points_3), true);
     // As pdt won't be modified anymore switch to 1-sheeted cover if possible
     GUDHI_CHECK(pdt.is_triangulation_in_1_sheet(),
-                std::invalid_argument("Uable to construct a triangulation within a single periodic domain."));
+                std::invalid_argument("Unable to construct a triangulation within a single periodic domain."));
+
+    pdt.convert_to_1_sheeted_covering();
 
     // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
     // Maybe need to set it to GENERAL mode
-    Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL);
+    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0,
+                                                                          Alpha_shape_3::GENERAL));
 
-    // filtration with alpha values from alpha shape
-    std::vector<Object> the_objects;
-    std::vector<Alpha_value_type> the_alpha_values;
+    Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, Alpha_value_type>(std::back_inserter(objects_),
+                                                                    std::back_inserter(alpha_values_));
 
-    Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                    std::back_inserter(the_alpha_values));
-
-    as.filtration_with_alpha_values(disp);
+    alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
 #ifdef DEBUG_TRACES
-    std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
+    std::cout << "filtration_with_alpha_values returns : " << objects_.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
   }
 
+
   template <typename SimplicialComplexForAlpha3d>
   void create_complex(SimplicialComplexForAlpha3d& complex) {
     using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value;
@@ -251,68 +274,86 @@ public:
   }
 
   /** \brief Inserts all Delaunay triangulation into the simplicial complex.
-   * It also computes the filtration values accordingly to the \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$.
-   *
-   * @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.
-   */
+ * It also computes the filtration values accordingly to the \ref createcomplexalgorithm
+ *
+ * \tparam SimplicialComplexForAlpha3d must meet `SimplicialComplexForAlpha3d` concept.
+ *
+ * @param[in] complex SimplicialComplexForAlpha3d to be created.
+ * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$.
+ *
+ * @return true if creation succeeds, false otherwise.
+ *
+ * @pre The simplicial complex must be empty (no vertices)
+ *
+ * Initialization can be launched once.
+ */
   template <typename SimplicialComplexForAlpha3d>
-  void create_complex(SimplicialComplexForAlpha3d& complex,
+  bool create_complex(SimplicialComplexForAlpha3d& complex,
                       typename SimplicialComplexForAlpha3d::Filtration_value max_alpha_square) {
+    if (complex.num_vertices() > 0) {
+      std::cerr << "Alpha_complex_3d create_complex - complex is not empty\n";
+      return false;  // ----- >>
+    }
+
     using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value;
-    using Vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle;
+    using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle;
+    using Alpha_shape_simplex_tree_map = std::map<Alpha_vertex_handle,
+                                                  Complex_vertex_handle>;
+    using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>;
+
 #ifdef DEBUG_TRACES
-    Alpha_shape_3::size_type count_vertices = 0;
-    Alpha_shape_3::size_type count_edges = 0;
-    Alpha_shape_3::size_type count_facets = 0;
-    Alpha_shape_3::size_type count_cells = 0;
+    std::size_t count_vertices = 0;
+    std::size_t count_edges = 0;
+    std::size_t count_facets = 0;
+    std::size_t count_cells = 0;
 #endif  // DEBUG_TRACES
 
-    // Loop on objects vector
-    Vertex_list vertex_list;
-    SimplicialComplexForAlpha3d simplex_tree;
-    std::map<Alpha_vertex_handle, Vertex_handle> map_cgal_simplex_tree;
-    typename std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-    for (auto object_iterator : the_objects) {
+    Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
+    auto the_alpha_value_iterator = alpha_values_.begin();
+    for (auto object_iterator : objects_) {
+      Vertex_list vertex_list;
+
       // Retrieve Alpha shape vertex list from object
       if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-        vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
+        for (auto i = 0; i < 4; i++) {
 #ifdef DEBUG_TRACES
-        count_cells++;
+          std::cout << "from cell[" << i << "]=" << (*cell)->vertex(i)->point() << std::endl;
 #endif  // DEBUG_TRACES
+          vertex_list.push_back((*cell)->vertex(i));
+        }
+        count_cells++;
       } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
-        vertex_list = from_facet<Vertex_list, Facet>(*facet);
+          for (auto i = 0; i < 4; i++) {
+            if ((*facet).second != i) {
 #ifdef DEBUG_TRACES
-        count_facets++;
+              std::cout << "from facet=[" << i << "]" << (*facet).first->vertex(i)->point() << std::endl;
 #endif  // DEBUG_TRACES
-      } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-        vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
+              vertex_list.push_back((*facet).first->vertex(i));
+            }
+          }
+        count_facets++;
+      } else if (const Edge *edge = CGAL::object_cast<Edge>(&object_iterator)) {
+            for (auto i : {(*edge).second, (*edge).third}) {
 #ifdef DEBUG_TRACES
-        count_edges++;
+              std::cout << "from edge[" << i << "]=" << (*edge).first->vertex(i)->point() << std::endl;
 #endif  // DEBUG_TRACES
-      } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
-        vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-#ifdef DEBUG_TRACES
+              vertex_list.push_back((*edge).first->vertex(i));
+            }
+        count_edges++;
+      } else if (const Alpha_vertex_handle *vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) {
         count_vertices++;
+#ifdef DEBUG_TRACES
+              std::cout << "from vertex=" << (*vertex)->point() << std::endl;
 #endif  // DEBUG_TRACES
+              vertex_list.push_back((*vertex));
       }
       // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-      std::vector<Vertex_handle> the_simplex;
+      Simplex_tree_vector_vertex the_simplex;
       for (auto the_alpha_shape_vertex : vertex_list) {
-        typename std::map<Alpha_vertex_handle, Vertex_handle>::iterator the_map_iterator =
-            map_cgal_simplex_tree.find(the_alpha_shape_vertex);
+        auto the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
         if (the_map_iterator == map_cgal_simplex_tree.end()) {
           // alpha shape not found
-          Vertex_handle vertex = map_cgal_simplex_tree.size();
+          Complex_vertex_handle vertex = map_cgal_simplex_tree.size();
 #ifdef DEBUG_TRACES
           std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
 #endif  // DEBUG_TRACES
@@ -320,7 +361,7 @@ public:
           map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
         } else {
           // alpha shape found
-          Vertex_handle vertex = the_map_iterator->second;
+          Complex_vertex_handle vertex = the_map_iterator->second;
 #ifdef DEBUG_TRACES
           std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
 #endif  // DEBUG_TRACES
@@ -332,84 +373,25 @@ public:
 #ifdef DEBUG_TRACES
       std::cout << "filtration = " << filtr << std::endl;
 #endif  // DEBUG_TRACES
-      simplex_tree.insert_simplex(the_simplex, filtr);
-      GUDHI_CHECK(the_alpha_value_iterator != the_alpha_values.end(), "CGAL provided more simplices than values");
+      complex.insert_simplex(the_simplex, filtr);
+      GUDHI_CHECK(the_alpha_value_iterator != alpha_values_.end(), "CGAL provided more simplices than values");
       ++the_alpha_value_iterator;
     }
 
 #ifdef DEBUG_TRACES
-    std::cout << "vertices \t\t" << count_vertices << std::endl;
+    std::cout << "vertices \t" << count_vertices << std::endl;
     std::cout << "edges \t\t" << count_edges << std::endl;
     std::cout << "facets \t\t" << count_facets << std::endl;
     std::cout << "cells \t\t" << count_cells << std::endl;
-
-    std::cout << "Information of the Simplex Tree: " << std::endl;
-    std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-    std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-    std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
-
-#ifdef DEBUG_TRACES
-    std::cout << "Iterator on vertices: " << std::endl;
-    for (auto vertex : simplex_tree.complex_vertex_range()) {
-      std::cout << vertex << " ";
-    }
-#endif  // DEBUG_TRACES
-
-  }
-
-private:
-  template <class Vertex_list, class Cell_handle>
-  Vertex_list from_cell(const Cell_handle& ch) {
-    Vertex_list the_list;
-    for (auto i = 0; i < 4; i++) {
-#ifdef DEBUG_TRACES
-      std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-      the_list.push_back(ch->vertex(i));
-    }
-    return the_list;
-  }
-
-  template <class Vertex_list, class Facet>
-  Vertex_list from_facet(const Facet& fct) {
-    Vertex_list the_list;
-    for (auto i = 0; i < 4; i++) {
-      if (fct.second != i) {
-#ifdef DEBUG_TRACES
-        std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-        the_list.push_back(fct.first->vertex(i));
-      }
-    }
-    return the_list;
-  }
-
-  template <class Vertex_list, class Edge_3>
-  Vertex_list from_edge(const Edge_3& edg) {
-    Vertex_list the_list;
-    for (auto i : {edg.second, edg.third}) {
-#ifdef DEBUG_TRACES
-      std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-      the_list.push_back(edg.first->vertex(i));
-    }
-    return the_list;
-  }
-
-  template <class Vertex_list, class Vertex_handle>
-  Vertex_list from_vertex(const Vertex_handle& vh) {
-    Vertex_list the_list;
-#ifdef DEBUG_TRACES
-    std::cout << "from vertex=" << vh->point() << std::endl;
 #endif  // DEBUG_TRACES
-    the_list.push_back(vh);
-    return the_list;
+    return true;
   }
 
 private:
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
+  // Needs to store alpha_shape_3_ptr_ as objects_ and alpha_shape_3_ptr_ are freed with alpha_shape_3_ptr_
+  std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_;
+  std::vector<CGAL::Object> objects_;
+  std::vector<Alpha_value_type> alpha_values_;
 
 };
 
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h
index 9433afdc..3b1981ca 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d_options.h
@@ -42,7 +42,6 @@ namespace alpha_complex {
 
 class Alpha_shapes_3d {
 private:
-  // Alpha_shape_3 templates type definitions
   using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
   using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>;
   using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>;
@@ -55,7 +54,6 @@ public:
 
   static const bool weighted = false;
   static const bool periodic = false;
-
 };
 
 class Exact_alpha_shapes_3d {
@@ -98,24 +96,24 @@ public:
 
 class Periodic_alpha_shapes_3d {
 private:
-  using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-  using Periodic_kernel = CGAL::Periodic_3_Delaunay_triangulation_traits_3<Kernel>;
+  // Traits
+  using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+  using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>;
 // Vertex type
   using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
-  using Vb = CGAL::Triangulation_vertex_base_3<Periodic_kernel, DsVb>;
-  using AsVb = CGAL::Alpha_shape_vertex_base_3<Periodic_kernel, Vb>;
+  using Vb = CGAL::Triangulation_vertex_base_3<PK, DsVb>;
+  using AsVb = CGAL::Alpha_shape_vertex_base_3<PK, Vb>;
 // Cell type
   using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>;
-  using Cb = CGAL::Triangulation_cell_base_3<Periodic_kernel, DsCb>;
-  using AsCb = CGAL::Alpha_shape_cell_base_3<Periodic_kernel, Cb>;
+  using Cb = CGAL::Triangulation_cell_base_3<PK, DsCb>;
+  using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>;
   using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>;
 
 public:
-  using Periodic_delaunay_triangulation_3 = CGAL::Periodic_3_Delaunay_triangulation_3<Periodic_kernel, Tds>;
+  using Periodic_delaunay_triangulation_3 = CGAL::Periodic_3_Delaunay_triangulation_3<PK, Tds>;
   using Alpha_shape_3 = CGAL::Alpha_shape_3<Periodic_delaunay_triangulation_3>;
-  using Point_3 = Periodic_kernel::Point_3;
-  using Alpha_value_type = Alpha_shape_3::FT;
-  using Iso_cuboid_3 = Periodic_kernel::Iso_cuboid_3;
+  using Point_3 = PK::Point_3;
+  using Iso_cuboid_3 = PK::Iso_cuboid_3;
 
   static const bool weighted = false;
   static const bool periodic = true;
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index 7ace6064..c4d986c4 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -8,10 +8,15 @@ if(CGAL_FOUND)
   add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp)
   target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
 
+  add_executable(old_alpha_complex_3d_persistence old_alpha_complex_3d_persistence.cpp)
+  target_link_libraries(old_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
   if (TBB_FOUND)
     target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES})
     target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES})
     target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+
+    target_link_libraries(old_alpha_complex_3d_persistence ${TBB_LIBRARIES})
   endif(TBB_FOUND)
   
   add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
index 8cda0b70..669b6e08 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -20,73 +20,24 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <boost/version.hpp>
 #include <boost/program_options.hpp>
-#include <boost/variant.hpp>
-
-#if BOOST_VERSION >= 105400
-#include <boost/container/static_vector.hpp>
-#endif
 
+#include <gudhi/Alpha_complex_3d.h>
+#include <gudhi/Alpha_complex_3d_options.h>
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Persistent_cohomology.h>
 #include <gudhi/Points_3D_off_io.h>
 
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Delaunay_triangulation_3.h>
-#include <CGAL/Alpha_shape_3.h>
-#include <CGAL/Alpha_shape_cell_base_3.h>
-#include <CGAL/Alpha_shape_vertex_base_3.h>
-#include <CGAL/iterator.h>
-
-#include <fstream>
-#include <cmath>
+#include <iostream>
 #include <string>
-#include <tuple>
-#include <map>
-#include <utility>
-#include <vector>
-#include <cstdlib>
-
-#include "alpha_complex_3d_helper.h"
-
-// Alpha_shape_3 templates type definitions
-using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>;
-using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>;
-using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
-using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
-using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>;
-
-// From file type definition
-using Point_3 = Kernel::Point_3;
-
-// filtration with alpha values needed type definition
-using Alpha_value_type = Alpha_shape_3::FT;
-using Object = CGAL::Object;
-using Dispatch =
-    CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>,
-                                   CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >,
-                                                      std::back_insert_iterator<std::vector<Alpha_value_type> > > >;
-using Cell_handle = Alpha_shape_3::Cell_handle;
-using Facet = Alpha_shape_3::Facet;
-using Edge_3 = Alpha_shape_3::Edge;
-using Vertex_handle = Alpha_shape_3::Vertex_handle;
-
-#if BOOST_VERSION >= 105400
-using Vertex_list = boost::container::static_vector<Alpha_shape_3::Vertex_handle, 4>;
-#else
-using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
-#endif
+#include <limits>  // for numeric_limits
 
 // gudhi type definition
-using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = ST::Filtration_value;
-using Simplex_tree_vertex = ST::Vertex_handle;
-using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::Alpha_shapes_3d>;
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = Simplex_tree::Filtration_value;
 using Persistent_cohomology =
-    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp>;
 
 void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
                      int &coeff_field_characteristic, Filtration_value &min_persistence);
@@ -100,108 +51,18 @@ int main(int argc, char **argv) {
   program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence);
 
   // Read the OFF file (input file name given as parameter) and triangulate points
-  Gudhi::Points_3D_off_reader<Point_3> off_reader(off_file_points);
+  Gudhi::Points_3D_off_reader<Alpha_complex_3d::Point_3> off_reader(off_file_points);
   // Check the read operation was correct
   if (!off_reader.is_valid()) {
     std::cerr << "Unable to read file " << off_file_points << std::endl;
     exit(-1);
   }
 
-  // Retrieve the points
-  std::vector<Point_3> lp = off_reader.get_point_cloud();
-
-  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
-  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
-#ifdef DEBUG_TRACES
-  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
-#endif  // DEBUG_TRACES
-
-  // filtration with alpha values from alpha shape
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
-
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                  std::back_inserter(the_alpha_values));
-
-  as.filtration_with_alpha_values(disp);
-#ifdef DEBUG_TRACES
-  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif  // DEBUG_TRACES
-
-  Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges = 0;
-  Alpha_shape_3::size_type count_facets = 0;
-  Alpha_shape_3::size_type count_cells = 0;
-
-  // Loop on objects vector
-  Vertex_list vertex_list;
-  ST simplex_tree;
-  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
-  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  for (auto object_iterator : the_objects) {
-    // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
-      count_cells++;
-    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
-      vertex_list = from_facet<Vertex_list, Facet>(*facet);
-      count_facets++;
-    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
-      count_edges++;
-    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
-      count_vertices++;
-      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-    }
-    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-    Simplex_tree_vector_vertex the_simplex;
-    for (auto the_alpha_shape_vertex : vertex_list) {
-      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end()) {
-        // alpha shape not found
-        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex.push_back(vertex);
-        map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
-      } else {
-        // alpha shape found
-        Simplex_tree_vertex vertex = the_map_iterator->second;
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex.push_back(vertex);
-      }
-    }
-    // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
-#ifdef DEBUG_TRACES
-    std::cout << "filtration = " << filtr << std::endl;
-#endif  // DEBUG_TRACES
-    simplex_tree.insert_simplex(the_simplex, filtr);
-    GUDHI_CHECK(the_alpha_value_iterator != the_alpha_values.end(), "CGAL provided more simplices than values");
-    ++the_alpha_value_iterator;
-  }
-
-#ifdef DEBUG_TRACES
-  std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t" << count_edges << std::endl;
-  std::cout << "facets \t\t" << count_facets << std::endl;
-  std::cout << "cells \t\t" << count_cells << std::endl;
+  Alpha_complex_3d alpha_complex(off_reader.get_point_cloud());
 
-  std::cout << "Information of the Simplex Tree: " << std::endl;
-  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
+  Simplex_tree simplex_tree;
 
-#ifdef DEBUG_TRACES
-  std::cout << "Iterator on vertices: " << std::endl;
-  for (auto vertex : simplex_tree.complex_vertex_range()) {
-    std::cout << vertex << " ";
-  }
-#endif  // DEBUG_TRACES
+  alpha_complex.create_complex(simplex_tree);
 
   // Sort the simplices in the order of the filtration
   simplex_tree.initialize_filtration();