Merge last trunk modifications

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

Former-commit-id: 213c7b76c5f8c3248b24a2ac9250ed59034d8076

16 files changed, 1932 insertions, 86 deletions

diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
index 2b8bff94..a51df127 100644
--- a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
+++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Vincent Rouvreau
  *
- *    Copyright (C) 2016  INRIA
+ *    Copyright (C) 2016 Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
diff --git a/src/Alpha_complex/doc/COPYRIGHT b/src/Alpha_complex/doc/COPYRIGHT
index dbad2380..5f1d97cc 100644
--- a/src/Alpha_complex/doc/COPYRIGHT
+++ b/src/Alpha_complex/doc/COPYRIGHT
@@ -4,7 +4,7 @@ computational topology.
 
 Author(s):       Vincent Rouvreau
 
-Copyright (C) 2015  INRIA
+Copyright (C) 2015 Inria
 
 This program is free software: you can redistribute it and/or modify it under
 the terms of the GNU General Public License as published by the Free Software
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h
index cf1a946a..7a375c9f 100644
--- a/src/Alpha_complex/doc/Intro_alpha_complex.h
+++ b/src/Alpha_complex/doc/Intro_alpha_complex.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Vincent Rouvreau
  *
- *    Copyright (C) 2015  INRIA
+ *    Copyright (C) 2015 Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -31,7 +31,7 @@ namespace alpha_complex {
 /**  \defgroup alpha_complex Alpha complex
  * 
  * \author    Vincent Rouvreau
- * 
+ *
  * @{
  * 
  * \section definition Definition
@@ -57,9 +57,13 @@ namespace alpha_complex {
  * href="http://doc.cgal.org/latest/Kernel_d/index.html#Chapter_dD_Geometry_Kernel">dD Geometry Kernel</a>
  * \cite cgal:s-gkd-15b from CGAL as template parameter.
  * 
- * \remark When the simplicial complex is constructed with an infinite value of alpha, the complex is a Delaunay
+ * \remark
+ * - When the simplicial complex is constructed with an infinite value of alpha, the complex is a Delaunay
  * complex.
- * 
+ * - For people only interested in the topology of the \ref alpha_complex (for instance persistence),
+ * \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii.
+ * \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii.
+ *
  * \section pointsexample Example from points
  * 
  * This example builds the Delaunay triangulation from the given points in a 2D static kernel, and creates a
@@ -89,63 +93,29 @@ namespace alpha_complex {
  * \image html "alpha_complex_doc.png" "Simplicial complex structure construction example"
  *
  * \subsection filtrationcomputation Filtration value computation algorithm
- *
- *
- *
- * <ul>
- *   <li style="list-style-type: none;">\f$ \textbf{for } i : dimension \rightarrow 0 \textbf{ do} \f$
- *     <ul>
- *       <li style="list-style-type: none;">\f$\textbf{for all } \sigma of dimension i \f$
- *         <ul>
- *           <li style="list-style-type: none;">\f$\textbf{if } filtration( \sigma ) is NaN \textbf{ then} \f$
- *             <ul>
- *               <li style="list-style-type: none;">\f$ filtration( \sigma ) = \alpha^2( \sigma ) \f$
- *               </li>
- *             </ul>
- *           </li>
- *           <li style="list-style-type: none;">\f$\textbf{end if}\f$
- *           </li>
- *           <li style="list-style-type: none;">\f$\textbf{for all } \tau face of \sigma \textbf{ do} \f$
- *             &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// propagate alpha filtration value
- *             <ul>
- *               <li style="list-style-type: none;">\f$\textbf{if } filtration( \tau ) is not NaN \textbf{ then} \f$
- *                 <ul>
- *                   <li style="list-style-type: none;">\f$ filtration( \tau ) = min ( filtration( \tau ), filtration( \sigma ) ) \f$
- *                   </li>
- *                 </ul>
- *               </li>
- *               <li style="list-style-type: none;">\f$\textbf{else}\f$
- *                 <ul>
- *                   <li style="list-style-type: none;">\f$\textbf{if } \tau is not Gabriel for \sigma \textbf{ then} \f$
- *                     <ul>
- *                       <li style="list-style-type: none;">\f$ filtration( \tau ) = filtration( \sigma ) \f$
- *                       </li>
- *                     </ul>
- *                   </li>
- *                   <li style="list-style-type: none;">\f$\textbf{end if}\f$
- *                   </li>
- *                 </ul>
- *               </li>
- *               <li style="list-style-type: none;">\f$\textbf{end if}\f$
- *               </li>
- *             </ul>
- *           </li>
- *           <li style="list-style-type: none;">\f$\textbf{end for}\f$
- *           </li>
- *         </ul>
- *       </li>
- *       <li style="list-style-type: none;">\f$\textbf{end for}\f$
- *       </li>
- *     </ul>
- *   </li>
- *   <li style="list-style-type: none;">\f$\textbf{end for}\f$
- *   </li>
- *   <li style="list-style-type: none;">\f$make\_filtration\_non\_decreasing()\f$
- *   </li>
- *   <li style="list-style-type: none;">\f$prune\_above\_filtration()\f$
- *   </li>
- * </ul>
- *
+ * <br>
+ * \f$
+ * \textbf{for } \text{i : dimension } \rightarrow 0 \textbf{ do}\\
+ * \quad \textbf{for all } \sigma \text{ of dimension i}\\
+ * \quad\quad \textbf{if } \text{filtration(} \sigma ) \text{ is NaN} \textbf{ then}\\
+ * \quad\quad\quad \text{filtration(} \sigma ) = \alpha^2( \sigma )\\
+ * \quad\quad \textbf{end if}\\
+ * \quad\quad \textbf{for all } \tau \text{ face of } \sigma \textbf{ do}\quad\quad
+ * \textit{// propagate alpha filtration value}\\
+ * \quad\quad\quad \textbf{if } \text{filtration(} \tau ) \text{ is not NaN} \textbf{ then}\\
+ * \quad\quad\quad\quad \text{filtration(} \tau \text{) = min( filtration(} \tau \text{), filtration(} \sigma
+ * \text{) )}\\
+ * \quad\quad\quad \textbf{else}\\
+ * \quad\quad\quad\quad \textbf{if } \textbf{if } \tau \text{ is not Gabriel for } \sigma \textbf{ then}\\
+ * \quad\quad\quad\quad\quad \text{filtration(} \tau \text{) = filtration(} \sigma \text{)}\\
+ * \quad\quad\quad\quad \textbf{end if}\\
+ * \quad\quad\quad \textbf{end if}\\
+ * \quad\quad \textbf{end for}\\
+ * \quad \textbf{end for}\\
+ * \textbf{end for}\\
+ * \text{make_filtration_non_decreasing()}\\
+ * \text{prune_above_filtration()}\\
+ * \f$
  *
  * \subsubsection dimension2 Dimension 2
  * 
@@ -195,8 +165,6 @@ namespace alpha_complex {
  * 
  * \include Alpha_complex/alphaoffreader_for_doc_32.txt
  * 
- * \copyright GNU General Public License v3.                         
- * \verbatim  Contact: gudhi-users@lists.gforge.inria.fr \endverbatim
  */
 /** @} */  // end defgroup alpha_complex
 
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index 5bf553e9..2fc62452 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -1,4 +1,3 @@
-cmake_minimum_required(VERSION 2.6)
 project(Alpha_complex_examples)
 
 # need CGAL 4.7
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 1ff95c3d..4c07eddb 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Vincent Rouvreau
  *
- *    Copyright (C) 2015  INRIA
+ *    Copyright (C) 2015 Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -34,6 +34,7 @@
 #include <CGAL/Epick_d.h>
 #include <CGAL/Spatial_sort_traits_adapter_d.h>
 #include <CGAL/property_map.h>  // for CGAL::Identity_property_map
+#include <CGAL/NT_converter.h>
 
 #include <iostream>
 #include <vector>
@@ -175,7 +176,7 @@ class Alpha_complex {
    *
    * @return The number of vertices.
    */
-  const std::size_t number_of_vertices() const {
+  std::size_t number_of_vertices() const {
     return vertex_handle_to_iterator_.size();
   }
 
@@ -268,8 +269,6 @@ class Alpha_complex {
       return false;  // ----- >>
     }
 
-    complex.set_dimension(triangulation_->maximal_dimension());
-
     // --------------------------------------------------------------------------------------------
     // Simplex_tree construction from loop on triangulation finite full cells list
     if (triangulation_->number_of_vertices() > 0) {
@@ -325,8 +324,9 @@ 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 = squared_radius(pointVector.begin(), pointVector.end());
+              alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end()));
             }
             complex.assign_filtration(f_simplex, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
index 7380547f..622fcae8 100644
--- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Vincent Rouvreau
  *
- *    Copyright (C) 2015  INRIA
+ *    Copyright (C) 2015 Inria
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -23,6 +23,7 @@
 #define BOOST_TEST_DYN_LINK
 #define BOOST_TEST_MODULE "alpha_complex"
 #include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
 
 #include <CGAL/Delaunay_triangulation.h>
 #include <CGAL/Epick_d.h>
@@ -36,7 +37,7 @@
 // to construct a simplex_tree from Delaunay_triangulation
 #include <gudhi/graph_simplicial_complex.h>
 #include <gudhi/Simplex_tree.h>
-#include <boost/mpl/list.hpp>
+#include <gudhi/Unitary_tests_utils.h>
 
 // Use dynamic_dimension_tag for the user to be able to set dimension
 typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d;
@@ -96,10 +97,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of
   BOOST_CHECK(simplex_tree_59.num_simplices() == 23);
 }
 
-bool are_almost_the_same(float a, float b) {
-  return std::fabs(a - b) < std::numeric_limits<float>::epsilon();
-}
-
 // Use static dimension_tag for the user not to be able to set dimension
 typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Kernel_4;
 typedef Kernel_4::Point_d Point_4;
@@ -159,23 +156,23 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   BOOST_CHECK(simplex_tree.num_simplices() == 15);
 
   std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
-  BOOST_CHECK(simplex_tree.dimension() == 4);
+  BOOST_CHECK(simplex_tree.dimension() == 3);
   std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() == 4);
 
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
     switch (simplex_tree.dimension(f_simplex)) {
       case 0:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 0.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 0.0);
         break;
       case 1:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 1.0/2.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 1.0/2.0);
         break;
       case 2:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 2.0/3.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 2.0/3.0);
         break;
       case 3:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 3.0/4.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 3.0/4.0);
         break;
       default:
         BOOST_CHECK(false);  // Shall not happen
