summaryrefslogtreecommitdiff
path: root/src/Alpha_shapes/include
diff options
context:
space:
mode:
Diffstat (limited to 'src/Alpha_shapes/include')
-rw-r--r--src/Alpha_shapes/include/gudhi/Alpha_shapes.h200
-rw-r--r--src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h213
2 files changed, 413 insertions, 0 deletions
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
new file mode 100644
index 00000000..b8efdb4d
--- /dev/null
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
@@ -0,0 +1,200 @@
+/* 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 <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <iostream>
+#include <iterator>
+#include <vector>
+#include <string>
+
+namespace Gudhi {
+
+namespace alphashapes {
+
+/** \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;
+
+ // 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;
+
+ private:
+ /** \brief Upper bound on the simplex tree of the simplicial complex.*/
+ Gudhi::Simplex_tree<> _st;
+
+ public:
+
+ Alpha_shapes(std::string off_file_name, int dimension) {
+ Delaunay_triangulation dt(dimension);
+ Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation>
+ off_reader(off_file_name, dt, true, true);
+ 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());
+ _st.set_filtration(0.0);
+ // triangulation points list
+ for (auto vit = triangulation.finite_vertices_begin();
+ vit != triangulation.finite_vertices_end(); ++vit) {
+ typeVectorVertex vertexVector;
+ Vertex_handle vertexHdl = std::distance(triangulation.finite_vertices_begin(), vit);
+ vertexVector.push_back(vertexHdl);
+
+ // Insert each point in the simplex tree
+ _st.insert_simplex(vertexVector, 0.0);
+
+#ifdef DEBUG_TRACES
+ std::cout << "P" << vertexHdl << ":";
+ for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
+ std::cout << *Coord << " ";
+ }
+ std::cout << std::endl;
+#endif // DEBUG_TRACES
+ }
+ // triangulation finite full cells list
+ for (auto cit = triangulation.finite_full_cells_begin();
+ cit != triangulation.finite_full_cells_end(); ++cit) {
+ typeVectorVertex vertexVector;
+ for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+ // Vertex handle is distance - 1
+ Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1;
+ vertexVector.push_back(vertexHdl);
+ }
+ // Insert each point in the simplex tree
+ _st.insert_simplex_and_subfaces(vertexVector, 0.0);
+
+#ifdef DEBUG_TRACES
+ std::cout << "C" << std::distance(triangulation.finite_full_cells_begin(), cit) << ":";
+ for (auto value : vertexVector) {
+ std::cout << value << ' ';
+ }
+ std::cout << 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_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_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
new file mode 100644
index 00000000..693b393e
--- /dev/null
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
@@ -0,0 +1,213 @@
+/* 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 "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_flag_visitor_reader {
+ Complex& complex_;
+ typedef typename Complex::Point Point;
+
+ const bool load_only_points_;
+
+ public:
+ explicit Delaunay_triangulation_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) :
+ complex_(complex),
+ load_only_points_(load_only_points) { }
+
+ void init(int dim, int num_vertices, int num_faces, int num_edges) {
+#ifdef DEBUG_TRACES
+ std::cout << "init" << std::endl;
+#endif // DEBUG_TRACES
+ }
+
+ void point(const std::vector<double>& point) {
+#ifdef DEBUG_TRACES
+ std::cout << "p ";
+ 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 alpha shapes, only points are read
+ }
+
+ void done() {
+#ifdef DEBUG_TRACES
+ std::cout << "done" << std::endl;
+#endif // DEBUG_TRACES
+ }
+};
+
+/**
+ *@brief Off reader visitor 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::Vertex_handle Vertex_handle;
+ // typedef typename Complex::Simplex_handle Simplex_handle;
+ typedef typename Complex::Point Point;
+
+ const bool load_only_points_;
+ std::vector<Point> points_;
+ // std::vector<Simplex_handle> maximal_faces_;
+
+ public:
+ explicit Delaunay_triangulation_off_visitor_reader(Complex& complex, bool load_only_points = false) :
+ complex_(complex),
+ load_only_points_(load_only_points) { }
+
+ void init(int dim, int num_vertices, int num_faces, int num_edges) {
+#ifdef DEBUG_TRACES
+ std::cout << "init - " << num_vertices << std::endl;
+#endif // DEBUG_TRACES
+ // maximal_faces_.reserve(num_faces);
+ points_.reserve(num_vertices);
+ }
+
+ void point(const std::vector<double>& point) {
+#ifdef DEBUG_TRACES
+ std::cout << "p ";
+ for (auto coordinate: point) {
+ std::cout << coordinate << " | ";
+ }
+ std::cout << std::endl;
+#endif // DEBUG_TRACES
+ points_.emplace_back(Point(point.size(), point.begin(), point.end()));
+ }
+
+ void maximal_face(const std::vector<int>& face) {
+ // For alpha shapes, only points are read
+ }
+
+ void done() {
+ complex_.insert(points_.begin(), points_.end());
+#ifdef DEBUG_TRACES
+ std::cout << "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, bool read_only_points = false,
+ bool is_flag = false) : valid_(false) {
+ std::ifstream stream(name_file);
+ if (stream.is_open()) {
+ if (is_flag) {
+ // For alpha shapes, only points are read
+ Delaunay_triangulation_off_flag_visitor_reader<Complex> off_visitor(read_complex, true);
+ Off_reader off_reader(stream);
+ valid_ = off_reader.read(off_visitor);
+ } else {
+ // For alpha shapes, only points are read
+ Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex, true);
+ Off_reader off_reader(stream);
+ valid_ = off_reader.read(off_visitor);
+ }
+ }
+ }
+
+ /**
+ * 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:
+ /**
+ * 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()) {
+ // 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";
+
+ // 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;
+ }
+
+ // Finite cells list
+ for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) {
+ stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; // Dimension
+ for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+ auto vertexHdl = *vit;
+ // auto vertexHdl = std::distance(save_complex.vertices_begin(), *vit) - 1;
+ // stream << std::distance(save_complex.vertices_begin(), *(vit)) - 1 << " ";
+ }
+ stream << std::endl;
+ }
+ stream.close();
+ } else {
+ std::cerr << "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_