Alpha complex construction from a list of CGAL points

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

Former-commit-id: 521269f793f9c16c2305db8c97678bea2bf95092

6 files changed, 192 insertions, 20 deletions

diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
index 0d7af117..ce278419 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
@@ -22,7 +22,7 @@ int main(int argc, char **argv) {
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from an OFF file
   // ----------------------------------------------------------------------------
-  Gudhi::alphacomplex::Alpha_complex<> alpha_complex_from_file(off_file_name);
+  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
   
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
new file mode 100644
index 00000000..fc0e2460
--- /dev/null
+++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
@@ -0,0 +1,71 @@
+// to construct a Delaunay_triangulation from a OFF file
+#include "gudhi/Delaunay_triangulation_off_io.h"
+#include "gudhi/Alpha_complex.h"
+
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+typedef Kernel::Point_d Point;
+typedef std::vector<Point> Vector_of_points;
+
+int main(int argc, char **argv) {
+
+  // ----------------------------------------------------------------------------
+  // Init of a list of points
+  // ----------------------------------------------------------------------------
+  Vector_of_points points;
+  std::vector<double> coords;
+  
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from the list of points
+  // ----------------------------------------------------------------------------
+  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_points(3, points.size(), points.begin(), points.end());
+
+  // ----------------------------------------------------------------------------
+  // Display information about the alpha complex
+  // ----------------------------------------------------------------------------
+  std::cout << "Alpha complex is of dimension " << alpha_complex_from_points.dimension() <<
+      " - " << alpha_complex_from_points.num_simplices() << " simplices - " <<
+      alpha_complex_from_points.num_vertices() << " vertices." << std::endl;
+  
+  std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+  for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) {
+    std::cout << "   ( ";
+    for (auto vertex : alpha_complex_from_points.simplex_vertex_range(f_simplex)) {
+      std::cout << vertex << " ";
+    }
+    std::cout << ") -> " << "[" << alpha_complex_from_points.filtration(f_simplex) << "] ";
+    std::cout << std::endl;
+  }
+  return 0;
+}
\ No newline at end of file
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index 04c0ba58..2e64e4db 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -17,6 +17,9 @@ if(CGAL_FOUND)
       #add_definitions(-DDEBUG_TRACES)
       add_executable ( alphaoffreader Alpha_complex_from_off.cpp )
       target_link_libraries(alphaoffreader  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
+
+      add_executable ( alphapoints Alpha_complex_from_points.cpp )
+      target_link_libraries(alphapoints  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
     else()
       message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
     endif()
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 97c30abb..138270ff 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -44,6 +44,7 @@
 #include <vector>
 #include <string>
 #include <limits> // NaN
+#include <iterator>     // std::iterator
 
 namespace Gudhi {
 
@@ -62,11 +63,6 @@ namespace alphacomplex {
  * Please refer to \ref alpha_complex for examples.
  *
  */
-template<typename IndexingTag = linear_indexing_tag,
-typename FiltrationValue = double,
-typename SimplexKey = int,  // must be a signed integer type
-typename VertexHandle = int  // must be a signed integer type, int convertible to it
->
 class Alpha_complex : public Simplex_tree<> {
  private:
   // From Simplex_tree
@@ -94,6 +90,9 @@ class Alpha_complex : public Simplex_tree<> {
 
   // Boost bimap type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.
   typedef boost::bimap< CGAL_vertex_iterator, Vertex_handle > Bimap_vertex;
+  
+  // size_type type from CGAL.
+  typedef Delaunay_triangulation::size_type size_type;
 
  private:
   /** \brief Boost bimap to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.*/
@@ -108,7 +107,7 @@ class Alpha_complex : public Simplex_tree<> {
    *
    * @param[in] off_file_name OFF file [path and] name.
    */
-  Alpha_complex(std::string& off_file_name)
+  Alpha_complex(const std::string& off_file_name)
       : triangulation(nullptr) {
     Gudhi::Delaunay_triangulation_off_reader<Delaunay_triangulation> off_reader(off_file_name);
     if (!off_reader.is_valid()) {
@@ -128,6 +127,27 @@ class Alpha_complex : public Simplex_tree<> {
     init();
   }
 
+  /** \brief Alpha_complex constructor from a list of points.
+   * Uses the Delaunay_triangulation_off_reader to construct the Delaunay triangulation required to initialize 
+   * the Alpha_complex.
+   *
+   * @param[in] dimension Dimension of points to be inserted.
+   * @param[in] size Number of points to be inserted.
+   * @param[in] firstPoint Iterator on the first point to be inserted.
+   * @param[in] last Point Iterator on the last point to be inserted.
+   */
+  template<typename ForwardIterator >
+  Alpha_complex(int dimension, size_type size, ForwardIterator firstPoint, ForwardIterator lastPoint)
+      : triangulation(nullptr) {
+    triangulation = new Delaunay_triangulation(dimension);
+    Delaunay_triangulation::size_type inserted = triangulation->insert<ForwardIterator>(firstPoint, lastPoint);
+    if (inserted != size) {
+      std::cerr << "Alpha_complex - insertion failed " << inserted << " != " << size<< std::endl;
+      exit(-1); // ----- >>
+    }
+    init();
+  }
+
   /** \brief Alpha_complex destructor from a Delaunay triangulation.
    *
    * @warning Deletes the Delaunay triangulation.
diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
index 86d4d9c3..9530314c 100644
--- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
@@ -28,15 +28,19 @@
 #include "gudhi/Delaunay_triangulation_off_io.h"
 #include "gudhi/Alpha_complex.h"
 
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+
 #include <cmath> // float comparison
+#include <limits>
 
 // Use dynamic_dimension_tag for the user to be able to set dimension
-typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K;
-typedef CGAL::Delaunay_triangulation<K> T;
-// The triangulation uses the default instantiation of the 
-// TriangulationDataStructure template parameter
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+typedef Kernel::Point_d Point;
+typedef std::vector<Point> Vector_of_points;
+// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
 
-BOOST_AUTO_TEST_CASE( S4_100_OFF_file ) {
+BOOST_AUTO_TEST_CASE(S4_100_OFF_file) {
   // ----------------------------------------------------------------------------
   //
   // Init of an alpha-complex from a OFF file
@@ -44,24 +48,24 @@ BOOST_AUTO_TEST_CASE( S4_100_OFF_file ) {
   // ----------------------------------------------------------------------------
   std::string off_file_name("S4_100.off");
   std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
-  
+
   Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
 
   const int DIMENSION = 4;
   std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.dimension() == DIMENSION);
-  
+
   const int NUMBER_OF_VERTICES = 100;
   std::cout << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.num_vertices() == NUMBER_OF_VERTICES);
-  
+
   const int NUMBER_OF_SIMPLICES = 6879;
   std::cout << "alpha_complex_from_file.num_simplices()=" << alpha_complex_from_file.num_simplices() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES);
 
 }
 
-BOOST_AUTO_TEST_CASE( S8_10_OFF_file ) {
+BOOST_AUTO_TEST_CASE(S8_10_OFF_file) {
   // ----------------------------------------------------------------------------
   //
   // Init of an alpha-complex from a OFF file
@@ -69,19 +73,91 @@ BOOST_AUTO_TEST_CASE( S8_10_OFF_file ) {
   // ----------------------------------------------------------------------------
   std::string off_file_name("S8_10.off");
   std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
-  
+
   Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
 
   const int DIMENSION = 8;
   std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.dimension() == DIMENSION);
-  
+
   const int NUMBER_OF_VERTICES = 10;
   std::cout << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.num_vertices() == NUMBER_OF_VERTICES);
-  
+
   const int NUMBER_OF_SIMPLICES = 1007;
   std::cout << "alpha_complex_from_file.num_simplices()=" << alpha_complex_from_file.num_simplices() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES);
+}
+
+bool are_almost_the_same(float a, float b) {
+  return std::fabs(a - b) < std::numeric_limits<float>::epsilon();
+}
+
+BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
+
+  // ----------------------------------------------------------------------------
+  // Init of a list of points
+  // ----------------------------------------------------------------------------
+  Vector_of_points points;
+  std::vector<double> coords;
+
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(0.0);
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords.clear();
+  coords.push_back(1.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  coords.push_back(0.0);
+  points.push_back(Point(coords.begin(), coords.end()));
+
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from the list of points
+  // ----------------------------------------------------------------------------
+  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_points(3, points.size(), points.begin(), points.end());
+
+  std::cout << "========== Alpha_complex_from_points ==========" << std::endl;
+
+  std::cout << "alpha_complex_from_points.dimension()=" << alpha_complex_from_points.dimension() << std::endl;
+  BOOST_CHECK(alpha_complex_from_points.dimension() == 3);
+  std::cout << "alpha_complex_from_points.num_simplices()=" << alpha_complex_from_points.num_simplices() << std::endl;
+  BOOST_CHECK(alpha_complex_from_points.num_simplices() == 15);
+  std::cout << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl;
+  BOOST_CHECK(alpha_complex_from_points.num_vertices() == 4);
+
+  for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) {
+    switch (alpha_complex_from_points.dimension(f_simplex)) {
+      case 0:
+        BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 0.0));
+        break;
+      case 1:
+        BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 1.0/2.0));
+        break;
+      case 2:
+        BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 2.0/3.0));
+        break;
+      case 3:
+        BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 3.0/4.0));
+        break;
+      default:
+        BOOST_CHECK(false); // Shall not happen
+        break;
+    }
+  }
 
 }
diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h
index 315aa0ac..770d2216 100644
--- a/src/common/doc/main_page.h
+++ b/src/common/doc/main_page.h
@@ -48,8 +48,10 @@ CGAL is a C++ library which provides easy access to efficient and reliable geome
 The following example requires the <a href="http://www.cgal.org/">Computational Geometry Algorithms Library</a> (CGAL) 
 and will not be built if CGAL is not installed:
   - Simplex_tree/simplex_tree_from_alpha_shapes_3
+  - Alpha_complex/Alpha_complex_from_off
+  - Alpha_complex/Alpha_complex_from_points
 
-Having CGAL version 4.4 or higher installed is recommended. The procedure to install this library according to
+Having CGAL version 4.7 or higher installed is recommended. The procedure to install this library according to
 your operating system is detailed here http://doc.cgal.org/latest/Manual/installation.html
 
 \subsection demos Demos and examples