@@ -232,17 +229,17 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   BOOST_CHECK(simplex_tree.num_simplices() == 10);
 
   std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
-  BOOST_CHECK(simplex_tree.dimension() == 4);
+  BOOST_CHECK(simplex_tree.dimension() == 1);
   std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() == 4);
 
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
     switch (simplex_tree.dimension(f_simplex)) {
       case 0:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 0.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 0.0);
         break;
       case 1:
-        BOOST_CHECK(are_almost_the_same(simplex_tree.filtration(f_simplex), 1.0/2.0));
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 1.0/2.0);
         break;
       default:
         BOOST_CHECK(false);  // Shall not happen
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index 9e0b3b3c..9255d3db 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -1,4 +1,3 @@
-cmake_minimum_required(VERSION 2.6)
 project(Alpha_complex_tests)
 
 if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
new file mode 100644
index 00000000..7ace6064
--- /dev/null
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -0,0 +1,64 @@
+project(Alpha_complex_utilities)
+
+if(CGAL_FOUND)
+  add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp)
+  target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp)
+  target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  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})
+
+  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})
+  endif(TBB_FOUND)
+  
+  add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
+  add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:exact_alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"  "-p" "2" "-m" "0.45")
+  add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights"  "-p" "2" "-m" "0.45")
+
+  install(TARGETS alpha_complex_3d_persistence DESTINATION bin)
+  install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin)
+  install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin)
+
+  if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
+    add_executable (alpha_complex_persistence alpha_complex_persistence.cpp)
+    target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+    add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp)
+    target_link_libraries(periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
+    if (TBB_FOUND)
+      target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES})
+      target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    endif(TBB_FOUND)
+    add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
+    add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:periodic_alpha_complex_3d_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "-p" "2" "-m" "0")
+
+    install(TARGETS alpha_complex_persistence DESTINATION bin)
+    install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin)
+
+  endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
+
+  if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+    add_executable(weighted_periodic_alpha_complex_3d_persistence weighted_periodic_alpha_complex_3d_persistence.cpp)
+    target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY})
+    if (TBB_FOUND)
+      target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES})
+    endif(TBB_FOUND)
+
+    add_test(NAME Alpha_complex_utilities_weigted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_periodic_alpha_complex_3d_persistence>
+        "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights"
+        "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "3" "1.0")
+
+    install(TARGETS weighted_periodic_alpha_complex_3d_persistence DESTINATION bin)
+
+  endif (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+
+endif(CGAL_FOUND)
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h
new file mode 100644
index 00000000..a72fd96d
--- /dev/null
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h
@@ -0,0 +1,74 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ALPHA_COMPLEX_3D_HELPER_H_
+#define ALPHA_COMPLEX_3D_HELPER_H_
+
+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;
+}
+
+#endif  // ALPHA_COMPLEX_3D_HELPER_H_
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..6e603155
--- /dev/null
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -0,0 +1,271 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    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/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 <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
+
+// 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 Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, 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);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  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);
+  // 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;
+
+  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
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+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) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    exit(-1);
+  }
+}
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
new file mode 100644
index 00000000..8e6c40b7
--- /dev/null
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -0,0 +1,138 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2016 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <boost/program_options.hpp>
+
+#include <CGAL/Epick_d.h>
+
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <iostream>
+#include <string>
+#include <limits>  // for numeric_limits
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Filtration_value = Simplex_tree::Filtration_value;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  Filtration_value alpha_square_max_value;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, coeff_field_characteristic,
+                  min_persistence);
+
+  // ----------------------------------------------------------------------------
+  // 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_points);
+
+  Simplex_tree simplex;
+  if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
+    // ----------------------------------------------------------------------------
+    // Display information about the alpha complex
+    // ----------------------------------------------------------------------------
+    std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+              << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
+
+    // Sort the simplices in the order of the filtration
+    simplex.initialize_filtration();
+
+    std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+    // Compute the persistence diagram of the complex
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
+        simplex);
+    // initializes the coefficient field for homology
+    pcoh.init_coefficients(coeff_field_characteristic);
+
+    pcoh.compute_persistent_cohomology(min_persistence);
+
+    // Output the diagram in filediag
+    if (output_file_diag.empty()) {
+      pcoh.output_diagram();
+    } else {
+      std::cout << "Result in file: " << output_file_diag << std::endl;
+      std::ofstream out(output_file_diag);
+      pcoh.output_diagram(out);
+      out.close();
+    }
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag,
+                     Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value)
+                                      ->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal alpha square value for the Alpha complex construction.")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of an Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    exit(-1);
+  }
+}
diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md
new file mode 100644
index 00000000..b1a33e4b
--- /dev/null
+++ b/src/Alpha_complex/utilities/alphacomplex.md
@@ -0,0 +1,165 @@
+---

