Alpha complex documetation - 1st part

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

Former-commit-id: b64e860b0fb6d60e706784605f699349ea17a869

1 files changed, 108 insertions, 99 deletions

diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index d25c05cb..44741e3b 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -52,17 +52,37 @@ namespace alphacomplex {
 #define Kinit(f) =k.f()
 
 /** \defgroup alpha_complex Alpha complex in dimension N
- *
- <DT>Implementations:</DT>
- Alpha complex in dimension N are a subset of Delaunay Triangulation in dimension N.
-
-
- * \author    Vincent Rouvreau
- * \version   1.0
- * \date      2015
- * \copyright GNU General Public License v3.
  * @{
+ * \author    Vincent Rouvreau
+ *
+ * \section Definition
+ * 
+ * Alpha_complex is a Simplex_tree constructed from each finite cell of a Delaunay Triangulation in dimension N.
+ * 
+ * The filtration value of each simplex is computed from the alpha value of the simplex if it is Gabriel or
+ * from the alpha value of the simplex coface that makes the simplex not Gabriel.
+ * 
+ * Please refer to \cite AlphaShapesDefinition for the alpha complex definition or to 
+ * \cite AlphaShapesIntroduction for alpha complex concept vulgarization.
+ * 
+ * \section Example
+ * 
+ * This example loads points from an OFF file, builds the Delaunay triangulation, and finally initialize the 
+ * alpha complex with it.
+ * Then, it is asked to display information about the alpha complex.
+ * 
+ * \include Alpha_complex_from_off.cpp
+ * 
+ * When launching:
+ * 
+ * \code $> ./alphaoffreader ../../data/points/alphashapedoc.off
+ * \endcode
+ *
+ * the program output is:
+ * 
+ * \include alphaoffreader_for_doc.txt
  */
+/** @} */  // end defgroup alpha_complex
 
 /**
  * \brief Alpha complex data structure.
@@ -74,89 +94,107 @@ namespace alphacomplex {
  *
  *
  */
-class Alpha_complex {
+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
-  /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/
+  // Type required to insert into a simplex_tree (with or without subfaces).
   typedef std::vector<Vertex_handle> Vector_vertex;
 
-  /** \brief Simplex_handle type from simplex_tree.*/
-  typedef typename Gudhi::Simplex_tree<>::Simplex_handle Simplex_handle;
-  /** \brief Simplex_result is the type returned from simplex_tree insert function.*/
+  // Simplex_result is the type returned from simplex_tree insert function.
   typedef typename std::pair<Simplex_handle, bool> Simplex_result;
 
-  /** \brief Filtration_simplex_range type from simplex_tree.*/
-  typedef typename Gudhi::Simplex_tree<>::Filtration_simplex_range Filtration_simplex_range;
-
-  /** \brief Simplex_vertex_range type from simplex_tree.*/
-  typedef typename Gudhi::Simplex_tree<>::Simplex_vertex_range Simplex_vertex_range;
-  
   // From CGAL
-  /** \brief Kernel for the Delaunay_triangulation. Dimension can be set dynamically.*/
+  // Kernel for the Delaunay_triangulation. Dimension can be set dynamically.
   typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
-  /** \brief Delaunay_triangulation type required to create an alpha-complex.*/
+
+  // Delaunay_triangulation type required to create an alpha-complex.
   typedef CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation;
 
   typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
   typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
 
-  /** \brief Type required to compute squared radius, or side of bounded sphere on a vector of points.*/
+  // Type required to compute squared radius, or side of bounded sphere on a vector of points.
   typedef std::vector<Kernel::Point_d> Vector_of_CGAL_points;
 
-  /** \brief Vertex_iterator type from CGAL.*/
+  // Vertex_iterator type from CGAL.
   typedef Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
 
-  /** \brief Boost bimap type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.*/
+  // 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;
-  
+
  private:
-  /** \brief Alpha complex is represented internally by a simplex tree.*/
-  Gudhi::Simplex_tree<> st_;
   /** \brief Boost bimap to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.*/
   Bimap_vertex cgal_simplextree;
   /** \brief Pointer on the CGAL Delaunay triangulation.*/
   Delaunay_triangulation* triangulation;
 
  public:
+  /** \brief Alpha_complex constructor from an OFF file name.
+   * Uses the Delaunay_triangulation_off_reader to construct the Delaunay triangulation required to initialize 
+   * the Alpha_complex.
+   *
+   * @param[in] off_file_name OFF file [path and] name.
+   */
   Alpha_complex(std::string& off_file_name)
-  : triangulation(nullptr) {
+      : triangulation(nullptr) {
     Gudhi::Delaunay_triangulation_off_reader<Delaunay_triangulation> off_reader(off_file_name);
     if (!off_reader.is_valid()) {
-      std::cerr << "Unable to read file " << off_file_name << std::endl;
+      std::cerr << "Alpha_complex - Unable to read file " << off_file_name << std::endl;
       exit(-1); // ----- >>
     }
     triangulation = off_reader.get_complex();
     init();
   }
 
+  /** \brief Alpha_complex constructor from a Delaunay triangulation.
+   *
+   * @param[in] triangulation_ptr Pointer on a Delaunay triangulation.
+   */
   Alpha_complex(Delaunay_triangulation* triangulation_ptr)
-  : triangulation(triangulation_ptr) {
+      : triangulation(triangulation_ptr) {
     init();
   }
 
+  /** \brief Alpha_complex destructor from a Delaunay triangulation.
+   *
+   * @warning Deletes the Delaunay triangulation.
+   */
   ~Alpha_complex() {
     delete triangulation;
   }
 
-  Filtration_simplex_range filtration_simplex_range() {
-    return st_.filtration_simplex_range();
-  }
-  
-  Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) {
-    return st_.simplex_vertex_range(sh);
-  }
-  
-  /** \brief Returns the filtration value of a simplex.
-   *
-   * Called on the null_simplex, returns INFINITY. */
-  Gudhi::Simplex_tree<>::Filtration_value filtration(Simplex_handle sh) {
-    return st_.filtration(sh);
-  }
-
  private:
-
+  /** \brief Initialize the Alpha_complex from the Delaunay triangulation.
+   *
+   * @warning Delaunay triangulation must be already constructed with at least one vertex and dimension must be more 
+   * than 0.
+   * 
+   * Initialization can be launched once.
+   */
   void init() {
-    st_.set_dimension(triangulation->maximal_dimension());
+    if (triangulation == nullptr) {
+      std::cerr << "Alpha_complex init - Cannot init from a NULL triangulation" << std::endl;
+      return; // ----- >>
+    }
+    if (triangulation->number_of_vertices() < 1) {
+      std::cerr << "Alpha_complex init - Cannot init from a triangulation without vertices" << std::endl;
+      return; // ----- >>
+    }
+    if (triangulation->maximal_dimension() < 1) {
+      std::cerr << "Alpha_complex init - Cannot init from a zero-dimension triangulation" << std::endl;
+      return; // ----- >>
+    }
+    if (num_vertices() > 0) {
+      std::cerr << "Alpha_complex init - Cannot init twice" << std::endl;
+      return; // ----- >>
+    }
+    
+    set_dimension(triangulation->maximal_dimension());
 
     // --------------------------------------------------------------------------------------------
     // bimap to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa
@@ -187,7 +225,7 @@ class Alpha_complex {
       std::cout << std::endl;
 #endif  // DEBUG_TRACES
       // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
-      Simplex_result insert_result = st_.insert_simplex_and_subfaces(vertexVector,
+      Simplex_result insert_result = insert_simplex_and_subfaces(vertexVector,
                                                                      std::numeric_limits<double>::quiet_NaN());
       if (!insert_result.second) {
         std::cerr << "Alpha_complex::init insert_simplex_and_subfaces failed" << std::endl;
@@ -198,16 +236,16 @@ class Alpha_complex {
     Filtration_value filtration_max = 0.0;
     // --------------------------------------------------------------------------------------------
     // ### For i : d -> 0
-    for (int decr_dim = st_.dimension(); decr_dim >= 0; decr_dim--) {
+    for (int decr_dim = dimension(); decr_dim >= 0; decr_dim--) {
       // ### Foreach Sigma of dim i
-      for (auto f_simplex : st_.skeleton_simplex_range(decr_dim)) {
-        int f_simplex_dim = st_.dimension(f_simplex);
+      for (auto f_simplex : skeleton_simplex_range(decr_dim)) {
+        int f_simplex_dim = dimension(f_simplex);
         if (decr_dim == f_simplex_dim) {
           Vector_of_CGAL_points pointVector;
 #ifdef DEBUG_TRACES
           std::cout << "Sigma of dim " << decr_dim << " is";
 #endif  // DEBUG_TRACES
-          for (auto vertex : st_.simplex_vertex_range(f_simplex)) {
+          for (auto vertex : simplex_vertex_range(f_simplex)) {
             pointVector.push_back((cgal_simplextree.right.at(vertex))->point());
 #ifdef DEBUG_TRACES
             std::cout << " " << vertex;
@@ -217,20 +255,20 @@ class Alpha_complex {
           std::cout << std::endl;
 #endif  // DEBUG_TRACES
           // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
-          if (isnan(st_.filtration(f_simplex))) {
+          if (isnan(filtration(f_simplex))) {
             Filtration_value alpha_complex_filtration = 0.0;
             // No need to compute squared_radius on a single point - alpha is 0.0
             if (f_simplex_dim > 0) {
               // squared_radius function initialization
               Kernel k;
               Squared_Radius squared_radius Kinit(compute_squared_radius_d_object);
-              
+
               alpha_complex_filtration = squared_radius(pointVector.begin(), pointVector.end());
             }
-            st_.assign_filtration(f_simplex, alpha_complex_filtration);
+            assign_filtration(f_simplex, alpha_complex_filtration);
             filtration_max = fmax(filtration_max, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
-            std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << st_.filtration(f_simplex) << std::endl;
+            std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << filtration(f_simplex) << std::endl;
 #endif  // DEBUG_TRACES
           }
           propagate_alpha_filtration(f_simplex, decr_dim);
@@ -242,30 +280,30 @@ class Alpha_complex {
 #ifdef DEBUG_TRACES
     std::cout << "filtration_max=" << filtration_max << std::endl;
 #endif  // DEBUG_TRACES
-    st_.set_filtration(filtration_max);
+    set_filtration(filtration_max);
   }
 
   template<typename Simplex_handle>
   void propagate_alpha_filtration(Simplex_handle f_simplex, int decr_dim) {
     // ### Foreach Tau face of Sigma
-    for (auto f_boundary : st_.boundary_simplex_range(f_simplex)) {
+    for (auto f_boundary : boundary_simplex_range(f_simplex)) {
 #ifdef DEBUG_TRACES
       std::cout << " | --------------------------------------------------" << std::endl;
       std::cout << " | Tau ";
-      for (auto vertex : st_.simplex_vertex_range(f_boundary)) {
+      for (auto vertex : simplex_vertex_range(f_boundary)) {
         std::cout << vertex << " ";
       }
       std::cout << "is a face of Sigma" << std::endl;
-      std::cout << " | isnan(filtration(Tau)=" << isnan(st_.filtration(f_boundary)) << std::endl;
+      std::cout << " | isnan(filtration(Tau)=" << isnan(filtration(f_boundary)) << std::endl;
 #endif  // DEBUG_TRACES
       // ### If filt(Tau) is not NaN
-      if (!isnan(st_.filtration(f_boundary))) {
+      if (!isnan(filtration(f_boundary))) {
         // ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
-        Filtration_value alpha_complex_filtration = fmin(st_.filtration(f_boundary), st_.filtration(f_simplex));
-        st_.assign_filtration(f_boundary, alpha_complex_filtration);
+        Filtration_value alpha_complex_filtration = fmin(filtration(f_boundary), filtration(f_simplex));
+        assign_filtration(f_boundary, alpha_complex_filtration);
         // No need to check for filtration_max, alpha_complex_filtration is a min of an existing filtration value
 #ifdef DEBUG_TRACES
-        std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << st_.filtration(f_boundary) << std::endl;
+        std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << filtration(f_boundary) << std::endl;
 #endif  // DEBUG_TRACES
         // ### Else
       } else {
@@ -275,11 +313,11 @@ class Alpha_complex {
           // insert the Tau points in a vector for is_gabriel function
           Vector_of_CGAL_points pointVector;
           Vertex_handle vertexForGabriel = Vertex_handle();
-          for (auto vertex : st_.simplex_vertex_range(f_boundary)) {
+          for (auto vertex : simplex_vertex_range(f_boundary)) {
             pointVector.push_back((cgal_simplextree.right.at(vertex))->point());
           }
           // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function
-          for (auto vertex : st_.simplex_vertex_range(f_simplex)) {
+          for (auto vertex : simplex_vertex_range(f_simplex)) {
             if (std::find(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertex))->point())
                 == pointVector.end()) {
               // vertex is not found in Tau
@@ -300,46 +338,17 @@ class Alpha_complex {
           if ((is_gabriel(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertexForGabriel))->point())
                == CGAL::ON_BOUNDED_SIDE)) {
             // ### filt(Tau) = filt(Sigma)
-            Filtration_value alpha_complex_filtration = st_.filtration(f_simplex);
-            st_.assign_filtration(f_boundary, alpha_complex_filtration);
+            Filtration_value alpha_complex_filtration = filtration(f_simplex);
+            assign_filtration(f_boundary, alpha_complex_filtration);
             // No need to check for filtration_max, alpha_complex_filtration is an existing filtration value
 #ifdef DEBUG_TRACES
-            std::cout << " | filt(Tau) = filt(Sigma) = " << st_.filtration(f_boundary) << std::endl;
+            std::cout << " | filt(Tau) = filt(Sigma) = " << filtration(f_boundary) << std::endl;
 #endif  // DEBUG_TRACES
           }
         }
       }
     }
   }
- public:
-
-  /** \brief Returns the number of vertices in the complex. */
-  size_t num_vertices() {
-    return st_.num_vertices();
-  }
-
-  /** \brief Returns the number of simplices in the complex.
-   *
-   * Does not count the empty simplex. */
-  const unsigned int& num_simplices() const {
-    return st_.num_simplices();
-  }
-
-  /** \brief Returns an upper bound on the dimension of the simplicial complex. */
-  int dimension() {
-    return st_.dimension();
-  }
-
-  /** \brief Returns an upper bound of the filtration values of the simplices. */
-  Filtration_value filtration() {
-    return st_.filtration();
-  }
-
-  friend std::ostream& operator<<(std::ostream& os, const Alpha_complex & alpha_complex) {
-    Gudhi::Simplex_tree<> st = alpha_complex.st_;
-    os << st << std::endl;
-    return os;
-  }
 };
 
 } // namespace alphacomplex