Merge pull request #355 from VincentRouvreau/alpha_complex_3d_python

Alpha complex is using 3d version in dim 3

9 files changed, 322 insertions, 175 deletions

diff --git a/.appveyor.yml b/.appveyor.yml
index d072a366..a257debc 100644
--- a/.appveyor.yml
+++ b/.appveyor.yml
@@ -11,23 +11,23 @@ configuration:
 environment:
   # update the vcpkg cache even if build fails
   APPVEYOR_SAVE_CACHE_ON_ERROR: true
+  PYTHON: "C:\\Python37-x64"
+  CMAKE_GMP_FLAGS: -DGMP_INCLUDE_DIR="c:/Tools/vcpkg/installed/x64-windows/include" -DGMP_LIBRARIES="c:/Tools/vcpkg/installed/x64-windows/lib/mpir.lib"
+  CMAKE_MPFR_FLAGS: -DMPFR_INCLUDE_DIR="c:/Tools/vcpkg/installed/x64-windows/include" -DMPFR_LIBRARIES="c:/Tools/vcpkg/installed/x64-windows/lib/mpfr.lib"
+  CMAKE_VCPKG_FLAGS: -DCMAKE_TOOLCHAIN_FILE=c:/Tools/vcpkg/scripts/buildsystems/vcpkg.cmake
 
   matrix:
       - target: Examples
         CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=ON -DWITH_GUDHI_TEST=OFF -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=OFF
-        PYTHON: "C:\\Python37-x64"
 
       - target: UnitaryTests
         CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=OFF
-        PYTHON: "C:\\Python37-x64"
 
       - target: Utilities
         CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_TEST=OFF -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=OFF
-        PYTHON: "C:\\Python37-x64"
 
       - target: Python
         CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_TEST=OFF -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=ON
-        PYTHON: "C:\\Python37-x64"
 
 
 cache:
@@ -37,17 +37,17 @@ cache:
 init:
   - echo %target%
 
-
+# tbb:x64-windows 
 install:
   - git submodule update --init
-  - vcpkg install tbb:x64-windows boost-disjoint-sets:x64-windows boost-serialization:x64-windows boost-date-time:x64-windows boost-system:x64-windows boost-filesystem:x64-windows boost-units:x64-windows boost-thread:x64-windows boost-program-options:x64-windows eigen3:x64-windows mpfr:x64-windows mpir:x64-windows cgal:x64-windows
+  - vcpkg install boost-disjoint-sets:x64-windows boost-serialization:x64-windows boost-date-time:x64-windows boost-system:x64-windows boost-filesystem:x64-windows boost-units:x64-windows boost-thread:x64-windows boost-program-options:x64-windows eigen3:x64-windows mpfr:x64-windows mpir:x64-windows cgal:x64-windows
   - SET PATH=c:\Tools\vcpkg\installed\x64-windows\bin;%PATH%
   - SET PATH=%PYTHON%;%PYTHON%\Scripts;%PYTHON%\Library\bin;%PATH%
   - SET PYTHONPATH=%PYTHON%\\Lib\\site-packages;%PYTHONPATH%
   - CALL "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" amd64
   - python --version
   - pip --version
-  - python -m pip install --upgrade pip
+  - python -m pip install --user --upgrade pip
   - python -m pip install --user -r .github/build-requirements.txt
   # No PyKeOps on windows, let's workaround this one.
   - for /F "tokens=*" %%A in (.github/test-requirements.txt) do python -m pip install --user  %%A
@@ -55,9 +55,10 @@ install:
 build_script:
   - mkdir build
   - cd build