+layout: page

+title: "Alpha complex"

+meta_title: "Alpha complex"

+teaser: ""

+permalink: /alphacomplex/

+---

+{::comment}

+Leave the lines above as it is required by the web site generator 'Jekyll'

+{:/comment}

+

+

+## alpha_complex_persistence ##

+

+This program computes the persistent homology with coefficient field Z/pZ of the dD alpha complex built from a dD point cloud.

+The output diagram contains one bar per line, written with the convention:

+

+```

+   p dim birth death

+```

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature,

+and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).

+

+**Usage**

+

+```

+   alpha_complex_persistence [options] <input OFF file>

+```

+

+where

+`<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-r [ --max-alpha-square-value ]` (default = inf) Maximal alpha square value for the Alpha complex construction.

+* `-p [ --field-charac ]` (default = 11)     Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+   alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off

+```

+

+N.B.:

+

+* Filtration values are alpha square values.

+

+

+## alpha_complex_3d_persistence ##

+This program computes the persistent homology with coefficient field Z/pZ of the 3D alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention:

+

+```

+p dim birth death

+```

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).

+

+**Usage**

+

+```

+   alpha_complex_3d_persistence [options] <input OFF file>

+```

+

+where `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45

+```

+

+N.B.:

+

+* `alpha_complex_3d_persistence` only accepts OFF files in dimension 3.

+* Filtration values are alpha square values.

+

+

+## exact_alpha_complex_3d_persistence ##

+

+Same as `alpha_complex_3d_persistence`, but using exact computation.

+It is slower, but it is necessary when points are on a grid for instance.

+

+

+

+## weighted_alpha_complex_3d_persistence ##

+

+Same as `alpha_complex_3d_persistence`, but using weighted points.

+

+**Usage**

+

+```

