diff options
| -rw-r--r-- | src/Alpha_complex/example/Alpha_complex_from_off.cpp | 25 | ||||
| -rw-r--r-- | src/Alpha_complex/example/Alpha_complex_from_points.cpp | 51 | ||||
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 272 | ||||
| -rw-r--r-- | src/Alpha_complex/test/Alpha_complex_unit_test.cpp | 174 |
4 files changed, 326 insertions, 196 deletions
diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp index a17155de..3ce4c7f4 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp @@ -1,37 +1,40 @@ -// to construct a Delaunay_triangulation from a OFF file -#include "gudhi/Delaunay_triangulation_off_io.h" -#include "gudhi/Alpha_complex.h" - #include <stdio.h> #include <stdlib.h> + #include <string> +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_complex.h" + void usage(char * const progName) { - std::cerr << "Usage: " << progName << " filename.off" << std::endl; - exit(-1); // ----- >> + std::cerr << "Usage: " << progName << " filename.off alpha_square_max_value" << std::endl; + std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0" << std::endl; + exit(-1); // ----- >> } int main(int argc, char **argv) { - if (argc != 2) { + if (argc != 3) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; usage(argv[0]); } std::string off_file_name(argv[1]); + double alpha_square_max_value = atof(argv[2]); // ---------------------------------------------------------------------------- // Init of an alpha complex from an OFF file // ---------------------------------------------------------------------------- typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel; - Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name); - + Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name, alpha_square_max_value); + // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- std::cout << "Alpha complex is of dimension " << alpha_complex_from_file.dimension() << " - " << alpha_complex_from_file.num_simplices() << " simplices - " << alpha_complex_from_file.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_file.filtration_simplex_range()) { std::cout << " ( "; @@ -42,4 +45,4 @@ int main(int argc, char **argv) { std::cout << std::endl; } return 0; -}
\ No newline at end of file +} diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp index 33680a40..b160d702 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp @@ -1,55 +1,40 @@ -// to construct a Delaunay_triangulation from a OFF file -#include "gudhi/Delaunay_triangulation_off_io.h" -#include "gudhi/Alpha_complex.h" - +#include <stdio.h> +#include <stdlib.h> #include <CGAL/Delaunay_triangulation.h> #include <CGAL/Epick_d.h> -#include <stdio.h> -#include <stdlib.h> #include <string> +#include <vector> + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_complex.h" 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); + + std::vector<double> coords = { 0.0, 0.0, 0.0, 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); + coords = { 0.0, 0.0, 1.0, 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); + coords = { 0.0, 1.0, 0.0, 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); + coords = { 1.0, 0.0, 0.0, 0.0 }; points.push_back(Point(coords.begin(), coords.end())); - + // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points // ---------------------------------------------------------------------------- - Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_points(3, points.size(), points.begin(), points.end()); + double max_alpha_square_value = 1e10; + Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_points(3, points.size(), points.begin(), points.end(), + max_alpha_square_value); // ---------------------------------------------------------------------------- // Display information about the alpha complex @@ -57,7 +42,7 @@ int main(int argc, char **argv) { 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 << " ( "; @@ -68,4 +53,4 @@ int main(int argc, char **argv) { std::cout << std::endl; } return 0; -}
\ No newline at end of file +} diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 5834e3df..06e69cf3 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -20,8 +20,8 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ -#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ +#ifndef ALPHA_COMPLEX_H_ +#define ALPHA_COMPLEX_H_ // to construct a simplex_tree from Delaunay_triangulation #include <gudhi/graph_simplicial_complex.h> @@ -31,8 +31,6 @@ #include <stdlib.h> #include <math.h> // isnan, fmax -#include <boost/bimap.hpp> - #include <CGAL/Delaunay_triangulation.h> #include <CGAL/Epick_d.h> #include <CGAL/algorithm.h> @@ -40,11 +38,11 @@ #include <CGAL/enum.h> #include <iostream> -#include <iterator> #include <vector> #include <string> #include <limits> // NaN -#include <iterator> // std::iterator +#include <map> +#include <utility> // std::pair namespace Gudhi { @@ -71,10 +69,6 @@ class Alpha_complex : public Simplex_tree<> { // Simplex_result is the type returned from simplex_tree insert function. typedef typename std::pair<Simplex_handle, bool> 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 typename CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation; @@ -92,38 +86,38 @@ class Alpha_complex : public Simplex_tree<> { // size_type type from CGAL. typedef typename Delaunay_triangulation::size_type size_type; - // Boost bimap type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa. - //typedef typename boost::bimap< CGAL_vertex_iterator, Vertex_handle > Bimap_vertex; - //typedef typename Bimap_vertex::value_type value_type; + // Double map type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa. typedef typename std::map< CGAL_vertex_iterator, Vertex_handle > Map_vertex_iterator_to_handle; typedef typename std::map< Vertex_handle, CGAL_vertex_iterator > Map_vertex_handle_to_iterator; private: /** \brief Map to switch from CGAL vertex iterator to simplex tree vertex handle.*/ - Map_vertex_iterator_to_handle vertex_iterator_to_handle; + Map_vertex_iterator_to_handle vertex_iterator_to_handle_; /** \brief Map to switch from simplex tree vertex handle to CGAL vertex iterator.*/ - Map_vertex_handle_to_iterator vertex_handle_to_iterator; + Map_vertex_handle_to_iterator vertex_handle_to_iterator_; /** \brief Pointer on the CGAL Delaunay triangulation.*/ - Delaunay_triangulation* triangulation; - /** \brief Kernel for triangulation functions access.*/ - Kernel kernel; + Delaunay_triangulation* triangulation_; + /** \brief Kernel for triangulation_ functions access.*/ + Kernel kernel_; + /** \brief Maximum value for alpha square.*/ + Filtration_value max_alpha_square_; 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) { + Alpha_complex(const std::string& off_file_name, Filtration_value max_alpha_square) + : triangulation_(nullptr), + max_alpha_square_(max_alpha_square) { Gudhi::Delaunay_triangulation_off_reader<Delaunay_triangulation> 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); // ----- >> + exit(-1); // ----- >> } - triangulation = off_reader.get_complex(); + triangulation_ = off_reader.get_complex(); init(); } @@ -131,8 +125,9 @@ class Alpha_complex : public Simplex_tree<> { * * @param[in] triangulation_ptr Pointer on a Delaunay triangulation. */ - Alpha_complex(Delaunay_triangulation* triangulation_ptr) - : triangulation(triangulation_ptr) { + Alpha_complex(Delaunay_triangulation* triangulation_ptr, Filtration_value max_alpha_square) + : triangulation_(triangulation_ptr), + max_alpha_square_(max_alpha_square) { init(); } @@ -146,13 +141,15 @@ class Alpha_complex : public Simplex_tree<> { * @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); - size_type inserted = triangulation->insert(firstPoint, lastPoint); + Alpha_complex(int dimension, size_type size, ForwardIterator firstPoint, ForwardIterator lastPoint, + Filtration_value max_alpha_square) + : triangulation_(nullptr), + max_alpha_square_(max_alpha_square) { + triangulation_ = new Delaunay_triangulation(dimension); + size_type inserted = triangulation_->insert(firstPoint, lastPoint); if (inserted != size) { std::cerr << "Alpha_complex - insertion failed " << inserted << " != " << size << std::endl; - exit(-1); // ----- >> + exit(-1); // ----- >> } init(); } @@ -162,7 +159,7 @@ class Alpha_complex : public Simplex_tree<> { * @warning Deletes the Delaunay triangulation. */ ~Alpha_complex() { - delete triangulation; + delete triangulation_; } /** \brief get_point returns the point corresponding to the vertex given as parameter. @@ -173,17 +170,16 @@ class Alpha_complex : public Simplex_tree<> { Point_d get_point(Vertex_handle vertex) { Point_d point; try { - if (vertex_handle_to_iterator[vertex] != nullptr) { - point = vertex_handle_to_iterator[vertex]->point(); + if (vertex_handle_to_iterator_[vertex] != nullptr) { + point = vertex_handle_to_iterator_[vertex]->point(); } - } catch (...) { + } catch (...) { std::cerr << "Alpha_complex - getPoint not found on vertex " << vertex << std::endl; } return point; } 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 @@ -192,108 +188,136 @@ class Alpha_complex : public Simplex_tree<> { * Initialization can be launched once. */ void init() { - if (triangulation == nullptr) { + if (triangulation_ == nullptr) { std::cerr << "Alpha_complex init - Cannot init from a NULL triangulation" << std::endl; - return; // ----- >> + return; // ----- >> } - if (triangulation->number_of_vertices() < 1) { + if (triangulation_->number_of_vertices() < 1) { std::cerr << "Alpha_complex init - Cannot init from a triangulation without vertices" << std::endl; - return; // ----- >> + return; // ----- >> } - if (triangulation->maximal_dimension() < 1) { + if (triangulation_->maximal_dimension() < 1) { std::cerr << "Alpha_complex init - Cannot init from a zero-dimension triangulation" << std::endl; - return; // ----- >> + return; // ----- >> } if (num_vertices() > 0) { std::cerr << "Alpha_complex init - Cannot init twice" << std::endl; - return; // ----- >> + return; // ----- >> } - set_dimension(triangulation->maximal_dimension()); + set_dimension(triangulation_->maximal_dimension()); // -------------------------------------------------------------------------------------------- - // bimap to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa + // double map 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) { + for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { #ifdef DEBUG_TRACES std::cout << "Vertex insertion - " << vertex_handle << " -> " << vit->point() << std::endl; #endif // DEBUG_TRACES - vertex_iterator_to_handle[vit] = vertex_handle; - vertex_handle_to_iterator[vertex_handle] = vit; + vertex_iterator_to_handle_[vit] = vertex_handle; + vertex_handle_to_iterator_[vertex_handle] = vit; vertex_handle++; } // -------------------------------------------------------------------------------------------- + Filtration_value filtration_max = 0.0; // -------------------------------------------------------------------------------------------- // 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; + for (auto cit = triangulation_->finite_full_cells_begin(); cit != triangulation_->finite_full_cells_end(); ++cit) { + Vector_vertex vertex_full_cell; #ifdef DEBUG_TRACES std::cout << "Simplex_tree insertion "; #endif // DEBUG_TRACES for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { if (*vit != nullptr) { #ifdef DEBUG_TRACES - std::cout << " " << vertex_iterator_to_handle[*vit]; + std::cout << " " << vertex_iterator_to_handle_[*vit]; #endif // DEBUG_TRACES // Vector of vertex construction for simplex_tree structure - vertexVector.push_back(vertex_iterator_to_handle[*vit]); + vertex_full_cell.push_back(vertex_iterator_to_handle_[*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<double>::quiet_NaN()); + + Simplex_tree<> full_cell; + full_cell.set_dimension(triangulation_->maximal_dimension()); + // Create a simplex tree containing only one of the full cells + Simplex_result insert_result = full_cell.insert_simplex_and_subfaces(vertex_full_cell); if (!insert_result.second) { std::cerr << "Alpha_complex::init insert_simplex_and_subfaces failed" << std::endl; + exit(-1); // ----->> } - } - // -------------------------------------------------------------------------------------------- - - 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; + bool skip_loop = false; + // +++ For i : d -> 0 + // This loop is skipped in case alpha²(Sigma) > max_alpha_square_ + for (int fc_decr_dim = full_cell.dimension(); (fc_decr_dim >= 0) && (!skip_loop); fc_decr_dim--) { + // +++ Foreach Sigma of dim i + // No need to skip this loop in case alpha²(Sigma) > max_alpha_square_ because of + // if (fc_decr_dim == f_simplex_dim) which means "only for a full cell" + for (auto fc_simplex : full_cell.skeleton_simplex_range(fc_decr_dim)) { + int f_simplex_dim = full_cell.dimension(fc_simplex); + if (fc_decr_dim == f_simplex_dim) { + Vector_of_CGAL_points pointVector; + Vector_vertex current_vertex; #ifdef DEBUG_TRACES - std::cout << "Sigma of dim " << decr_dim << " is"; + std::cout << "Sigma of dim " << fc_decr_dim << " is"; #endif // DEBUG_TRACES - for (auto vertex : simplex_vertex_range(f_simplex)) { - pointVector.push_back(get_point(vertex)); + for (auto vertex : full_cell.simplex_vertex_range(fc_simplex)) { + pointVector.push_back(get_point(vertex)); + current_vertex.push_back(vertex); #ifdef DEBUG_TRACES - std::cout << " " << vertex; + std::cout << " " << vertex; #endif // DEBUG_TRACES - } + } #ifdef DEBUG_TRACES - std::cout << std::endl; + 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 - Squared_Radius squared_radius = kernel.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); + Simplex_handle sigma_handle = find(current_vertex); + // +++ If filt(Sigma) is NaN : filt(Sigma) = alpha²(Sigma) + if ((sigma_handle == null_simplex()) || isnan(filtration(sigma_handle))) { + Filtration_value alpha_complex_filtration = compute_alpha_square(pointVector.begin(), pointVector.end(), + f_simplex_dim); + if (alpha_complex_filtration <= max_alpha_square_) { + // Only insert Sigma in Simplex tree if alpha²(Sigma) <= max_alpha_square_ + if (sigma_handle == null_simplex()) { +#ifdef DEBUG_TRACES + std::cout << "Alpha_complex::init Sigma not found" << std::endl; +#endif // DEBUG_TRACES + insert_result = insert_simplex(current_vertex, std::numeric_limits<double>::quiet_NaN()); + if (!insert_result.second) { + std::cerr << "Alpha_complex::init insert_simplex failed" << std::endl; + exit(-1); // ----->> + } + // Sigma is the new inserted simplex handle + sigma_handle = insert_result.first; + } #ifdef DEBUG_TRACES - std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << filtration(f_simplex) << std::endl; + std::cout << "Alpha_complex::init filtration = " << alpha_complex_filtration << std::endl; #endif // DEBUG_TRACES + assign_filtration(sigma_handle, alpha_complex_filtration); + filtration_max = fmax(filtration_max, alpha_complex_filtration); + } else { + // if alpha²(Sigma) > max_alpha_square_ go to the next full cell + skip_loop = true; +#ifdef DEBUG_TRACES + std::cout << "Alpha_complex::init skip loop on this full cell" << std::endl; +#endif // DEBUG_TRACES + break; + } + } // --- If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) + if (filtration(sigma_handle) <= max_alpha_square_) { + // Propagate alpha filtration value in Simplex tree if alpha²(Sigma) <= max_alpha_square_ + // in case Sigma is not found AND not inserted (alpha_complex_filtration > max_alpha_square_), + // filtration(null_simplex()) returns INFINITY => no propagation + propagate_alpha_filtration(full_cell, fc_simplex, fc_decr_dim, sigma_handle); + } } - propagate_alpha_filtration(f_simplex, decr_dim); - } - } + } // --- Foreach Sigma of dim i + } // --- For i : d -> 0 } // -------------------------------------------------------------------------------------------- @@ -303,41 +327,60 @@ class Alpha_complex : public Simplex_tree<> { set_filtration(filtration_max); } - template<typename Simplex_handle> - void propagate_alpha_filtration(Simplex_handle f_simplex, int decr_dim) { + template<typename ForwardIterator > + Filtration_value compute_alpha_square(ForwardIterator firstPoint, ForwardIterator lastPoint, int f_simplex_dim) { + Filtration_value alpha_square_value = 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 + Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); + + alpha_square_value = squared_radius(firstPoint, lastPoint); + } + return alpha_square_value; + } + + void propagate_alpha_filtration(Simplex_tree& full_cell, Simplex_handle fc_simplex, int fc_decr_dim, + Simplex_handle sigma_handle) { // ### Foreach Tau face of Sigma - for (auto f_boundary : boundary_simplex_range(f_simplex)) { + for (auto f_boundary : full_cell.boundary_simplex_range(fc_simplex)) { #ifdef DEBUG_TRACES std::cout << " | --------------------------------------------------" << std::endl; std::cout << " | Tau "; - for (auto vertex : simplex_vertex_range(f_boundary)) { +#endif // DEBUG_TRACES + Vector_vertex tau_vertex; + for (auto vertex : full_cell.simplex_vertex_range(f_boundary)) { + tau_vertex.push_back(vertex); +#ifdef DEBUG_TRACES std::cout << vertex << " "; +#endif // DEBUG_TRACES } +#ifdef DEBUG_TRACES std::cout << "is a face of Sigma" << std::endl; - std::cout << " | isnan(filtration(Tau)=" << isnan(filtration(f_boundary)) << std::endl; #endif // DEBUG_TRACES + Simplex_handle tau_handle = find(tau_vertex); // ### If filt(Tau) is not NaN - if (!isnan(filtration(f_boundary))) { + + if ((tau_handle != null_simplex()) && (!isnan(filtration(tau_handle)))) { // ### 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); + Filtration_value alpha_complex_filtration = fmin(filtration(tau_handle), filtration(sigma_handle)); + assign_filtration(tau_handle, 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; + std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << alpha_complex_filtration << 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) { + if (fc_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)) { + for (auto vertex : full_cell.simplex_vertex_range(f_boundary)) { pointVector.push_back(get_point(vertex)); } // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function - for (auto vertex : simplex_vertex_range(f_simplex)) { + for (auto vertex : simplex_vertex_range(sigma_handle)) { if (std::find(pointVector.begin(), pointVector.end(), get_point(vertex)) == pointVector.end()) { // vertex is not found in Tau vertexForGabriel = vertex; @@ -346,7 +389,7 @@ class Alpha_complex : public Simplex_tree<> { } } // is_gabriel function initialization - Is_Gabriel is_gabriel = kernel.side_of_bounded_sphere_d_object(); + Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), get_point(vertexForGabriel)) != CGAL::ON_BOUNDED_SIDE; #ifdef DEBUG_TRACES @@ -354,23 +397,34 @@ class Alpha_complex : public Simplex_tree<> { #endif // DEBUG_TRACES // ### If Tau is not Gabriel of Sigma if (false == is_gab) { + if (tau_handle == null_simplex()) { +#ifdef DEBUG_TRACES + std::cout << " | Tau not found" << std::endl; +#endif // DEBUG_TRACES + // in case Tau is not yet created + Simplex_result insert_result = insert_simplex(tau_vertex, std::numeric_limits<double>::quiet_NaN()); + if (!insert_result.second) { + std::cerr << "Alpha_complex::propagate_alpha_filtration insert_simplex failed" << std::endl; + exit(-1); // ----->> + } + // Sigma is the new inserted simplex handle + tau_handle = insert_result.first; + } // ### filt(Tau) = filt(Sigma) - Filtration_value alpha_complex_filtration = filtration(f_simplex); - assign_filtration(f_boundary, alpha_complex_filtration); + assign_filtration(tau_handle, filtration(sigma_handle)); // 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; + std::cout << " | filt(Tau) = filt(Sigma) = " << filtration(sigma_handle) << std::endl; #endif // DEBUG_TRACES } } } } } - }; -} // namespace alphacomplex +} // namespace alphacomplex -} // namespace Gudhi +} // namespace Gudhi -#endif // SRC_ALPHA_COMPLEX_INCLUDE_GUDHI_ALPHA_COMPLEX_H_ +#endif // ALPHA_COMPLEX_H_ diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index b7aceb5e..aa9500e7 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -24,15 +24,17 @@ #include <boost/test/included/unit_test.hpp> #include <boost/system/error_code.hpp> #include <boost/chrono/thread_clock.hpp> -// 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 <cmath> // float comparison +#include <cmath> // float comparison #include <limits> +#include <string> +#include <vector> + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_complex.h" // Use dynamic_dimension_tag for the user to be able to set dimension typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d; @@ -45,9 +47,11 @@ 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; + double max_alpha_square_value = 1e10; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; - Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name); + Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name, max_alpha_square_value); const int DIMENSION = 4; std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl; @@ -59,8 +63,51 @@ BOOST_AUTO_TEST_CASE(S4_100_OFF_file) { 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); + // TODO(VR) : BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + + // + TODO(VR) : in wait of num_simplices fix in Simplex_tree [ + int num_simplices = 0; + for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { + num_simplices++; + } + std::cout << "num_simplices=" << num_simplices << std::endl; + BOOST_CHECK(num_simplices == NUMBER_OF_SIMPLICES); + // - TODO(VR) : in wait of num_simplices fix in Simplex_tree ] +} + +BOOST_AUTO_TEST_CASE(S4_100_OFF_file_filtered) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("S4_100.off"); + double max_alpha_square_value = 0.999; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name, max_alpha_square_value); + + 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 = 13; // Versus 100, because of filtered alpha value + 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 = 90; // Versus 6879, because of filtered alpha value + std::cout << "alpha_complex_from_file.num_simplices()=" << alpha_complex_from_file.num_simplices() << std::endl; + // TODO(VR) : BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + // + TODO(VR) : in wait of num_simplices fix in Simplex_tree [ + int num_simplices = 0; + for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { + num_simplices++; + } + std::cout << "num_simplices=" << num_simplices << std::endl; + BOOST_CHECK(num_simplices == NUMBER_OF_SIMPLICES); + // - TODO(VR) : in wait of num_simplices fix in Simplex_tree ] } BOOST_AUTO_TEST_CASE(S8_10_OFF_file) { @@ -70,9 +117,11 @@ 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; + double max_alpha_square_value = 1e10; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; - Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name); + Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name, max_alpha_square_value); const int DIMENSION = 8; std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl; @@ -84,7 +133,51 @@ BOOST_AUTO_TEST_CASE(S8_10_OFF_file) { 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); + // TODO(VR) : BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + + // + TODO(VR) : in wait of num_simplices fix in Simplex_tree [ + int num_simplices = 0; + for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { + num_simplices++; + } + std::cout << "num_simplices=" << num_simplices << std::endl; + BOOST_CHECK(num_simplices == NUMBER_OF_SIMPLICES); + // - TODO(VR) : in wait of num_simplices fix in Simplex_tree ] +} + +BOOST_AUTO_TEST_CASE(S8_10_OFF_file_filtered) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("S8_10.off"); + double max_alpha_square_value = 1.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name, max_alpha_square_value); + + 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 = 895; // Versus 1007, because of filtered alpha value + std::cout << "alpha_complex_from_file.num_simplices()=" << alpha_complex_from_file.num_simplices() << std::endl; + // TODO(VR) : BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + + // + TODO(VR) : in wait of num_simplices fix in Simplex_tree [ + int num_simplices = 0; + for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { + num_simplices++; + } + std::cout << "num_simplices=" << num_simplices << std::endl; + BOOST_CHECK(num_simplices == NUMBER_OF_SIMPLICES); + // - TODO(VR) : in wait of num_simplices fix in Simplex_tree ] } bool are_almost_the_same(float a, float b) { @@ -100,57 +193,53 @@ typedef std::vector<Point> Vector_of_points; bool is_point_in_list(Vector_of_points points_list, Point point) { for (auto& point_in_list : points_list) { if (point_in_list == point) { - return true; // point found + return true; // point found } } - return false; // point not found + return false; // point not found } BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { - // ---------------------------------------------------------------------------- // Init of a list of points // ---------------------------------------------------------------------------- Vector_of_points points; - std::vector<double> coords; - - points.clear(); - coords.clear(); - coords.push_back(0.0); - coords.push_back(0.0); - coords.push_back(0.0); - coords.push_back(1.0); + std::vector<double> coords = { 0.0, 0.0, 0.0, 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); + coords = { 0.0, 0.0, 1.0, 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); + coords = { 0.0, 1.0, 0.0, 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); + coords = { 1.0, 0.0, 0.0, 0.0 }; points.push_back(Point(coords.begin(), coords.end())); // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points // ---------------------------------------------------------------------------- - Gudhi::alphacomplex::Alpha_complex<Kernel_s> alpha_complex_from_points(3, points.size(), points.begin(), points.end()); + double max_alpha_square_value = 1e10; + Gudhi::alphacomplex::Alpha_complex<Kernel_s> alpha_complex_from_points(3, points.size(), points.begin(), points.end(), + max_alpha_square_value); std::cout << "========== Alpha_complex_from_points ==========" << std::endl; + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; // TODO(VR) : in wait of num_simplices fix in Simplex_tree + for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) { + num_simplices++; // TODO(VR) : in wait of num_simplices fix in Simplex_tree + 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; + } + std::cout << "alpha_complex_from_points.dimension()=" << alpha_complex_from_points.dimension() << std::endl; BOOST_CHECK(alpha_complex_from_points.dimension() == 4); - 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); + // TODO(VR) : std::cout << "alpha_complex_from_points.num_simplices()=" << alpha_complex_from_points.num_simplices() + // << std::endl; + // TODO(VR) : BOOST_CHECK(alpha_complex_from_points.num_simplices() == 15); + BOOST_CHECK(num_simplices == 15); // TODO(VR) : in wait of num_simplices fix in Simplex_tree 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); @@ -169,11 +258,11 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { 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 + BOOST_CHECK(false); // Shall not happen break; } } - + Point p1 = alpha_complex_from_points.get_point(1); std::cout << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl; BOOST_CHECK(4 == p1.dimension()); @@ -205,5 +294,4 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { Point p1234 = alpha_complex_from_points.get_point(1234); std::cout << "alpha_complex_from_points.get_point(1234)=" << p1234.dimension() << std::endl; BOOST_CHECK(!is_point_in_list(points, p1234)); - } |