/* This file is part of the Gudhi Library. The Gudhi library
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
* Author(s): Vincent Rouvreau
*
* Copyright (C) 2015 INRIA Saclay (France)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_
#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_
// to construct a simplex_tree from Delaunay_triangulation
#include
#include
#include
#include
#include // isnan, fmax
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include // NaN
#include // std::iterator
namespace Gudhi {
namespace alphacomplex {
#define Kinit(f) =k.f()
/**
* \brief Alpha complex data structure.
*
* \details
* The data structure can be constructed from a CGAL Delaunay triangulation (for more informations on CGAL Delaunay
* triangulation, please refer to the corresponding chapter in page http://doc.cgal.org/latest/Triangulation/) or from
* an OFF file (cf. Delaunay_triangulation_off_reader).
*
* Please refer to \ref alpha_complex for examples.
*
*/
class Alpha_complex : public Simplex_tree<> {
private:
// From Simplex_tree
// Type required to insert into a simplex_tree (with or without subfaces).
typedef std::vector Vector_vertex;
// Simplex_result is the type returned from simplex_tree insert function.
typedef typename std::pair Simplex_result;
// From CGAL
// Kernel for the Delaunay_triangulation. Dimension can be set dynamically.
typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
// Delaunay_triangulation type required to create an alpha-complex.
typedef CGAL::Delaunay_triangulation Delaunay_triangulation;
typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
// Type required to compute squared radius, or side of bounded sphere on a vector of points.
typedef std::vector Vector_of_CGAL_points;
// Vertex_iterator type from CGAL.
typedef Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
// 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.*/
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(const std::string& off_file_name)
: triangulation(nullptr) {
Gudhi::Delaunay_triangulation_off_reader off_reader(off_file_name);
if (!off_reader.is_valid()) {
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) {
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
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(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.
*/
~Alpha_complex() {
delete triangulation;
}
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() {
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
// Start to insert at handle = 0 - default integer value
Vertex_handle vertex_handle = Vertex_handle();
// Loop on triangulation vertices list
for (CGAL_vertex_iterator vit = triangulation->vertices_begin(); vit != triangulation->vertices_end(); ++vit) {
cgal_simplextree.insert(Bimap_vertex::value_type(vit, vertex_handle));
vertex_handle++;
}
// --------------------------------------------------------------------------------------------
// --------------------------------------------------------------------------------------------
// Simplex_tree construction from loop on triangulation finite full cells list
for (auto cit = triangulation->finite_full_cells_begin(); cit != triangulation->finite_full_cells_end(); ++cit) {
Vector_vertex vertexVector;
#ifdef DEBUG_TRACES
std::cout << "Simplex_tree insertion ";
#endif // DEBUG_TRACES
for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
#ifdef DEBUG_TRACES
std::cout << " " << cgal_simplextree.left.at(*vit);
#endif // DEBUG_TRACES
// Vector of vertex construction for simplex_tree structure
vertexVector.push_back(cgal_simplextree.left.at(*vit));
}
#ifdef DEBUG_TRACES
std::cout << std::endl;
#endif // DEBUG_TRACES
// Insert each simplex and its subfaces in the simplex tree - filtration is NaN
Simplex_result insert_result = insert_simplex_and_subfaces(vertexVector,
std::numeric_limits::quiet_NaN());
if (!insert_result.second) {
std::cerr << "Alpha_complex::init insert_simplex_and_subfaces failed" << std::endl;
}
}
// --------------------------------------------------------------------------------------------
Filtration_value filtration_max = 0.0;
// --------------------------------------------------------------------------------------------
// ### For i : d -> 0
for (int decr_dim = dimension(); decr_dim >= 0; decr_dim--) {
// ### Foreach Sigma of dim i
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 : simplex_vertex_range(f_simplex)) {
pointVector.push_back((cgal_simplextree.right.at(vertex))->point());
#ifdef DEBUG_TRACES
std::cout << " " << vertex;
#endif // DEBUG_TRACES
}
#ifdef DEBUG_TRACES
std::cout << std::endl;
#endif // DEBUG_TRACES
// ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
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());
}
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) =" << filtration(f_simplex) << std::endl;
#endif // DEBUG_TRACES
}
propagate_alpha_filtration(f_simplex, decr_dim);
}
}
}
// --------------------------------------------------------------------------------------------
#ifdef DEBUG_TRACES
std::cout << "filtration_max=" << filtration_max << std::endl;
#endif // DEBUG_TRACES
set_filtration(filtration_max);
}
template
void propagate_alpha_filtration(Simplex_handle f_simplex, int decr_dim) {
// ### Foreach Tau face of Sigma
for (auto f_boundary : boundary_simplex_range(f_simplex)) {
#ifdef DEBUG_TRACES
std::cout << " | --------------------------------------------------" << std::endl;
std::cout << " | Tau ";
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(filtration(f_boundary)) << std::endl;
#endif // DEBUG_TRACES
// ### If filt(Tau) is not NaN
if (!isnan(filtration(f_boundary))) {
// ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
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)) = " << filtration(f_boundary) << std::endl;
#endif // DEBUG_TRACES
// ### Else
} else {
// No need to compute is_gabriel for dimension <= 2
// i.e. : Sigma = (3,1) => Tau = 1
if (decr_dim > 1) {
// insert the Tau points in a vector for is_gabriel function
Vector_of_CGAL_points pointVector;
Vertex_handle vertexForGabriel = Vertex_handle();
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 : 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
vertexForGabriel = vertex;
// No need to continue loop
break;
}
}
// is_gabriel function initialization
Kernel k;
Is_Gabriel is_gabriel Kinit(side_of_bounded_sphere_d_object);
#ifdef DEBUG_TRACES
bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertexForGabriel))->point())
!= CGAL::ON_BOUNDED_SIDE;
std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl;
#endif // DEBUG_TRACES
// ### If Tau is not Gabriel of Sigma
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 = 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) = " << filtration(f_boundary) << std::endl;
#endif // DEBUG_TRACES
}
}
}
}
}
};
} // namespace alphacomplex
} // namespace Gudhi
#endif // SRC_ALPHA_COMPLEX_INCLUDE_GUDHI_ALPHA_COMPLEX_H_