+   weighted_alpha_complex_3d_persistence [options] <input OFF file> <weights input file>

+```

+

+where

+

+* `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+* `<input weights file>` is the path to the file containing the weights of the points (one value per line).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

+

+**Example**

+

+```

+   weighted_alpha_complex_3d_persistence ../../data/points/tore3D_300.off ../../data/points/tore3D_300.weights -p 2 -m 0.45

+```

+

+

+N.B.:

+

+* Weights values are explained on CGAL [Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0)

+and [Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation.

+* Filtration values are alpha square values.

+

+

+## periodic_alpha_complex_3d_persistence ##

+Same as `alpha_complex_3d_persistence`, but using periodic alpha shape 3d.

+Refer to the [CGAL's 3D Periodic Triangulations User Manual](https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html) for more details.

+

+**Usage**

+

+```

+   periodic_alpha_complex_3d_persistence [options] <input OFF file> <cuboid file>

+```

+

+where

+

+* `<input OFF file>` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).

+* `<cuboid file>` is the path to the file describing the periodic domain. It must be in the format described

+[here]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsIsoCuboid).

+

+**Allowed options**

+

+* `-h [ --help ]` Produce help message

+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output.

+* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals

+

+

+**Example**

+

+```

+periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt -p 3 -m 1.0

+```

