Code review : Homogenize alpha shape examples command line

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

Former-commit-id: 8c692434eb48c525047bc34a6fab7d7c1af64d56

10 files changed, 406 insertions, 252 deletions

diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index 18699650..cefc1ad6 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -3,11 +3,11 @@ 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})
+  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})
+  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})
+  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})
@@ -16,11 +16,11 @@ if(CGAL_FOUND)
   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" "2" "0.45")
+      "${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" "2" "0.45")
+      "${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" "2" "0.45")
+      "${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)
@@ -28,11 +28,10 @@ if(CGAL_FOUND)
 
   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})
+    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})
+    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})
@@ -41,7 +40,7 @@ if(CGAL_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" "2" "0")
+        "${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)
diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README
index 6e7d16b7..c3dd170b 100644
--- a/src/Alpha_complex/utilities/README
+++ b/src/Alpha_complex/utilities/README
@@ -8,22 +8,24 @@ This program computes the persistent homology with coefficient field Z/pZ of the
 where `dim` is the dimension of the homological feature, `b` and `d` 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 <input OFF file> <p> <min_persistence>` 
-where
-`<input OFF file>` is the path to the input point cloud in OFF format.
-`<p>` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer.
-`<min_persistence>` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
+`alpha_complex_3d_persistence [options] <OFF input file>`
+
+**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 2 0.45`
+`alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45`
 
 outputs:
 ```
 Simplex_tree dim: 3
-2  0 0 inf 
-2  1 0.0682162 1.0001 
-2  1 0.0934117 1.00003 
-2  2 0.56444 1.03938 
+2  0 0 inf
+2  1 0.0682162 1.0001
+2  1 0.0934117 1.00003
+2  2 0.56444 1.03938
 ```
 
 Here we retrieve expected Betti numbers on a tore 3D:
@@ -48,21 +50,34 @@ Same as `alpha_complex_3d_persistence`, but using exact computation. It is slowe
 Same as `alpha_complex_3d_persistence`, but using weighted points.
 
 **Usage**
-`weighted_alpha_complex_3d_persistence <input OFF file> <input weights file> <p> <min_persistence>` 
-where
-`<input OFF file>` is the path to the input point cloud in OFF format.
-`<input weights file>` is the path to the file containing the weights of the points (one value per line).
-`<p>` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer.
-`<min_persistence>` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
+`weighted_alpha_complex_3d_persistence [options] <OFF input file> <weights input file>`
 
+**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`
 
-## `periodic_alpha_complex_3d_persistence` ##
-This program computes the persistent homology with coefficient field Z/pZ of the 3D periodic alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention:
+outputs:
+```
+Simplex_tree dim: 3
+2  0 -1 inf
+2  1 -0.931784 0.000103311
+2  1 -0.906588 2.60165e-05
+2  2 -0.43556 0.0393798
+```
 
-`p dim b d`
+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.
 
-where `dim` is the dimension of the homological feature, `b` and `d` 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).
+
+## `periodic_alpha_complex_3d_persistence` ##
+Same as `alpha_complex_3d_persistence`, but using periodic alpha shape 3d.
 
 **Usage**
 `periodic_alpha_complex_3d_persistence <input OFF file> <cuboid file> <p> <min_persistence>`
@@ -72,21 +87,31 @@ where
 `<p>` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer.
 `<min_persistence>` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
 
+**Usage**
+`./periodic_alpha_complex_3d_persistence [options] input-file cuboid-file`
+
+**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 3 1.0`
+`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`
 
 outputs:
 ```
 Periodic Delaunay computed.
 Simplex_tree dim: 3
-3  0 0 inf 
-3  1 0.0025 inf 
-3  1 0.0025 inf 
-3  1 0.0025 inf 
-3  2 0.005 inf 
-3  2 0.005 inf 
-3  2 0.005 inf 
-3  3 0.0075 inf 
+3  0 0 inf
+3  1 0.0025 inf
+3  1 0.0025 inf
+3  1 0.0025 inf
+3  2 0.005 inf
+3  2 0.005 inf
+3  2 0.005 inf
+3  3 0.0075 inf
 ```
 
 Here we retrieve expected Betti numbers on an 3D iso-oriented cuboids:
@@ -98,14 +123,12 @@ Betti numbers[3] = 1
 ```
 
 N.B.: 
-* `periodic_alpha_complex_3d_persistence` only accepts OFF files in dimension 3.
-* In this example, the periodic cube is hard coded to { x = [0,1]; y = [0,1]; z = [0,1] }
+* Cuboid file must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsIsoCuboid).
 * Filtration values are alpha square values.
 
 
 
 
-
 ## `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:
 
@@ -119,7 +142,7 @@ where `dim` is the dimension of the homological feature, `b` and `d` are respect
 **Allowed options**
 
 * `-h [ --help ]` Produce help message
-* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout.
+* `-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.
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h
index 7865e4ec..6b3b7d5d 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h
@@ -23,7 +23,7 @@
 #ifndef ALPHA_COMPLEX_3D_HELPER_H_
 #define ALPHA_COMPLEX_3D_HELPER_H_
 
-template<class Vertex_list, class Cell_handle>
+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++) {
@@ -35,7 +35,7 @@ Vertex_list from_cell(const Cell_handle& ch) {
   return the_list;
 }
 
-template<class Vertex_list, class Facet>
+template <class Vertex_list, class Facet>
 Vertex_list from_facet(const Facet& fct) {
   Vertex_list the_list;
   for (auto i = 0; i < 4; i++) {
@@ -49,7 +49,7 @@ Vertex_list from_facet(const Facet& fct) {
   return the_list;
 }
 
-template<class Vertex_list, class Edge_3>
+template <class Vertex_list, class Edge_3>
 Vertex_list from_edge(const Edge_3& edg) {
   Vertex_list the_list;
   for (auto i = 0; i < 4; i++) {
@@ -63,7 +63,7 @@ Vertex_list from_edge(const Edge_3& edg) {
   return the_list;
 }
 
-template<class Vertex_list, class Vertex_handle>
+template <class Vertex_list, class Vertex_handle>
 Vertex_list from_vertex(const Vertex_handle& vh) {
   Vertex_list the_list;
 #ifdef DEBUG_TRACES
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
index f63ff0f6..f2085ab2 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -20,6 +20,7 @@
  *    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>
@@ -56,10 +57,10 @@ 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 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;
@@ -70,42 +71,28 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
 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 Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
 using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(const std::string& progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  exit(-1);
-}
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
 
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 4) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[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);
 
-  int coeff_field_characteristic = atoi(argv[2]);
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
 
-  Filtration_value min_persistence = 0.0;
-  int returnedScanValue = sscanf(argv[3], "%f", &min_persistence);
-  if ((returnedScanValue == EOF) || (min_persistence < -1.0)) {
-    std::cerr << "Error: " << argv[3] << " is not correct\n";
-    usage(argv[0]);
-  }
+  program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence);
 
-  // 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);
+  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 " << offInputFile << std::endl;
-    usage(argv[0]);
+    std::cerr << "Unable to read file " << off_file_points << std::endl;
+    exit(-1);
   }
 
   // Retrieve the triangulation
@@ -143,28 +130,28 @@ int main(int argc, char * const argv[]) {
   Filtration_value filtration_max = 0.0;
   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)) {
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
       vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
       count_cells++;
       if (dim_max < 3) {
         // Cell is of dim 3
         dim_max = 3;
       }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
       vertex_list = from_facet<Vertex_list, Facet>(*facet);
       count_facets++;
       if (dim_max < 2) {
         // Facet is of dim 2
         dim_max = 2;
       }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
       vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
       count_edges++;
       if (dim_max < 1) {
         // Edge_3 is of dim 1
         dim_max = 1;
       }
-    } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
       count_vertices++;
       vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
     }
@@ -190,7 +177,7 @@ int main(int argc, char * const argv[]) {
       }
     }
     // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
 #ifdef DEBUG_TRACES
     std::cout << "filtration = " << filtr << std::endl;
 #endif  // DEBUG_TRACES
@@ -211,7 +198,6 @@ int main(int argc, char * const argv[]) {
   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;
@@ -236,7 +222,58 @@ int main(int argc, char * const argv[]) {
 
   pcoh.compute_persistent_cohomology(min_persistence);
 
-  pcoh.output_diagram();
+  // 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;
+    std::abort();
+  }
+}
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index 9e84e91f..2105220a 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -14,12 +14,9 @@
 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);
+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;
@@ -28,13 +25,13 @@ int main(int argc, char **argv) {
   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);
+  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 >;
+  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;
@@ -42,17 +39,16 @@ int main(int argc, char **argv) {
     // ----------------------------------------------------------------------------
     // 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;
+    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);
+    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);
 
@@ -72,30 +68,26 @@ int main(int argc, char **argv) {
   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) {
+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 ");
+  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");
+  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);
@@ -104,8 +96,7 @@ void program_options(int argc, char * argv[]
   all.add(visible).add(hidden);
 
   po::variables_map vm;
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), 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")) {
diff --git a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp
index 8ce68406..7919c7fd 100644
--- a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp
@@ -20,6 +20,7 @@
  *    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>
@@ -40,7 +41,7 @@
 #include <list>
 #include <vector>
 
-#include "../utilities/alpha_complex_3d_helper.h"
+#include "alpha_complex_3d_helper.h"
 
 // Alpha_shape_3 templates type definitions
 using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
@@ -57,10 +58,10 @@ 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 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;
@@ -71,42 +72,28 @@ using Vertex_list = std::list<Vertex_handle>;
 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 Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
 using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  exit(-1);
-}
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
 
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 4) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[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);
 
-  int coeff_field_characteristic = atoi(argv[2]);
+int main(int argc, char **argv) {
+  std::string off_file_points;
+  std::string output_file_diag;
+  int coeff_field_characteristic;
+  Filtration_value min_persistence;
 
-  Filtration_value min_persistence = 0.0;
-  int returnedScanValue = sscanf(argv[3], "%f", &min_persistence);
-  if ((returnedScanValue == EOF) || (min_persistence < -1.0)) {
-    std::cerr << "Error: " << argv[3] << " is not correct\n";
-    usage(argv[0]);
-  }
+  program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence);
 
-  // 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);
+  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 " << offInputFile << std::endl;
-    usage(argv[0]);
+    std::cerr << "Unable to read file " << off_file_points << std::endl;
+    exit(-1);
   }
 
   // Retrieve the triangulation
@@ -144,28 +131,28 @@ int main(int argc, char * const argv[]) {
   Filtration_value filtration_max = 0.0;
   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)) {
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
       vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
       count_cells++;
       if (dim_max < 3) {
         // Cell is of dim 3
         dim_max = 3;
       }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
       vertex_list = from_facet<Vertex_list, Facet>(*facet);
       count_facets++;
       if (dim_max < 2) {
         // Facet is of dim 2
         dim_max = 2;
       }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
       vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
       count_edges++;
       if (dim_max < 1) {
         // Edge_3 is of dim 1
         dim_max = 1;
       }
-    } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
+    } else if (const Vertex_handle *vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
       count_vertices++;
       vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
     }
@@ -192,7 +179,7 @@ int main(int argc, char * const argv[]) {
     }
     // 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());
+    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
@@ -213,7 +200,6 @@ int main(int argc, char * const argv[]) {
   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;
@@ -242,3 +228,46 @@ int main(int argc, char * const argv[]) {
 
   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;
+    std::abort();
+  }
+}
diff --git a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp
index 8140a3c5..10b62f75 100644
--- a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp
@@ -20,6 +20,7 @@
  *    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>
@@ -63,10 +64,10 @@ 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 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;
@@ -77,54 +78,41 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
 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 Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
 using Alpha_shape_simplex_tree_pair = std::pair<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(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_iso_cuboid_3_file coeff_field_characteristic[" <<
-               "integer > 0] min_persistence[float >= -1.0]\n" <<
-               "  path_to_OFF_file is the path to your points cloud in OFF format.\n" <<
-               "  path_to_iso_cuboid_3_file is the path to the iso cuboid file with the following format :\n" <<
-               "    x_min y_min z_min x_max y_max z_max\n" <<
-               "  In this example, the periodic cube will be " <<
-               "{ x = [x_min,x_max]; y = [y_min,y_max]; z = [z_min,z_max] }.\n" <<
-               "  For more information, please refer to\n" <<
-               "    https://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Iso__cuboid__3.html\n";
-
-  exit(-1);
-}
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
 
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 5) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[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);
 
-  int coeff_field_characteristic = atoi(argv[3]);
-  Filtration_value min_persistence = strtof(argv[4], nullptr);
+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 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);
+  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 " << offInputFile << std::endl;
-    usage(argv[0]);
+    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(argv[2]);
+  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 " << argv[2] << std::endl;
-    usage(argv[0]);
+    std::cerr << "Unable to read file " << cuboid_file << std::endl;
+    exit(-1);
   }
 
   // Retrieve the triangulation
@@ -168,29 +156,29 @@ int main(int argc, char * const argv[]) {
   Filtration_value filtration_max = 0.0;
   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)) {
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
       vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
       count_cells++;
       if (dim_max < 3) {
         // Cell is of dim 3
         dim_max = 3;
       }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
       vertex_list = from_facet<Vertex_list, Facet>(*facet);
       count_facets++;
       if (dim_max < 2) {
         // Facet is of dim 2
         dim_max = 2;
       }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
       vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
       count_edges++;
       if (dim_max < 1) {
         // Edge_3 is of dim 1
         dim_max = 1;
       }
-    } else if (const Alpha_shape_3::Vertex_handle * vertex =
-               CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
+    } else if (const Alpha_shape_3::Vertex_handle *vertex =
+                   CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
       count_vertices++;
       vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
     }
@@ -216,7 +204,7 @@ int main(int argc, char * const argv[]) {
       }
     }
     // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
 #ifdef DEBUG_TRACES
     std::cout << "filtration = " << filtr << std::endl;
 #endif  // DEBUG_TRACES
@@ -237,7 +225,6 @@ int main(int argc, char * const argv[]) {
   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;
@@ -266,3 +253,50 @@ int main(int argc, char * const argv[]) {
 
   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;
+    std::abort();
+  }
+}
diff --git a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp
index a4ecf9da..293170f7 100644
--- a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp
@@ -20,6 +20,7 @@
  *    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>
@@ -42,7 +43,7 @@
 #include <vector>
 #include <cstdlib>
 
-#include "../utilities/alpha_complex_3d_helper.h"
+#include "alpha_complex_3d_helper.h"
 
 // Traits
 using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
@@ -60,10 +61,10 @@ using Weighted_point_3 = Gt::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 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;
@@ -74,48 +75,38 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
 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 Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
 using Alpha_shape_simplex_tree_pair = std::pair<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(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  std::cerr << "  path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n";
-  std::cerr << "  Weights values are explained on CGAL documentation:\n";
-  std::cerr << "    https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n";
-  std::cerr << "    https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\n";
-  exit(-1);
-}
+using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>;
+using Persistent_cohomology =
+    Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>;
 
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 5) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[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);
 
-  int coeff_field_characteristic = atoi(argv[3]);
-  Filtration_value min_persistence = strtof(argv[4], nullptr);
+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 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);
+  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 " << offInputFile << std::endl;
-    usage(argv[0]);
+    std::cerr << "Unable to read OFF file " << off_file_points << std::endl;
+    exit(-1);
   }
 
   // Retrieve the triangulation
   std::vector<Point_3> lp = off_reader.get_point_cloud();
 
   // Read weights information from file
-  std::ifstream weights_ifstr(argv[2]);
+  std::ifstream weights_ifstr(weight_file);
   std::vector<Weighted_point_3> wp;
   if (weights_ifstr.good()) {
     double weight = 0.0;
@@ -127,12 +118,12 @@ int main(int argc, char * const argv[]) {
       index++;
     }
     if (index != lp.size()) {
-      std::cerr << "Bad number of weights in file " << argv[2] << std::endl;
-      usage(argv[0]);
+      std::cerr << "Bad number of weights in file " << weight_file << std::endl;
+      exit(-1);
     }
   } else {
-    std::cerr << "Unable to read file " << argv[2] << std::endl;
-    usage(argv[0]);
+    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.
@@ -167,29 +158,29 @@ int main(int argc, char * const argv[]) {
   Filtration_value filtration_max = 0.0;
   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)) {
+    if (const Cell_handle *cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
       vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
       count_cells++;
       if (dim_max < 3) {
         // Cell is of dim 3
         dim_max = 3;
       }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
+    } else if (const Facet *facet = CGAL::object_cast<Facet>(&object_iterator)) {
       vertex_list = from_facet<Vertex_list, Facet>(*facet);
       count_facets++;
       if (dim_max < 2) {
         // Facet is of dim 2
         dim_max = 2;
       }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
+    } else if (const Edge_3 *edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
       vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
       count_edges++;
       if (dim_max < 1) {
         // Edge_3 is of dim 1
         dim_max = 1;
       }
-    } else if (const Alpha_shape_3::Vertex_handle * vertex =
-               CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
+    } else if (const Alpha_shape_3::Vertex_handle *vertex =
+                   CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
       count_vertices++;
       vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
     }
@@ -215,7 +206,7 @@ int main(int argc, char * const argv[]) {
       }
     }
     // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
+    Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator);
 #ifdef DEBUG_TRACES
     std::cout << "filtration = " << filtr << std::endl;
 #endif  // DEBUG_TRACES
@@ -236,7 +227,6 @@ int main(int argc, char * const argv[]) {
   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;
@@ -265,3 +255,50 @@ int main(int argc, char * const argv[]) {
 
   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;
+    std::abort();
+  }
+}
diff --git a/src/Bottleneck_distance/utilities/CMakeLists.txt b/src/Bottleneck_distance/utilities/CMakeLists.txt
index cdf76c85..063b6ae3 100644
--- a/src/Bottleneck_distance/utilities/CMakeLists.txt
+++ b/src/Bottleneck_distance/utilities/CMakeLists.txt
@@ -3,6 +3,9 @@ project(Bottleneck_distance_utilities)
 
 if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
   add_executable (bottleneck_read_file bottleneck_read_file.cpp)
+  if (TBB_FOUND)
+    target_link_libraries(bottleneck_read_file ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
 
   add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file
       COMMAND $<TARGET_FILE:bottleneck_read_file>
diff --git a/src/Rips_complex/utilities/README b/src/Rips_complex/utilities/README
index 93ec3658..ca10424d 100644
--- a/src/Rips_complex/utilities/README
+++ b/src/Rips_complex/utilities/README
@@ -13,6 +13,7 @@ where `dim` is the dimension of the homological feature, `b` and `d` are respect
 **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-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.
 * `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.
 * `-p [ --field-charac ]` (default = 11)     Characteristic p of the coefficient field Z/pZ for computing homology.
@@ -45,7 +46,7 @@ outputs:
 
 
 ## `rips_distance_matrix_persistence` ##
-Same as `rips_persistence` but taking an distance matrix as input.
+Same as `rips_persistence` but taking a distance matrix as input.
 
 **Example**
 `rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0`