diff options
Diffstat (limited to 'src/Alpha_complex')
-rw-r--r-- | src/Alpha_complex/example/CMakeLists.txt | 31 | ||||
-rw-r--r-- | src/Alpha_complex/example/Delaunay_triangulation_off_rw.cpp | 79 | ||||
-rw-r--r-- | src/Alpha_complex/example/Simplex_tree_from_delaunay_triangulation.cpp | 77 | ||||
-rwxr-xr-x | src/Alpha_complex/example/alphashapedoc.off | 10 | ||||
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_shapes.h | 311 | ||||
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h | 192 | ||||
-rw-r--r-- | src/Alpha_complex/test/Alpha_shapes_unit_test.cpp | 126 | ||||
-rw-r--r-- | src/Alpha_complex/test/CMakeLists.txt | 31 | ||||
-rw-r--r-- | src/Alpha_complex/test/README | 12 | ||||
-rw-r--r-- | src/Alpha_complex/test/S4_100.off | 102 | ||||
-rw-r--r-- | src/Alpha_complex/test/S8_10.off | 12 |
11 files changed, 983 insertions, 0 deletions
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt new file mode 100644 index 00000000..582a9322 --- /dev/null +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIAlphaShapesExample) + +# need CGAL 4.7 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.7.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + + add_definitions(-DDEBUG_TRACES) + + add_executable ( dtoffrw Delaunay_triangulation_off_rw.cpp ) + target_link_libraries(dtoffrw ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) + + add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp ) + target_link_libraries(stfromdt ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Alpha_complex/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_complex/example/Delaunay_triangulation_off_rw.cpp new file mode 100644 index 00000000..fe889ec0 --- /dev/null +++ b/src/Alpha_complex/example/Delaunay_triangulation_off_rw.cpp @@ -0,0 +1,79 @@ +/* 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) 2014 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 <http://www.gnu.org/licenses/>. + */ + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> +#include <vector> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// 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 instanciation of the +// TriangulationDataStructure template parameter + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " inputFile.off dimension outputFile.off" << std::endl; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if (argc != 4) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + int dimension = 0; + int returnedScanValue = sscanf(argv[2], "%d", &dimension); + if ((returnedScanValue == EOF) || (dimension <= 0)) { + std::cerr << "Error: " << argv[2] << " is not correct" << std::endl; + usage(argv[0]); + } + + T dt(dimension); + std::string offFileName(argv[1]); + Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(offFileName, dt); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offFileName << std::endl; + exit(-1); // ----- >> + } + + std::cout << "number_of_finite_full_cells= " << dt.number_of_finite_full_cells() << std::endl; + + std::string outFileName(argv[3]); + std::string offOutputFile(outFileName); + Gudhi::alphashapes::Delaunay_triangulation_off_writer<T> off_writer(offOutputFile, dt); + + return 0; +}
\ No newline at end of file diff --git a/src/Alpha_complex/example/Simplex_tree_from_delaunay_triangulation.cpp b/src/Alpha_complex/example/Simplex_tree_from_delaunay_triangulation.cpp new file mode 100644 index 00000000..0a24fb56 --- /dev/null +++ b/src/Alpha_complex/example/Simplex_tree_from_delaunay_triangulation.cpp @@ -0,0 +1,77 @@ +/* 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) 2014 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 <http://www.gnu.org/licenses/>. + */ + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_shapes.h" + +// to construct a simplex_tree from Delaunay_triangulation +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Simplex_tree.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// 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 instanciation of the +// TriangulationDataStructure template parameter + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " filename.off" << std::endl; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if (argc != 2) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + std::string off_file_name(argv[1]); + + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a OFF file + // + // ---------------------------------------------------------------------------- + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name); + + std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl; + std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl; + std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl; + std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl; + //std::cout << alpha_shapes_from_file << std::endl; + + return 0; +}
\ No newline at end of file diff --git a/src/Alpha_complex/example/alphashapedoc.off b/src/Alpha_complex/example/alphashapedoc.off new file mode 100755 index 00000000..bb790193 --- /dev/null +++ b/src/Alpha_complex/example/alphashapedoc.off @@ -0,0 +1,10 @@ +nOFF +2 7 0 0 +1.0 1.0 +7.0 0.0 +4.0 6.0 +9.0 6.0 +0.0 14.0 +2.0 19.0 +9.0 17.0 + diff --git a/src/Alpha_complex/include/gudhi/Alpha_shapes.h b/src/Alpha_complex/include/gudhi/Alpha_shapes.h new file mode 100644 index 00000000..f23df51a --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_shapes.h @@ -0,0 +1,311 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ +#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ + +// to construct a Delaunay_triangulation from a OFF file +#include <gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h> + +// to construct a simplex_tree from Delaunay_triangulation +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> + +#include <stdio.h> +#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> +#include <CGAL/assertions.h> +#include <CGAL/enum.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <string> +#include <limits> +#include <map> + +namespace Gudhi { + +namespace alphashapes { + +#define Kinit(f) =k.f() + +/** \defgroup alpha_shapes Alpha shapes in dimension N + * + <DT>Implementations:</DT> + Alpha shapes 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. + * @{ + */ + +/** + * \brief Alpha shapes data structure. + * + * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation + * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$. + * + * Details may be found in \cite boissonnatmariasimplextreealgorithmica. + * + * \implements FilteredComplex + * + */ +class Alpha_shapes { + private: + // From Simplex_tree + /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/ + typedef std::vector<Vertex_handle> typeVectorVertex; + + typedef typename Gudhi::Simplex_tree<>::Simplex_handle Simplex_handle; + typedef typename std::pair<Simplex_handle, bool> Simplex_result; + + // From CGAL + /** \brief 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-shape. + */ + 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 insert into a simplex_tree (with or without subfaces).*/ + typedef std::vector<Kernel::Point_d> typeVectorPoint; + + private: + /** \brief Upper bound on the simplex tree of the simplicial complex.*/ + Gudhi::Simplex_tree<> st_; + + public: + + Alpha_shapes(std::string& off_file_name) { + // Construct a default Delaunay_triangulation (dim=0) - dim will be set in visitor reader init function + Delaunay_triangulation dt(2); + Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation> off_reader(off_file_name, dt); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_name << std::endl; + exit(-1); // ----- >> + } +#ifdef DEBUG_TRACES + std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl; + std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl; + std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl; +#endif // DEBUG_TRACES + init<Delaunay_triangulation>(dt); + } + + template<typename T> + Alpha_shapes(T& triangulation) { + init<T>(triangulation); + } + + ~Alpha_shapes() { } + + private: + + template<typename T> + void init(T& triangulation) { + st_.set_dimension(triangulation.maximal_dimension()); + + // -------------------------------------------------------------------------------------------- + // bimap to retrieve vertex handles from points and vice versa + typedef boost::bimap< Kernel::Point_d, Vertex_handle > bimap_points_vh; + bimap_points_vh points_to_vh; + // Start to insert at handle = 0 - default integer value + Vertex_handle vertex_handle = Vertex_handle(); + // Loop on triangulation vertices list + for (auto vit = triangulation.vertices_begin(); vit != triangulation.vertices_end(); ++vit) { + points_to_vh.insert(bimap_points_vh::value_type(vit->point(), 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) { + typeVectorVertex 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 << " " << points_to_vh.left.at((*vit)->point()); +#endif // DEBUG_TRACES + // Vector of vertex construction for simplex_tree structure + vertexVector.push_back(points_to_vh.left.at((*vit)->point())); + } +#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 = st_.insert_simplex_and_subfaces(vertexVector, + std::numeric_limits<double>::quiet_NaN()); + if (!insert_result.second) { + std::cerr << "Alpha_shapes::init insert_simplex_and_subfaces failed" << std::endl; + } + } + // -------------------------------------------------------------------------------------------- + + Filtration_value filtration_max = 0.0; + + Kernel k; + Squared_Radius squared_radius Kinit(compute_squared_radius_d_object); + Is_Gabriel is_gabriel Kinit(side_of_bounded_sphere_d_object); + // -------------------------------------------------------------------------------------------- + // ### For i : d -> 0 + for (int decr_dim = st_.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); + if (decr_dim == f_simplex_dim) { + typeVectorPoint pointVector; +#ifdef DEBUG_TRACES + std::cout << "Sigma of dim " << decr_dim << " is"; +#endif // DEBUG_TRACES + for (auto vertex : st_.simplex_vertex_range(f_simplex)) { + pointVector.push_back(points_to_vh.right.at(vertex)); +#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(st_.filtration(f_simplex))) { + Filtration_value alpha_shapes_filtration = 0.0; + // No need to compute squared_radius on a single point - alpha is 0.0 + if (f_simplex_dim > 0) { + alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end()); + } + st_.assign_filtration(f_simplex, alpha_shapes_filtration); + filtration_max = fmax(filtration_max, alpha_shapes_filtration); +#ifdef DEBUG_TRACES + std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << st_.filtration(f_simplex) << std::endl; +#endif // DEBUG_TRACES + } + + // ### Foreach Tau face of Sigma + for (auto f_boundary : st_.boundary_simplex_range(f_simplex)) { +#ifdef DEBUG_TRACES + std::cout << " | --------------------------------------------------" << std::endl; + std::cout << " | Tau "; + for (auto vertex : st_.simplex_vertex_range(f_boundary)) { + std::cout << vertex << " "; + } + std::cout << "is a face of Sigma" << std::endl; +#endif // DEBUG_TRACES + // insert the Tau points in a vector for is_gabriel function + typeVectorPoint pointVector; + Vertex_handle vertexForGabriel = Vertex_handle(); + for (auto vertex : st_.simplex_vertex_range(f_boundary)) { + pointVector.push_back(points_to_vh.right.at(vertex)); + } + // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function + for (auto vertex : st_.simplex_vertex_range(f_simplex)) { + if (std::find(pointVector.begin(), pointVector.end(), points_to_vh.right.at(vertex)) == pointVector.end()) { + // vertex is not found in Tau + vertexForGabriel = vertex; + // No need to continue loop + break; + } + } +#ifdef DEBUG_TRACES + std::cout << " | isnan(filtration(Tau)=" << isnan(st_.filtration(f_boundary)) << std::endl; + bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), points_to_vh.right.at(vertexForGabriel)) + != CGAL::ON_BOUNDED_SIDE; + std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; +#endif // DEBUG_TRACES + // ### If filt(Tau) is not NaN + // ### or Tau is not Gabriel of Sigma + if (!isnan(st_.filtration(f_boundary)) || + (is_gabriel(pointVector.begin(), pointVector.end(), points_to_vh.right.at(vertexForGabriel)) == CGAL::ON_BOUNDED_SIDE) + ) { + // ### filt(Tau) = fmin(filt(Tau), filt(Sigma)) + Filtration_value alpha_shapes_filtration = fmin(st_.filtration(f_boundary), st_.filtration(f_simplex)); + st_.assign_filtration(f_boundary, alpha_shapes_filtration); + filtration_max = fmax(filtration_max, alpha_shapes_filtration); +#ifdef DEBUG_TRACES + std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << st_.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + } + } + } + } + } + // -------------------------------------------------------------------------------------------- + +#ifdef DEBUG_TRACES + std::cout << "filtration_max=" << filtration_max << std::endl; +#endif // DEBUG_TRACES + st_.set_filtration(filtration_max); + } + + 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_shapes & alpha_shape) { + // TODO: Program terminated with signal SIGABRT, Aborted - Maybe because of copy constructor + Gudhi::Simplex_tree<> st = alpha_shape.st_; + os << st << std::endl; + return os; + } +}; + +} // namespace alphashapes + +} // namespace Gudhi + +#endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ diff --git a/src/Alpha_complex/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_complex/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h new file mode 100644 index 00000000..a4e5e2fe --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h @@ -0,0 +1,192 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ +#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ +#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ + +#include <string> +#include <vector> +#include <fstream> +#include <map> + +#include "gudhi/Off_reader.h" + +namespace Gudhi { + +namespace alphashapes { + +/** + *@brief Off reader visitor with flag that can be passed to Off_reader to read a Delaunay_triangulation_complex. + */ +template<typename Complex> +class Delaunay_triangulation_off_visitor_reader { + Complex& complex_; + typedef typename Complex::Point Point; + + public: + + explicit Delaunay_triangulation_off_visitor_reader(Complex& complex) : + complex_(complex) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { +#ifdef DEBUG_TRACES + std::cout << "Delaunay_triangulation_off_visitor_reader::init - dim=" << dim << " - num_vertices=" << + num_vertices << " - num_faces=" << num_faces << " - num_edges=" << num_edges << std::endl; +#endif // DEBUG_TRACES + if (num_faces > 0) { + std::cerr << "Delaunay_triangulation_off_visitor_reader::init faces are not taken into account from OFF " << + "file for Delaunay triangulation - faces are computed." << std::endl; + } + if (num_edges > 0) { + std::cerr << "Delaunay_triangulation_off_visitor_reader::init edges are not taken into account from OFF " << + "file for Delaunay triangulation - edges are computed." << std::endl; + } + //complex_.set_current_dimension(dim); + } + + void point(const std::vector<double>& point) { +#ifdef DEBUG_TRACES + std::cout << "Delaunay_triangulation_off_visitor_reader::point "; + for (auto coordinate : point) { + std::cout << coordinate << " | "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + complex_.insert(Point(point.size(), point.begin(), point.end())); + } + + void maximal_face(const std::vector<int>& face) { + // For Delaunay Triangulation, only points are read +#ifdef DEBUG_TRACES + std::cout << "Delaunay_triangulation_off_visitor_reader::face "; + for (auto vertex : face) { + std::cout << vertex << " | "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + } + + void done() { +#ifdef DEBUG_TRACES + std::cout << "Delaunay_triangulation_off_visitor_reader::done" << std::endl; +#endif // DEBUG_TRACES + } +}; + +/** + *@brief Class that allows to load a Delaunay_triangulation_complex from an off file. + */ +template<typename Complex> +class Delaunay_triangulation_off_reader { + public: + + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex) : valid_(false) { + std::ifstream stream(name_file); + if (stream.is_open()) { + Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } else { + std::cerr << "Delaunay_triangulation_off_reader::Delaunay_triangulation_off_reader could not open file " << + name_file << std::endl; + } + + } + + /** + * return true if reading did not meet problems. + */ + bool is_valid() const { + return valid_; + } + + private: + bool valid_; +}; + +template<typename Complex> +class Delaunay_triangulation_off_writer { + public: + typedef typename Complex::Point Point; + + /** + * name_file : file where the off will be written + * save_complex : complex that be outputted in the file + * for now only save triangles. + */ + Delaunay_triangulation_off_writer(const std::string & name_file, const Complex& save_complex) { + std::ofstream stream(name_file); + if (stream.is_open()) { + if (save_complex.current_dimension() == 3) { + // OFF header + stream << "OFF" << std::endl; + // no endl on next line - don't know why... + stream << save_complex.number_of_vertices() << " " << save_complex.number_of_finite_full_cells() << " 0"; + } else { + // nOFF header + stream << "nOFF" << std::endl; + // no endl on next line - don't know why... + stream << save_complex.current_dimension() << " " << save_complex.number_of_vertices() << " " << + save_complex.number_of_finite_full_cells() << " 0"; + + } + + // bimap to retrieve vertex handles from points and vice versa + std::map< Point, int > points_to_vh; + // Start to insert at default handle value + int vertex_handle = int(); + + // Points list + for (auto vit = save_complex.vertices_begin(); vit != save_complex.vertices_end(); ++vit) { + for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { + stream << *Coord << " "; + } + stream << std::endl; + points_to_vh[vit->point()] = vertex_handle; + vertex_handle++; + } + + for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) { + std::vector<int> vertexVector; + stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + stream << points_to_vh[(*vit)->point()] << " "; + } + stream << std::endl; + } + stream.close(); + } else { + std::cerr << "Delaunay_triangulation_off_writer::Delaunay_triangulation_off_writer could not open file " << + name_file << std::endl; + } + } +}; + +} // namespace alphashapes + +} // namespace Gudhi + +#endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ diff --git a/src/Alpha_complex/test/Alpha_shapes_unit_test.cpp b/src/Alpha_complex/test/Alpha_shapes_unit_test.cpp new file mode 100644 index 00000000..d5db3bfa --- /dev/null +++ b/src/Alpha_complex/test/Alpha_shapes_unit_test.cpp @@ -0,0 +1,126 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#define BOOST_TEST_MODULE alpha_shapes +#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/Alpha_shapes/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_shapes.h" + +// to construct a simplex_tree from Delaunay_triangulation +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Simplex_tree.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// 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 instanciation of the +// TriangulationDataStructure template parameter + +BOOST_AUTO_TEST_CASE( OFF_file ) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("S4_100.off"); + std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl; + + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, 4); + + const int DIMENSION = 4; + std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.dimension() == DIMENSION); + + const double FILTRATION = 0.0; + std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.filtration() == FILTRATION); + + const int NUMBER_OF_VERTICES = 100; + std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.num_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_SIMPLICES = 6779; + std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + +} + +BOOST_AUTO_TEST_CASE( Delaunay_triangulation ) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a Delauny triangulation + // + // ---------------------------------------------------------------------------- + T dt(8); + std::string off_file_name("S8_10.off"); + std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl; + + Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt); + std::cout << "off_reader.is_valid()=" << off_reader.is_valid() << std::endl; + BOOST_CHECK(off_reader.is_valid()); + + const int NUMBER_OF_VERTICES = 10; + std::cout << "dt.number_of_vertices()=" << dt.number_of_vertices() << std::endl; + BOOST_CHECK(dt.number_of_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_FULL_CELLS = 30; + std::cout << "dt.number_of_full_cells()=" << dt.number_of_full_cells() << std::endl; + BOOST_CHECK(dt.number_of_full_cells() == NUMBER_OF_FULL_CELLS); + + const int NUMBER_OF_FINITE_FULL_CELLS = 6; + std::cout << "dt.number_of_finite_full_cells()=" << dt.number_of_finite_full_cells() << std::endl; + BOOST_CHECK(dt.number_of_finite_full_cells() == NUMBER_OF_FINITE_FULL_CELLS); + + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt); + + const int DIMENSION = 8; + std::cout << "alpha_shapes_from_dt.dimension()=" << alpha_shapes_from_dt.dimension() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.dimension() == DIMENSION); + + const double FILTRATION = 0.0; + std::cout << "alpha_shapes_from_dt.filtration()=" << alpha_shapes_from_dt.filtration() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.filtration() == FILTRATION); + + std::cout << "alpha_shapes_from_dt.num_vertices()=" << alpha_shapes_from_dt.num_vertices() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.num_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_SIMPLICES = 997; + std::cout << "alpha_shapes_from_dt.num_simplices()=" << alpha_shapes_from_dt.num_simplices() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.num_simplices() == NUMBER_OF_SIMPLICES); +} + diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt new file mode 100644 index 00000000..3cf97b71 --- /dev/null +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIAlphaShapesTest) + +# need CGAL 4.7 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.7.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + add_definitions(-DDEBUG_TRACES) + add_executable ( AlphaShapesUnitTest Alpha_shapes_unit_test.cpp ) + target_link_libraries(AlphaShapesUnitTest ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_test(AlphaShapesUnitTest ${CMAKE_CURRENT_BINARY_DIR}/AlphaShapesUnitTest) + + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Alpha_shapes/include/gudhi") diff --git a/src/Alpha_complex/test/README b/src/Alpha_complex/test/README new file mode 100644 index 00000000..244a2b84 --- /dev/null +++ b/src/Alpha_complex/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./AlphaShapesUnitTest --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Alpha_complex/test/S4_100.off b/src/Alpha_complex/test/S4_100.off new file mode 100644 index 00000000..0a5dc58c --- /dev/null +++ b/src/Alpha_complex/test/S4_100.off @@ -0,0 +1,102 @@ +OFF +100 0 0 +0.562921 -0.735261 -0.256472 0.277007 +-0.803733 -0.0527915 -0.315125 0.501918 +-0.24946 -0.354982 -0.410773 -0.801887 +0.916381 -0.0512295 0.371049 0.141223 +0.182222 0.940836 -0.171362 -0.228599 +-0.787145 -0.129213 0.568102 -0.202402 +0.0187866 -0.22093 -0.832882 -0.507095 +-0.4702 -0.533814 0.353776 0.607286 +-0.159798 -0.504771 0.263586 0.806346 +0.295546 0.162541 0.452931 0.825279 +0.242043 -0.107437 -0.913612 -0.308521 +0.875759 -0.113035 -0.469189 -0.0114505 +0.547877 -0.762247 -0.256972 -0.229729 +-0.172302 0.521057 0.412013 -0.727363 +-0.724729 -0.0574074 -0.0290602 0.686023 +0.700434 -0.102636 0.687285 0.162779 +-0.681386 0.0946893 0.610047 0.393178 +-0.847553 -0.357132 0.383743 0.0827718 +0.72297 -0.161631 -0.608517 0.284424 +0.757394 0.141549 0.196065 -0.606528 +0.78094 0.00901997 0.434536 0.448586 +0.14166 -0.619339 0.614589 -0.467582 +0.473105 -0.537832 -0.0103746 0.697711 +-0.208004 0.536218 0.818027 0.00605288 +0.743694 -0.628926 0.188072 0.126488 +-0.462228 -0.278147 0.35514 -0.76345 +-0.17361 0.249211 0.758567 -0.57648 +0.416958 -0.254924 -0.576373 -0.654946 +-0.590751 -0.286089 -0.424896 0.623402 +-0.639538 -0.739693 -0.203745 0.0482932 +0.0731787 0.132121 0.864022 0.480266 +-0.149644 -0.164724 -0.249746 -0.94239 +-0.348592 0.0120379 -0.928656 -0.126239 +0.395328 -0.54513 0.149976 -0.723917 +-0.974164 0.14707 0.157191 0.068302 +-0.166425 0.119943 -0.627304 -0.751269 +0.031947 -0.358518 -0.708301 -0.607251 +0.93781 -0.155368 -0.30951 0.0240237 +0.276094 -0.753155 -0.597088 -0.00387459 +-0.642876 -0.200551 -0.263517 -0.690687 +0.178711 0.604987 -0.262989 0.729993 +-0.520347 0.497922 -0.676144 0.155374 +-0.703999 0.500219 -0.484381 0.139789 +-0.131013 0.835735 0.506779 0.166004 +-0.536116 -0.566557 0.229226 0.582279 +-0.334105 0.158252 0.926091 0.0754059 +-0.0362677 0.296076 0.897108 0.325915 +-0.57486 0.798575 0.15324 -0.0912754 +0.498602 0.0186805 0.72824 0.469801 +-0.960329 0.0473356 0.261005 -0.0860505 +0.899134 -0.381392 -0.214508 0.00921711 +0.570576 0.567224 0.393019 -0.445237 +-0.761763 -0.614589 -0.0546476 -0.197513 +0.188584 0.289531 0.174031 0.922129 +-0.458506 -0.583876 0.639297 -0.2004 +0.785343 -0.21571 0.0794082 -0.574804 +0.0819036 0.65961 -0.247426 0.704973 +0.573125 0.49706 0.373026 0.534145 +-0.513286 -0.626226 0.208535 -0.548536 +0.460558 0.468686 0.507832 -0.55707 +0.716158 -0.488201 0.388209 -0.313164 +0.881074 0.152441 0.380128 -0.236589 +0.885793 0.0386389 0.161009 -0.433537 +-0.365162 0.298384 0.292846 0.831784 +0.364934 0.632269 -0.197205 -0.654346 +-0.31469 -0.429991 0.665304 -0.522923 +-0.734198 0.462914 -0.135691 -0.477756 +-0.422885 0.674444 -0.364143 -0.483419 +0.829218 -0.154622 -0.381147 0.378439 +0.887881 0.310479 -0.109528 0.321363 +-0.354398 -0.693974 0.456019 -0.429941 +-0.492045 -0.160008 0.044387 0.854587 +0.0595532 0.158421 0.412577 -0.895062 +-0.211441 0.491794 -0.153521 0.83058 +-0.33558 -0.504711 0.353831 -0.71236 +-0.735211 -0.197714 0.525626 0.379593 +0.465818 -0.424245 0.769469 -0.104627 +-0.641071 -0.286339 -0.704442 -0.103923 +-0.00446569 0.0249849 -0.194417 -0.980591 +-0.610081 -0.252448 0.176698 -0.729966 +-0.0859217 -0.154471 0.715027 0.676382 +0.091315 0.0723382 -0.855023 -0.505337 +0.165362 0.200983 -0.428242 -0.865373 +-0.587465 0.303019 -0.152442 0.734729 +0.454946 -0.319828 0.437063 -0.706902 +-0.384368 0.277509 0.879225 -0.0470385 +0.523335 -0.330233 -0.208592 0.757335 +0.895086 0.0448492 0.268089 -0.353466 +-0.0272491 -0.567336 -0.72254 -0.39411 +-0.0745014 -0.121818 -0.882466 0.448179 +0.382304 -0.240135 0.851109 -0.267941 +-0.418057 -0.852847 -0.3128 0.00606452 +-0.554046 0.304237 0.272381 -0.725453 +0.155115 -0.0894732 -0.245017 -0.952838 +0.114459 -0.130722 0.953669 0.245614 +0.0913002 -0.462466 0.244433 0.847374 +-0.198849 0.0785111 0.131441 -0.967997 +-0.303154 -0.686484 0.639333 0.167604 +0.521455 0.256835 -0.0584503 -0.811606 +-0.109787 0.870544 0.161523 0.451676 diff --git a/src/Alpha_complex/test/S8_10.off b/src/Alpha_complex/test/S8_10.off new file mode 100644 index 00000000..1d67e10f --- /dev/null +++ b/src/Alpha_complex/test/S8_10.off @@ -0,0 +1,12 @@ +OFF +10 0 0 +0.440036 -0.574754 -0.200485 0.216537 -0.501251 -0.0329236 -0.196529 0.313023 +-0.129367 -0.184089 -0.213021 -0.415848 0.783529 -0.0438025 0.317256 0.120749 +0.132429 0.683748 -0.124536 -0.166133 -0.540695 -0.0887576 0.390234 -0.139031 +0.0137399 -0.161581 -0.609142 -0.370872 -0.320669 -0.364053 0.24127 0.41416 +-0.115313 -0.36425 0.190208 0.581871 0.204605 0.112527 0.313562 0.571337 +0.168272 -0.0746917 -0.635156 -0.214488 0.629498 -0.0812499 -0.337255 -0.00823068 +0.369896 -0.514626 -0.173493 -0.1551 -0.127105 0.384377 0.303936 -0.536566 +-0.49013 -0.0388242 -0.0196532 0.463953 0.515962 -0.0756047 0.506276 0.119908 +-0.434258 0.060347 0.388793 0.250579 -0.653127 -0.275207 0.295714 0.0637842 +0.596172 -0.133284 -0.501793 0.234541 0.428452 0.0800735 0.110912 -0.343109 |