-  - cmake -G "Visual Studio 15 2017 Win64" %CMAKE_FLAGS% -DGMP_INCLUDE_DIR="c:/Tools/vcpkg/installed/x64-windows/include" -DGMP_LIBRARIES="c:/Tools/vcpkg/installed/x64-windows/lib/mpir.lib" -DGMP_LIBRARIES_DIR="c:/Tools/vcpkg/installed/x64-windows/lib" -DCMAKE_TOOLCHAIN_FILE=c:/Tools/vcpkg/scripts/buildsystems/vcpkg.cmake ..
+  - cmake -G "Visual Studio 15 2017 Win64" %CMAKE_FLAGS% %CMAKE_GMP_FLAGS% %CMAKE_MPFR_FLAGS% %CMAKE_VCPKG_FLAGS% ..
   - if [%target%]==[Python] (
       cd src/python &
+      type setup.py &
       MSBuild Cython.sln /m /p:Configuration=Release /p:Platform=x64 &
       ctest -j 1 --output-on-failure -C Release
     ) else (
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 40ddc08e..b04efd52 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -61,7 +61,7 @@ jobs:
             cd build;            
             cmake -DCMAKE_BUILD_TYPE=Release -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=ON -DPython_ADDITIONAL_VERSIONS=3 ..;
             cd python;
-            python3 setup.py build_ext -j 2 --inplace;
+            python3 setup.py build_ext --inplace;
             make sphinx;
             cp -R sphinx /tmp/sphinx;
             python3 setup.py install;
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index c19ebb79..f56e12d0 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -14,6 +14,9 @@
 
 #include <boost/version.hpp>
 #include <boost/variant.hpp>
+#include <boost/range/size.hpp>
+#include <boost/range/combine.hpp>
+#include <boost/range/adaptor/transformed.hpp>
 
 #include <gudhi/Debug_utils.h>
 #include <gudhi/Alpha_complex_options.h>
@@ -277,8 +280,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
   Alpha_complex_3d(const InputPointRange& points) {
     static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");
 
-    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(
-        new Alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL));
+    alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(
+        std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL);
   }
 
   /** \brief Alpha_complex constructor from a list of points and associated weights.
@@ -299,20 +302,15 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
   Alpha_complex_3d(const InputPointRange& points, WeightRange weights) {
     static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d");
     static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d");
-    GUDHI_CHECK((weights.size() == points.size()),
+    // FIXME: this test is only valid if we have a forward range
+    GUDHI_CHECK(boost::size(weights) == boost::size(points),
                 std::invalid_argument("Points number in range different from weights range number"));
 
-    std::vector<Weighted_point_3> weighted_points_3;
+    auto weighted_points_3 = boost::range::combine(points, weights)
+      | boost::adaptors::transformed([](auto const&t){return Weighted_point_3(boost::get<0>(t), boost::get<1>(t));});
 
-    std::size_t index = 0;
-    weighted_points_3.reserve(points.size());
-    while ((index < weights.size()) && (index < points.size())) {
-      weighted_points_3.push_back(Weighted_point_3(points[index], weights[index]));
-      index++;
-    }
-
-    alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(
-        new Alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL));
+    alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(
+        std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL);
   }
 
   /** \brief Alpha_complex constructor from a list of points and an iso-cuboid coordinates.
@@ -356,7 +354,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
 
     // 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_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL));
+    alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL);
   }
 
   /** \brief Alpha_complex constructor from a list of points, associated weights and an iso-cuboid coordinates.
@@ -388,31 +386,27 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
                    FT z_min, FT x_max, FT y_max, FT z_max) {
     static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d");
     static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d");
-    GUDHI_CHECK((weights.size() == points.size()),
+    // FIXME: this test is only valid if we have a forward range
+    GUDHI_CHECK(boost::size(weights) == boost::size(points),
                 std::invalid_argument("Points number in range different from weights range number"));
     // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it.
     GUDHI_CHECK(
         (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min),
         std::invalid_argument("The size of the cuboid in every directions is not the same."));
 
-    std::vector<Weighted_point_3> weighted_points_3;
-
-    std::size_t index = 0;
-    weighted_points_3.reserve(points.size());
-
 #ifdef GUDHI_DEBUG
     // Defined in GUDHI_DEBUG to avoid unused variable warning for GUDHI_CHECK
     FT maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
 #endif
 
-    while ((index < weights.size()) && (index < points.size())) {
-      GUDHI_CHECK((weights[index] < maximal_possible_weight) && (weights[index] >= 0),
-                  std::invalid_argument("Invalid weight at index " + std::to_string(index + 1) +
-                                        ". Must be positive and less than maximal possible weight = 1/64*cuboid length "
-                                        "squared, which is not an acceptable input."));
-      weighted_points_3.push_back(Weighted_point_3(points[index], weights[index]));
-      index++;
-    }
+    auto weighted_points_3 = boost::range::combine(points, weights)
+      | boost::adaptors::transformed([=](auto const&t){
+          auto w = boost::get<1>(t);
+          GUDHI_CHECK((w < maximal_possible_weight) && (w >= 0),
+              std::invalid_argument("Invalid weight " + std::to_string(w) +
+                ". Must be non-negative and less than maximal possible weight = 1/64*cuboid length squared."));
+          return Weighted_point_3(boost::get<0>(t), w);
+          });
 
     // Define the periodic cube
     Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
@@ -426,7 +420,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
 
     // 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_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL));
+    alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL);
   }
 
   /** \brief Inserts all Delaunay triangulation into the simplicial complex.
@@ -471,6 +465,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
 #ifdef DEBUG_TRACES
     std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
+    if (objects.size() == 0) {
+      std::cerr << "Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane\n";
+      return false;  // ----- >>
+    }
 
     using Alpha_value_iterator = typename std::vector<FT>::const_iterator;
     Alpha_value_iterator alpha_value_iterator = alpha_values.begin();
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 04f64ce3..4f26481e 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -14,6 +14,15 @@ function( add_GUDHI_PYTHON_lib THE_LIB )
   endif(EXISTS ${THE_LIB})
 endfunction( add_GUDHI_PYTHON_lib )
 
+function( add_GUDHI_PYTHON_lib_dir THE_LIB_DIR )
+  # deals when it is not set - error on windows
+  if(EXISTS ${THE_LIB_DIR})
+    set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${THE_LIB_DIR}', " PARENT_SCOPE)
+  else()
+    message("add_GUDHI_PYTHON_lib_dir - '${THE_LIB_DIR}' does not exist")
+  endif()
+endfunction( add_GUDHI_PYTHON_lib_dir )
+
 # THE_TEST is the python test file name (without .py extension) containing tests functions
 function( add_gudhi_py_test THE_TEST )
   if(PYTEST_FOUND)
@@ -154,7 +163,7 @@ if(PYTHONINTERP_FOUND)
       else(CGAL_HEADER_ONLY)
         add_gudhi_debug_info("CGAL version ${CGAL_VERSION}")
         add_GUDHI_PYTHON_lib("${CGAL_LIBRARY}")
-        set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${CGAL_LIBRARIES_DIR}', ")
+        add_GUDHI_PYTHON_lib_dir("${CGAL_LIBRARIES_DIR}")
         message("** Add CGAL ${CGAL_LIBRARIES_DIR}")
         # If CGAL is not header only, CGAL library may link with boost system,
         if(CMAKE_BUILD_TYPE MATCHES Debug)
@@ -162,7 +171,7 @@ if(PYTHONINTERP_FOUND)
         else()
           add_GUDHI_PYTHON_lib("${Boost_SYSTEM_LIBRARY_RELEASE}")
         endif()
-        set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${Boost_LIBRARY_DIRS}', ")
+        add_GUDHI_PYTHON_lib_dir("${Boost_LIBRARY_DIRS}")
         message("** Add Boost ${Boost_LIBRARY_DIRS}")
       endif(CGAL_HEADER_ONLY)
       # GMP and GMPXX are not required, but if present, CGAL will link with them.
@@ -170,13 +179,17 @@ if(PYTHONINTERP_FOUND)
         add_gudhi_debug_info("GMP_LIBRARIES = ${GMP_LIBRARIES}")
         set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMP', ")
         add_GUDHI_PYTHON_lib("${GMP_LIBRARIES}")
-        set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${GMP_LIBRARIES_DIR}', ")
+        if(NOT GMP_LIBRARIES_DIR)
+          get_filename_component(GMP_LIBRARIES_DIR ${GMP_LIBRARIES} PATH)
+          message("GMP_LIBRARIES_DIR from GMP_LIBRARIES set to ${GMP_LIBRARIES_DIR}")
+        endif(NOT GMP_LIBRARIES_DIR)
+        add_GUDHI_PYTHON_lib_dir("${GMP_LIBRARIES_DIR}")
         message("** Add gmp ${GMP_LIBRARIES_DIR}")
         if(GMPXX_FOUND)
           add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}")
           set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ")
           add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}")
-          set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${GMPXX_LIBRARIES_DIR}', ")
+          add_GUDHI_PYTHON_lib_dir("${GMPXX_LIBRARIES_DIR}")
           message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}")
         endif(GMPXX_FOUND)
       endif(GMP_FOUND)
@@ -187,9 +200,10 @@ if(PYTHONINTERP_FOUND)
         # In case CGAL is not header only, all MPFR variables are set except MPFR_LIBRARIES_DIR - Just set it
         if(NOT MPFR_LIBRARIES_DIR)
           get_filename_component(MPFR_LIBRARIES_DIR ${MPFR_LIBRARIES} PATH)
+          message("MPFR_LIBRARIES_DIR from MPFR_LIBRARIES set to ${MPFR_LIBRARIES_DIR}")
         endif(NOT MPFR_LIBRARIES_DIR)
-        set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${MPFR_LIBRARIES_DIR}', ")
-        message("** Add mpfr ${MPFR_LIBRARIES}")
+        add_GUDHI_PYTHON_lib_dir("${MPFR_LIBRARIES_DIR}")
+        message("** Add mpfr ${MPFR_LIBRARIES_DIR}")
       endif(MPFR_FOUND)
     endif(CGAL_FOUND)
 
@@ -216,7 +230,7 @@ if(PYTHONINTERP_FOUND)
         add_GUDHI_PYTHON_lib("${TBB_RELEASE_LIBRARY}")
         add_GUDHI_PYTHON_lib("${TBB_MALLOC_RELEASE_LIBRARY}")
       endif()
-      set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${TBB_LIBRARY_DIRS}', ")
+      add_GUDHI_PYTHON_lib_dir("${TBB_LIBRARY_DIRS}")
       message("** Add tbb ${TBB_LIBRARY_DIRS}")
       set(GUDHI_PYTHON_INCLUDE_DIRS "${GUDHI_PYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ")
     endif()
diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst
index 097470c1..fffcb3db 100644
--- a/src/python/doc/alpha_complex_user.rst
+++ b/src/python/doc/alpha_complex_user.rst
@@ -34,8 +34,8 @@ Remarks
   the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the
   output is a valid filtration (faces have a filtration value no larger than their cofaces).
 * For performances reasons, it is advised to use Alpha_complex with `CGAL <installation.html#cgal>`_ :math:`\geq` 5.0.0.
-
-For performances reasons, it is advised to use CGAL :math:`\geq` 5.0.0.
+* The vertices in the output simplex tree are not guaranteed to match the order of the input points. One can use
+  :func:`~gudhi.AlphaComplex.get_point` to get the initial point back.
 
 Example from points
 -------------------
@@ -178,49 +178,22 @@ In the following example, a threshold of :math:`\alpha^2 = 32.0` is used.
 Example from OFF file
 ^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the Delaunay triangulation from the points given by an OFF file, and initializes the alpha complex
-with it.
-
+This example builds the alpha complex from 300 random points on a 2-torus.
 
-Then, it is asked to display information about the alpha complex:
+Then, it computes the persistence diagram and displays it:
 
-.. testcode::
+.. plot::
+   :include-source:
 
+    import matplotlib.pyplot as plt
     import gudhi
     alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \
-        '/data/points/alphacomplexdoc.off')
-    simplex_tree = alpha_complex.create_simplex_tree(max_alpha_square=32.0)
+        '/data/points/tore3D_300.off')
+    simplex_tree = alpha_complex.create_simplex_tree()
     result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \
         repr(simplex_tree.num_simplices()) + ' simplices - ' + \
         repr(simplex_tree.num_vertices()) + ' vertices.'
     print(result_str)
-    fmt = '%s -> %.2f'
-    for filtered_value in simplex_tree.get_filtration():
-        print(fmt % tuple(filtered_value))
-
-the program output is:
-
-.. testoutput::
-
-   Alpha complex is of dimension 2 - 20 simplices - 7 vertices.
-   [0] -> 0.00
-   [1] -> 0.00
-   [2] -> 0.00
-   [3] -> 0.00
-   [4] -> 0.00
-   [5] -> 0.00
-   [6] -> 0.00
-   [2, 3] -> 6.25
-   [4, 5] -> 7.25
-   [0, 2] -> 8.50
-   [0, 1] -> 9.25
-   [1, 3] -> 10.00
-   [1, 2] -> 11.25
-   [1, 2, 3] -> 12.50
-   [0, 1, 2] -> 13.00
-   [5, 6] -> 13.25
-   [2, 4] -> 20.00
-   [4, 6] -> 22.74
-   [4, 5, 6] -> 22.74
-   [3, 6] -> 30.25
-
+    diag = simplex_tree.persistence()
+    gudhi.plot_persistence_diagram(diag)
+    plt.show()
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index a356384d..ea128743 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -20,6 +20,7 @@ import os
 
 from gudhi.simplex_tree cimport *
 from gudhi.simplex_tree import SimplexTree
+from gudhi import read_points_from_off_file
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
@@ -27,11 +28,9 @@ __license__ = "GPL v3"
 
 cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
     cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface":
-        Alpha_complex_interface(vector[vector[double]] points, bool fast_version) nogil except +
-        # bool from_file is a workaround for cython to find the correct signature
-        Alpha_complex_interface(string off_file, bool fast_version, bool from_file) nogil except +
+        Alpha_complex_interface(vector[vector[double]] points, bool fast_version, bool exact_version) nogil except +
         vector[double] get_point(int vertex) nogil except +
-        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool exact_version, bool default_filtration_value) nogil except +
+        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except +
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
@@ -54,7 +53,6 @@ cdef class AlphaComplex:
     """
 
     cdef Alpha_complex_interface * this_ptr
-    cdef bool exact
 
     # Fake constructor that does nothing but documenting the constructor
     def __init__(self, points=None, off_file='', precision='safe'):
@@ -76,21 +74,20 @@ cdef class AlphaComplex:
     def __cinit__(self, points = None, off_file = '', precision = 'safe'):
         assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'"
         cdef bool fast = precision == 'fast'
-        self.exact = precision == 'exact'
+        cdef bool exact = precision == 'exact'
 
         cdef vector[vector[double]] pts
         if off_file:
             if os.path.isfile(off_file):
-                self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True)
+                points = read_points_from_off_file(off_file = off_file)
             else:
                 print("file " + off_file + " not found.")
-        else:
-            if points is None:
-                # Empty Alpha construction
-                points=[]
-            pts = points
-            with nogil:
-                self.this_ptr = new Alpha_complex_interface(pts, fast)
+        if points is None:
+            # Empty Alpha construction
+            points=[]
+        pts = points
+        with nogil:
+            self.this_ptr = new Alpha_complex_interface(pts, fast, exact)
 
     def __dealloc__(self):
         if self.this_ptr != NULL:
@@ -102,7 +99,7 @@ cdef class AlphaComplex:
         return self.this_ptr != NULL
 
     def get_point(self, vertex):
-        """This function returns the point corresponding to a given vertex.
+        """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`.
 
         :param vertex: The vertex.
         :type vertex: int
@@ -128,5 +125,5 @@ cdef class AlphaComplex:
         cdef bool compute_filtration = default_filtration_value == True
         with nogil:
             self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr,
-                                              mas, self.exact, compute_filtration)
+                                              mas, compute_filtration)
         return stree
diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h
new file mode 100644
index 00000000..d699ad9b
--- /dev/null
+++ b/src/python/include/Alpha_complex_factory.h
@@ -0,0 +1,139 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef INCLUDE_ALPHA_COMPLEX_FACTORY_H_
+#define INCLUDE_ALPHA_COMPLEX_FACTORY_H_
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Alpha_complex_3d.h>
+#include <gudhi/Alpha_complex_options.h>
+#include <CGAL/Epeck_d.h>
+#include <CGAL/Epick_d.h>
+
+#include <boost/range/adaptor/transformed.hpp>
+
+#include "Simplex_tree_interface.h"
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <memory>  // for std::unique_ptr
+
+namespace Gudhi {
+
+namespace alpha_complex {
+
+template <typename CgalPointType>
+std::vector<double> pt_cgal_to_cython(CgalPointType const& point) {
+  std::vector<double> vd;
+  vd.reserve(point.dimension());
+  for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++)
+    vd.push_back(CGAL::to_double(*coord));
+  return vd;
+}
+
+template <typename CgalPointType>
+static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) {
+  return CgalPointType(vec.size(), vec.begin(), vec.end());
+}
+
+class Abstract_alpha_complex {
+ public:
+  virtual std::vector<double> get_point(int vh) = 0;
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) = 0;
+};
+
+class Exact_Alphacomplex_dD : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using Point = typename Kernel::Point_d;
+
+ public:
+  Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    Point const& point = alpha_complex_.get_point(vh);
+    return pt_cgal_to_cython(point);
+  }
+
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  Alpha_complex<Kernel> alpha_complex_;
+};
+
+class Inexact_Alphacomplex_dD : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  using Point = typename Kernel::Point_d;
+
+ public:
+  Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    Point const& point = alpha_complex_.get_point(vh);
+    return pt_cgal_to_cython(point);
+  }
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  Alpha_complex<Kernel> alpha_complex_;
+};
+
+template <complexity Complexity>
+class Alphacomplex_3D : public Abstract_alpha_complex {
+ private:
+  using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3;
+
+  static Point pt_cython_to_cgal_3(std::vector<double> const& vec) {
+    return Point(vec[0], vec[1], vec[2]);
+  }
+
+ public:
+  Alphacomplex_3D(const std::vector<std::vector<double>>& points)
+    : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    Point const& point = alpha_complex_.get_point(vh);
+    return pt_cgal_to_cython(point);
+  }
+
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                           bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square);
+  }
+
+ private:
+  Alpha_complex_3d<Complexity, false, false> alpha_complex_;
+};
+
+
+}  // namespace alpha_complex
+
+}  // namespace Gudhi
+
+#endif  // INCLUDE_ALPHA_COMPLEX_FACTORY_H_
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index 3ac5db1f..23be194d 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -11,12 +11,8 @@
 #ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_H_
 #define INCLUDE_ALPHA_COMPLEX_INTERFACE_H_
 
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Alpha_complex.h>
-#include <CGAL/Epeck_d.h>
-#include <CGAL/Epick_d.h>
-
-#include <boost/range/adaptor/transformed.hpp>
+#include "Alpha_complex_factory.h"
+#include <gudhi/Alpha_complex_options.h>
 
 #include "Simplex_tree_interface.h"
 
@@ -30,67 +26,51 @@ namespace Gudhi {
 namespace alpha_complex {
 
 class Alpha_complex_interface {
- private:
-  using Exact_kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
-  using Inexact_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
-  using Point_exact_kernel = typename Exact_kernel::Point_d;
-  using Point_inexact_kernel = typename Inexact_kernel::Point_d;
-
-  template <typename CgalPointType>
-  std::vector<double> pt_cgal_to_cython(CgalPointType& point) {
-    std::vector<double> vd;
-    for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++)
-      vd.push_back(CGAL::to_double(*coord));
-    return vd;
-  }
-
-  template <typename CgalPointType>
-  static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) {
-    return CgalPointType(vec.size(), vec.begin(), vec.end());
-  }
-
  public:
-  Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version)
-      : fast_version_(fast_version) {
-    if (fast_version_) {
-      ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(
-          boost::adaptors::transform(points, pt_cython_to_cgal<Point_inexact_kernel>));
-    } else {
-      ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(
-          boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>));
-    }
-  }
-
-  Alpha_complex_interface(const std::string& off_file_name, bool fast_version, bool from_file = true)
-      : fast_version_(fast_version) {
-    if (fast_version_)
-      ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(off_file_name);
-    else
-      ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(off_file_name);
+  Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version)
+  : points_(points),
+    fast_version_(fast_version),
+    exact_version_(exact_version) {
   }
 
   std::vector<double> get_point(int vh) {
-    if (fast_version_) {
-      Point_inexact_kernel const& point = ac_inexact_ptr_->get_point(vh);
-      return pt_cgal_to_cython(point);
-    } else {
-      Point_exact_kernel const& point = ac_exact_ptr_->get_point(vh);
-      return pt_cgal_to_cython(point);
-    }
+    return alpha_ptr_->get_point(vh);
   }
 
-  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool exact_version,
+  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
                            bool default_filtration_value) {
-    if (fast_version_)
-      ac_inexact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value);
-    else
-      ac_exact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value);
+    if (points_.size() > 0) {
+      std::size_t dimension = points_[0].size();
+      if (dimension == 3 && !default_filtration_value) {
+        if (fast_version_)
+          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_);
+        else if (exact_version_)
+          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_);
+        else
+          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_);
+        if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) {
+          // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2
+          dimension--;
+          alpha_ptr_.reset();
+        }
+      }
+      // Not ** else ** because we have to take into account if 3d fails
+      if (dimension != 3 || default_filtration_value) {
+        if (fast_version_) {
+          alpha_ptr_ = std::make_unique<Inexact_Alphacomplex_dD>(points_, exact_version_);
+        } else {
+          alpha_ptr_ = std::make_unique<Exact_Alphacomplex_dD>(points_, exact_version_);
+        }
+        alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value);
+      }
+    }
   }
 
  private:
+  std::unique_ptr<Abstract_alpha_complex> alpha_ptr_;
+  std::vector<std::vector<double>> points_;
   bool fast_version_;
-  std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_;
-  std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_;
+  bool exact_version_;
 };
 
 }  // namespace alpha_complex
diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index 943ad2c4..a4ee260b 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -8,7 +8,7 @@
       - YYYY/MM Author: Description of the modification
 """
 
-from gudhi import AlphaComplex, SimplexTree
+import gudhi as gd
 import math
 import numpy as np
 import pytest
@@ -25,17 +25,16 @@ __license__ = "MIT"
 
 
 def _empty_alpha(precision):
-    alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision)
+    alpha_complex = gd.AlphaComplex(points=[[0, 0]], precision = precision)
     assert alpha_complex.__is_defined() == True
 
 def test_empty_alpha():
-    _empty_alpha('fast')
-    _empty_alpha('safe')
-    _empty_alpha('exact')
+    for precision in ['fast', 'safe', 'exact']:
+        _empty_alpha(precision)
 
 def _infinite_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    alpha_complex = AlphaComplex(points=point_list, precision = precision)
+    alpha_complex = gd.AlphaComplex(points=point_list, precision = precision)
     assert alpha_complex.__is_defined() == True
 
     simplex_tree = alpha_complex.create_simplex_tree()
@@ -84,13 +83,12 @@ def _infinite_alpha(precision):
         assert False
 
 def test_infinite_alpha():
-    _infinite_alpha('fast')
-    _infinite_alpha('safe')
-    _infinite_alpha('exact')
+    for precision in ['fast', 'safe', 'exact']:
+        _infinite_alpha(precision)
 
 def _filtered_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
+    filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision)
 
     simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25)
 
@@ -128,9 +126,8 @@ def _filtered_alpha(precision):
     assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)]
 
 def test_filtered_alpha():
-    _filtered_alpha('fast')
-    _filtered_alpha('safe')
-    _filtered_alpha('exact')
+    for precision in ['fast', 'safe', 'exact']:
+        _filtered_alpha(precision)
 
 def _safe_alpha_persistence_comparison(precision):
     #generate periodic signal
@@ -144,10 +141,10 @@ def _safe_alpha_persistence_comparison(precision):
     embedding2 = [[signal[i], delayed[i]] for i in range(len(time))]
     
     #build alpha complex and simplex tree
-    alpha_complex1 = AlphaComplex(points=embedding1, precision = precision)
+    alpha_complex1 = gd.AlphaComplex(points=embedding1, precision = precision)
     simplex_tree1 = alpha_complex1.create_simplex_tree()
     
-    alpha_complex2 = AlphaComplex(points=embedding2, precision = precision)
+    alpha_complex2 = gd.AlphaComplex(points=embedding2, precision = precision)
     simplex_tree2 = alpha_complex2.create_simplex_tree()
     
     diag1 = simplex_tree1.persistence()
@@ -165,7 +162,7 @@ def test_safe_alpha_persistence_comparison():
 
 def _delaunay_complex(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
+    filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision)
 
     simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True)
 
@@ -197,6 +194,54 @@ def _delaunay_complex(precision):
         assert math.isnan(filtered_value[1])
 
 def test_delaunay_complex():
-    _delaunay_complex('fast')
-    _delaunay_complex('safe')
-    _delaunay_complex('exact')
+    for precision in ['fast', 'safe', 'exact']:
+        _delaunay_complex(precision)
+
+def _3d_points_on_a_plane(precision, default_filtration_value):
+    alpha = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/alphacomplexdoc.off',
+                         precision = precision)
+
+    simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value)
+    assert simplex_tree.dimension() == 2
+    assert simplex_tree.num_vertices() == 7
+    assert simplex_tree.num_simplices() == 25
+
+def test_3d_points_on_a_plane():
+    for default_filtration_value in [True, False]:
+        for precision in ['fast', 'safe', 'exact']:
+            _3d_points_on_a_plane(precision, default_filtration_value)
+
+def _3d_tetrahedrons(precision):
+    points = 10*np.random.rand(10, 3)
+    alpha = gd.AlphaComplex(points=points, precision = precision)
+    st_alpha = alpha.create_simplex_tree(default_filtration_value = False)
+    # New AlphaComplex for get_point to work
+    delaunay = gd.AlphaComplex(points=points, precision = precision)
+    st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True)
+
+    delaunay_tetra = []
+    for sk in st_delaunay.get_skeleton(4):
+        if len(sk[0]) == 4:
+            tetra = [delaunay.get_point(sk[0][0]),
+                     delaunay.get_point(sk[0][1]),
+                     delaunay.get_point(sk[0][2]),
+                     delaunay.get_point(sk[0][3]) ]
+            delaunay_tetra.append(sorted(tetra, key=lambda tup: tup[0]))
+
+    alpha_tetra = []
+    for sk in st_alpha.get_skeleton(4):
+        if len(sk[0]) == 4:
+            tetra = [alpha.get_point(sk[0][0]),
+                     alpha.get_point(sk[0][1]),
+                     alpha.get_point(sk[0][2]),
+                     alpha.get_point(sk[0][3]) ]
+            alpha_tetra.append(sorted(tetra, key=lambda tup: tup[0]))
+
+    # Check the tetrahedrons from one list are in the second one
+    assert len(alpha_tetra) == len(delaunay_tetra)
+    for tetra_from_del in delaunay_tetra:
+        assert tetra_from_del in alpha_tetra
+
+def test_3d_tetrahedrons():
+    for precision in ['fast', 'safe', 'exact']:
+        _3d_tetrahedrons(precision)