+

+N.B.:

+

+* Cuboid file must be in the format described [here]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsIsoCuboid).

+* Filtration values are alpha square values.

diff --git a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..61f49bb1
--- /dev/null
+++ b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,265 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#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 <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 Exact_tag = CGAL::Tag_true;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, CGAL::Default, Exact_tag>;
+using Fb = CGAL::Alpha_shape_cell_base_3<Kernel, CGAL::Default, Exact_tag>;
+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, Exact_tag>;
+
+// 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;
+using Vertex_list = std::vector<Vertex_handle>;
+
+// 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 Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, 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);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  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);
+  // 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
+    // you can also use the_alpha_value_iterator->exact()
+    Filtration_value filtr = /*std::sqrt*/ CGAL::to_double(the_alpha_value_iterator->exact());
+#ifdef DEBUG_TRACES
+    std::cout << "filtration = " << filtr << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree.insert_simplex(the_simplex, filtr);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#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;
+
+  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
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+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) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    exit(-1);
+  }
+}
diff --git a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..a261c5a3
--- /dev/null
+++ b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,302 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *                     Pawel Dlotko - 2017 - Swansea University, UK
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#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/Periodic_3_Delaunay_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_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 <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// 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<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<PK, DsCb>;
+using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>;
+using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>;
+using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3<PK, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<P3DT3>;
+using Point_3 = PK::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;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// 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 Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string cuboid_file;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, cuboid_file, 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);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Read iso_cuboid_3 information from file
+  std::ifstream iso_cuboid_str(cuboid_file);
+  double x_min, y_min, z_min, x_max, y_max, z_max;
+  if (iso_cuboid_str.good()) {
+    iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max;
+  } else {
+    std::cerr << "Unable to read file " << cuboid_file << std::endl;
+    exit(-1);
+  }
+  // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
+  if ((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::cerr << "The size of the cuboid in every directions is not the same." << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Define the periodic cube
+  P3DT3 pdt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
+  // Heuristic for inserting large point sets (if pts is reasonably large)
+  pdt.insert(lp.begin(), lp.end(), true);
+  // As pdt won't be modified anymore switch to 1-sheeted cover if possible
+  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;
+
+  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);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#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;
+
+  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
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ")(
+      "cuboid-file", po::value<std::string>(&cuboid_file),
+      "Name of file describing the periodic domain. Format is: min_hx min_hy min_hz\nmax_hx max_hy max_hz");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+  pos.add("cuboid-file", 2);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file") || !vm.count("cuboid-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a periodic 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file cuboid-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    exit(-1);
+  }
+}
diff --git a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..aa7ddee2
--- /dev/null
+++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,316 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <boost/program_options.hpp>
+#include <boost/variant.hpp>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_3D_off_io.h>
+
+#include <CGAL/config.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Regular_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>
+
+// For CGAL < 4.11
+#if CGAL_VERSION_NR < 1041100000
+#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#endif  // CGAL_VERSION_NR < 1041100000
+
+#include <fstream>
+#include <cmath>
+#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;
+
+// For CGAL < 4.11
+#if CGAL_VERSION_NR < 1041100000
+using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>;
+using Fb = CGAL::Alpha_shape_cell_base_3<Gt>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
+using Triangulation_3 = CGAL::Regular_triangulation_3<Gt, Tds>;
+
+// From file type definition
+using Point_3 = Gt::Bare_point;
+using Weighted_point_3 = Gt::Weighted_point;
+
+// For CGAL >= 4.11
+#else   // CGAL_VERSION_NR < 1041100000
+using Rvb = CGAL::Regular_triangulation_vertex_base_3<Kernel>;
+using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Rvb>;
+using Rcb = CGAL::Regular_triangulation_cell_base_3<Kernel>;
+using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Rcb>;
+using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;
+using Triangulation_3 = CGAL::Regular_triangulation_3<Kernel, Tds>;
+
+// From file type definition
+using Point_3 = Triangulation_3::Bare_point;
+using Weighted_point_3 = Triangulation_3::Weighted_point;
+#endif  // CGAL_VERSION_NR < 1041100000
+
+using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_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;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// 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 Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence);
+
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string weight_file;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, weight_file, 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);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
+    exit(-1);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Read weights information from file
+  std::ifstream weights_ifstr(weight_file);
+  std::vector<Weighted_point_3> wp;
+  if (weights_ifstr.good()) {
+    double weight = 0.0;
+    std::size_t index = 0;
+    wp.reserve(lp.size());
+    // Attempt read the weight in a double format, return false if it fails
+    while ((weights_ifstr >> weight) && (index < lp.size())) {
+      wp.push_back(Weighted_point_3(lp[index], weight));
+      index++;
+    }
+    if (index != lp.size()) {
+      std::cerr << "Bad number of weights in file " << weight_file << std::endl;
+      exit(-1);
+    }
+  } else {
+    std::cerr << "Unable to read weights file " << weight_file << std::endl;
+    exit(-1);
+  }
+
+  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
+  Alpha_shape_3 as(wp.begin(), wp.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);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#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;
+
+  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
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (output_file_diag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::ofstream out(output_file_diag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file,
+                     std::string &output_file_diag, int &coeff_field_characteristic,
+                     Filtration_value &min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ")(
+      "weight-file", po::value<std::string>(&weight_file),
+      "Name of file containing a point weights. Format is one weigt per line:   W1\n...\nWn ");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+  pos.add("weight-file", 2);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file") || !vm.count("weight-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a weighted 3D Alpha complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file weight-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    exit(-1);
+  }
+}
diff --git a/src/Alpha_complex/utilities/weighted_periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_periodic_alpha_complex_3d_persistence.cpp
new file mode 100644
index 00000000..d030c88c
--- /dev/null
+++ b/src/Alpha_complex/utilities/weighted_periodic_alpha_complex_3d_persistence.cpp
@@ -0,0 +1,288 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *                     Pawel Dlotko - 2017 - Swansea University, UK
+ *
+ *    Copyright (C) 2014 Inria
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <boost/variant.hpp>
+
+#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/Periodic_3_regular_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_regular_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 <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <vector>
+#include <cstdlib>
+
+#include "alpha_complex_3d_helper.h"
+
+// Traits
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using PK = CGAL::Periodic_3_regular_triangulation_traits_3<Kernel>;
+
+// Vertex type
+using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
+using Vb = CGAL::Regular_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::Regular_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>;
+using P3RT3 = CGAL::Periodic_3_regular_triangulation_3<PK, Tds>;
+using Alpha_shape_3 = CGAL::Alpha_shape_3<P3RT3>;
+
+using Point_3 = P3RT3::Bare_point;
+using Weighted_point_3 = P3RT3::Weighted_point;
+
+// 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;
+using Vertex_list = std::vector<Alpha_shape_3::Vertex_handle>;
+
+// 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 Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
+
+void usage(const std::string& progName) {
+  std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_weight_file path_to_the_cuboid_file "
+                                        "coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n";
+  exit(-1);
+}
+
+int main(int argc, char* const argv[]) {
+  // program args management
+  if (argc != 6) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
+    usage(argv[0]);
+  }
+
+  int coeff_field_characteristic = atoi(argv[4]);
+  Filtration_value min_persistence = strtof(argv[5], nullptr);
+
+  // Read points from file
+  std::string offInputFile(argv[1]);
+  // Read the OFF file (input file name given as parameter) and triangulate points
+  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
+  // Check the read operation was correct
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << offInputFile << std::endl;
+    usage(argv[0]);
+  }
+
+  // Retrieve the points
+  std::vector<Point_3> lp = off_reader.get_point_cloud();
+
+  // Read iso_cuboid_3 information from file
+  std::ifstream iso_cuboid_str(argv[3]);
+  double x_min, y_min, z_min, x_max, y_max, z_max;
+  if (iso_cuboid_str.is_open()) {
+    if (!(iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max)) {
+      std::cerr << argv[3] << " - Bad file format." << std::endl;
+      usage(argv[0]);
+    }
+
+  } else {
+    std::cerr << "Unable to read file " << argv[3] << std::endl;
+    usage(argv[0]);
+  }
+  // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
+  if ((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::cerr << "The size of the cuboid in every directions is not the same." << std::endl;
+    exit(-1);
+  }
+
+  double maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
+
+  // Read weights information from file
+  std::ifstream weights_ifstr(argv[2]);
+  std::vector<Weighted_point_3> wp;
+  if (weights_ifstr.is_open()) {
+    double weight = 0.0;
+    std::size_t index = 0;
+    wp.reserve(lp.size());
+    // Attempt read the weight in a double format, return false if it fails
+    while ((weights_ifstr >> weight) && (index < lp.size())) {
+      if ((weight >= maximal_possible_weight) || (weight < 0)) {
+        std::cerr << "At line " << (index + 1) << ", the weight (" << weight
+                  << ") is negative or more than or equal to maximal possible weight (" << maximal_possible_weight
+                  << ") = 1/64*cuboid length squared, which is not an acceptable input." << std::endl;
+        exit(-1);
+      }
+
+      wp.push_back(Weighted_point_3(lp[index], weight));
+      index++;
+    }
+    if (index != lp.size()) {
+      std::cerr << "Bad number of weights in file " << argv[2] << std::endl;
+      usage(argv[0]);
+    }
+  } else {
+    std::cerr << "Unable to read file " << argv[2] << std::endl;
+    usage(argv[0]);
+  }
+
+  // Define the periodic cube
+  P3RT3 prt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
+  // Heuristic for inserting large point sets (if pts is reasonably large)
+  prt.insert(wp.begin(), wp.end(), true);
+  // As prt won't be modified anymore switch to 1-sheeted cover if possible
+  if (prt.is_triangulation_in_1_sheet()) {
+    prt.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 << "Weighted 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(prt, 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 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);
+    if (the_alpha_value_iterator != the_alpha_values.end())
+      ++the_alpha_value_iterator;
+    else
+      std::cout << "This shall not happen" << std::endl;
+  }
+
+#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;
+
+  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
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree, true);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  pcoh.output_diagram();
+
+  return 0;
+}