diff options
author | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-06-12 10:47:38 +0000 |
---|---|---|
committer | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-06-12 10:47:38 +0000 |
commit | 24ed9e6896b0aafcaeaf19e5d2970915282fa7b7 (patch) | |
tree | 2bc02bec0eac7e145d63a56fd7fd2ea0bf565db3 /src | |
parent | 844d6205eb6705935417b0ab45b0c71230bd9ed6 (diff) |
delaunay off reader/writer fix.
alpha complex algo seems ok. tests are Nok.
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/alphashapes@611 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 37cea3b33d783a77a23fa3a54a1bee106508d981
Diffstat (limited to 'src')
-rw-r--r-- | src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp | 38 | ||||
-rw-r--r-- | src/Alpha_shapes/include/gudhi/Alpha_shapes.h | 109 | ||||
-rw-r--r-- | src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h | 36 | ||||
-rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 150 | ||||
-rw-r--r-- | src/common/include/gudhi/Off_reader.h | 4 |
5 files changed, 242 insertions, 95 deletions
diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp index 6c9fcdab..fe889ec0 100644 --- a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp +++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp @@ -44,12 +44,12 @@ typedef CGAL::Delaunay_triangulation<K> T; // TriangulationDataStructure template parameter void usage(char * const progName) { - std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl; + std::cerr << "Usage: " << progName << " inputFile.off dimension outputFile.off" << std::endl; exit(-1); // ----- >> } int main(int argc, char **argv) { - if (argc != 3) { + if (argc != 4) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; usage(argv[0]); } @@ -68,37 +68,11 @@ int main(int argc, char **argv) { 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::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; - - // Points list - /*for (T::Vertex_iterator vit = dt.vertices_begin(); vit != dt.vertices_end(); ++vit) { - for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { - std::cout << *Coord << " "; - } - std::cout << std::endl; - } - std::cout << std::endl;*/ - - /*int i = 0, j = 0; - typedef T::Full_cell_iterator Full_cell_iterator; - typedef T::Facet Facet; - - for (Full_cell_iterator cit = dt.full_cells_begin(); cit != dt.full_cells_end(); ++cit) { - if (!dt.is_infinite(cit)) { - j++; - continue; - } - Facet fct(cit, cit->index(dt.infinite_vertex())); - i++; - } - std::cout << "There are " << i << " facets on the convex hull." << std::endl; - std::cout << "There are " << j << " facets not on the convex hull." << std::endl; -*/ - - std::string offOutputFile("out.off"); + std::string outFileName(argv[3]); + std::string offOutputFile(outFileName); Gudhi::alphashapes::Delaunay_triangulation_off_writer<T> off_writer(offOutputFile, dt); return 0; diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h index 841c883a..2bc8b221 100644 --- a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h +++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h @@ -32,6 +32,9 @@ #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> @@ -208,12 +211,14 @@ class Alpha_shapes { Squared_Radius squared_radius Kinit(compute_squared_radius_d_object); Is_Gabriel is_gabriel Kinit(side_of_bounded_sphere_d_object); - std::map<Kernel::Point_d, Vertex_handle> points_to_vh; + // 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[vit->point()] = vertex_handle; + points_to_vh.insert(bimap_points_vh::value_type(vit->point(), vertex_handle)); vertex_handle++; } @@ -223,10 +228,10 @@ class Alpha_shapes { typeVectorPoint pointVector; for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { #ifdef DEBUG_TRACES - std::cout << "points_to_vh=" << points_to_vh[(*vit)->point()] << std::endl; + std::cout << "points_to_vh=" << points_to_vh.left.at((*vit)->point()) << std::endl; #endif // DEBUG_TRACES // Vector of vertex construction for simplex_tree structure - vertexVector.push_back(points_to_vh[(*vit)->point()]); + vertexVector.push_back(points_to_vh.left.at((*vit)->point())); // Vector of points for alpha_shapes filtration value computation pointVector.push_back((*vit)->point()); } @@ -244,17 +249,99 @@ class Alpha_shapes { filtration_max = fmax(filtration_max, alpha_shapes_filtration); } } - - // Loop on triangulation finite full cells list - for (auto f_simplex : st_.skeleton_simplex_range(st_.dimension() - 1)) { - std::cout << "vertex = [" << st_.filtration(f_simplex) << "] "; + + // ### 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); +#ifdef DEBUG_TRACES + std::cout << "f_simplex_dim= " << f_simplex_dim << " - decr_dim= " << decr_dim << std::endl; +#endif // DEBUG_TRACES + if (decr_dim == f_simplex_dim) { + typeVectorPoint pointVector; +#ifdef DEBUG_TRACES + std::cout << "vertex "; +#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 << (int) 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); +#ifdef DEBUG_TRACES + std::cout << "From NaN to alpha_shapes_filtration=" << 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 << "Sigma "; + for (auto vertex : st_.simplex_vertex_range(f_simplex)) { + std::cout << (int) vertex << " "; + } + std::cout << " - Tau "; + for (auto vertex : st_.simplex_vertex_range(f_boundary)) { + std::cout << (int) vertex << " "; + } + std::cout << 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; + } + } + // ### 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)) + ) { + // ### 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); +#ifdef DEBUG_TRACES + std::cout << "From Boundary to alpha_shapes_filtration=" << st_.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + } + } + } + } + } + +#ifdef DEBUG_TRACES + std::cout << "The complex contains " << st_.num_simplices() << " simplices" << std::endl; + std::cout << " - dimension " << st_.dimension() << " - filtration " << st_.filtration() << std::endl; + std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for (auto f_simplex : st_.filtration_simplex_range()) { + std::cout << " " << "[" << st_.filtration(f_simplex) << "] "; for (auto vertex : st_.simplex_vertex_range(f_simplex)) { std::cout << (int) vertex << " "; } std::cout << std::endl; - if (st_.filtration(f_simplex) == filtration_unknown) - st_.assign_filtration(f_simplex, filtration_max); // TODO(VR) Compute filtration value from simplex } +#endif // DEBUG_TRACES #ifdef DEBUG_TRACES std::cout << "filtration_max=" << filtration_max << std::endl; @@ -286,7 +373,7 @@ class Alpha_shapes { return st_.filtration(); } - friend std::ostream& operator<<(std::ostream& os, const Alpha_shapes& alpha_shape) { + 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; 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 index 1413ad89..a4e5e2fe 100644 --- 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 @@ -25,7 +25,7 @@ #include <string> #include <vector> #include <fstream> -#include <iterator> // std::distance +#include <map> #include "gudhi/Off_reader.h" @@ -130,6 +130,7 @@ class Delaunay_triangulation_off_reader { template<typename Complex> class Delaunay_triangulation_off_writer { public: + typedef typename Complex::Point Point; /** * name_file : file where the off will be written @@ -153,44 +154,29 @@ class Delaunay_triangulation_off_writer { } + // 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.full_cells_begin(); cit != save_complex.full_cells_end(); ++cit) { + 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) { - // Vector of vertex construction for simplex_tree structure - // Vertex handle is distance - 1 - //int vertexHdl = std::distance(save_complex.vertices_begin(), *vit); - // infinite cell is -1 for infinite - //vertexVector.push_back(vertexHdl); - // Vector of points for alpha_shapes filtration value computation + stream << points_to_vh[(*vit)->point()] << " "; } stream << std::endl; } - - - - - - /* - // Finite cells list - for (auto cit = save_complex.full_cells_begin(); cit != save_complex.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 << "Delaunay_triangulation_off_writer::Delaunay_triangulation_off_writer could not open file " << diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index b79e3c8f..32fb2f43 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -37,7 +37,6 @@ #include <vector> namespace Gudhi { - /** \defgroup simplex_tree Filtered Complexes * * A simplicial complex \f$\mathbf{K}\f$ @@ -84,8 +83,8 @@ namespace Gudhi { * */ template<typename IndexingTag = linear_indexing_tag, - typename FiltrationValue = double, typename SimplexKey = int // must be a signed integer type - , typename VertexHandle = int // must be a signed integer type, int convertible to it +typename FiltrationValue = double, typename SimplexKey = int // must be a signed integer type +, typename VertexHandle = int // must be a signed integer type, int convertible to it // , bool ContiguousVertexHandles = true //true is Vertex_handles are exactly the set [0;n) > class Simplex_tree { @@ -245,6 +244,7 @@ class Simplex_tree { Filtration_simplex_range filtration_simplex_range() { return filtration_simplex_range(Indexing_tag()); } + /** \brief Returns a range over the vertices of a simplex. * * The order in which the vertices are visited is the decreasing order for < on Vertex_handles, @@ -316,12 +316,14 @@ class Simplex_tree { Simplex_key key(Simplex_handle sh) { return sh->second.key(); } + /** \brief Returns the simplex associated to a key. * * The filtration must be initialized. */ Simplex_handle simplex(Simplex_key key) { return filtration_vect_[key]; } + /** \brief Returns the filtration value of a simplex. * * Called on the null_simplex, returns INFINITY. */ @@ -330,12 +332,23 @@ class Simplex_tree { return sh->second.filtration(); } else { return INFINITY; - } // filtration(); } + } + } + + /** \brief Sets the filtration value of a simplex. + * + * No action if called on the null_simplex*/ + void assign_filtration(Simplex_handle sh, Filtration_value fv) { + if (sh != null_simplex()) { + sh->second.assign_filtration(fv); + } } + /** \brief Returns an upper bound of the filtration values of the simplices. */ Filtration_value filtration() { return threshold_; } + /** \brief Returns a Simplex_handle different from all Simplex_handles * associated to the simplices in the simplicial complex. * @@ -343,20 +356,24 @@ class Simplex_tree { Simplex_handle null_simplex() { return Dictionary_it(NULL); } + /** \brief Returns a key different for all keys associated to the * simplices of the simplicial complex. */ Simplex_key null_key() { return -1; } + /** \brief Returns a Vertex_handle different from all Vertex_handles associated * to the vertices of the simplicial complex. */ Vertex_handle null_vertex() { return null_vertex_; } + /** \brief Returns the number of vertices in the complex. */ size_t num_vertices() { return root_.members_.size(); } + /** \brief Returns the number of simplices in the complex. * * Does not count the empty simplex. */ @@ -376,6 +393,7 @@ class Simplex_tree { } return dim - 1; } + /** \brief Returns an upper bound on the dimension of the simplicial complex. */ int dimension() { return dimension_; @@ -423,7 +441,6 @@ class Simplex_tree { Simplex_handle find_vertex(Vertex_handle v) { return root_.members_.begin() + v; } -//{ return root_.members_.find(v); } /** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex. * @@ -479,7 +496,6 @@ class Simplex_tree { return res_insert; } - /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of * Vertex_handles, in the simplicial complex. * @@ -487,33 +503,35 @@ class Simplex_tree { * @param[in] filtration the filtration value assigned to the new N-simplex. */ template<class RandomAccessVertexRange> - void insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex, - Filtration_value filtration = 0.0) { + std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex, + Filtration_value filtration = 0.0) { + std::pair<Simplex_handle, bool> returned; if (Nsimplex.size() > 1) { for (unsigned int NIndex = 0; NIndex < Nsimplex.size(); NIndex++) { // insert N (N-1)-Simplex RandomAccessVertexRange NsimplexMinusOne; for (unsigned int NListIter = 0; NListIter < Nsimplex.size() - 1; NListIter++) { // (N-1)-Simplex creation - NsimplexMinusOne.push_back( Nsimplex[(NIndex + NListIter) % Nsimplex.size()]); + NsimplexMinusOne.push_back(Nsimplex[(NIndex + NListIter) % Nsimplex.size()]); } // (N-1)-Simplex recursive call - insert_simplex_and_subfaces(NsimplexMinusOne, filtration); + returned = insert_simplex_and_subfaces(NsimplexMinusOne, filtration); } // N-Simplex insert - std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration); + returned = insert_simplex(Nsimplex, filtration); if (returned.second == true) { num_simplices_++; } } else if (Nsimplex.size() == 1) { // 1-Simplex insert - End of recursivity - std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration); + returned = insert_simplex(Nsimplex, filtration); if (returned.second == true) { num_simplices_++; } } else { // Nothing to insert - empty vector } + return returned; } /** \brief Assign a value 'key' to the key of the simplex @@ -540,17 +558,6 @@ class Simplex_tree { return sh->second.children(); } -// void display_simplex(Simplex_handle sh) -// { -// std::cout << " " << "[" << filtration(sh) << "] "; -// for( auto vertex : simplex_vertex_range(sh) ) -// { std::cout << vertex << " "; } -// } - - // void print(Simplex_handle sh, std::ostream& os = std::cout) - // { for(auto v : simplex_vertex_range(sh)) {os << v << " ";} - // os << std::endl; } - public: /** Returns a pointer to the root nodes of the simplex tree. */ Siblings * root() { @@ -562,10 +569,12 @@ class Simplex_tree { void set_filtration(Filtration_value fil) { threshold_ = fil; } + /** Set a number of simplices for the simplicial complex. */ void set_num_simplices(const unsigned int& num_simplices) { num_simplices_ = num_simplices; } + /** Set a dimension for the simplicial complex. */ void set_dimension(int dimension) { dimension_ = dimension; @@ -623,6 +632,7 @@ class Simplex_tree { } return ((it1 == rg1.end()) && (it2 != rg2.end())); } + /** \brief StrictWeakOrdering, for the simplices, defined by the filtration. * * It corresponds to the partial order @@ -631,8 +641,7 @@ class Simplex_tree { * to be smaller. The filtration function must be monotonic. */ struct is_before_in_filtration { explicit is_before_in_filtration(Simplex_tree * st) - : st_(st) { - } + : st_(st) { } bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const { if (st_->filtration(sh1) != st_->filtration(sh2)) { @@ -708,6 +717,7 @@ class Simplex_tree { } } } + /** \brief Expands the Simplex_tree containing only its one skeleton * until dimension max_dim. * @@ -731,6 +741,7 @@ class Simplex_tree { } private: + /** \brief Recursive expansion of the simplex tree.*/ void siblings_expansion(Siblings * siblings, // must contain elements int k) { @@ -769,6 +780,7 @@ class Simplex_tree { } } } + /** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2) * and assigns the maximal possible Filtration_value to the Nodes. */ void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection, @@ -800,6 +812,7 @@ class Simplex_tree { } } } + /** Maximum over 3 values.*/ Filtration_value maximum(Filtration_value a, Filtration_value b, Filtration_value c) { @@ -824,6 +837,92 @@ class Simplex_tree { os << filtration(sh) << " \n"; } } + //---------------------------------------------------------------------------------------------- + //---------------------------------------------------------------------------------------------- + private: + + /** Recursive search of cofaces + */ + template <class RandomAccessVertexRange> + void rec_coface(RandomAccessVertexRange &vertices, Siblings *curr_sib, Dictionary *curr_res, std::vector<Dictionary>& cofaces, unsigned int length, unsigned long codimension) { + for (auto sib = curr_sib->members().begin(); sib != curr_sib->members().end() && (vertices.empty() || sib->first <= vertices[vertices.size() - 1]); ++sib) { + bool continueRecursion = (codimension == length || curr_res->size() <= codimension); // dimension of actual simplex <= codimension + if (vertices.empty()) { + if (curr_res->size() >= length && continueRecursion) + // If we reached the end of the vertices, and the simplex has more vertices than the given simplex, we found a coface + { + curr_res->emplace(sib->first, sib->second); + bool egalDim = (codimension == length || curr_res->size() == codimension); // dimension of actual simplex == codimension + if (egalDim) + cofaces.push_back(*curr_res); + if (has_children(sib)) + rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension); + curr_res->erase(curr_res->end() - 1); + } + } else if (continueRecursion) { + if (sib->first == vertices[vertices.size() - 1]) // If curr_sib matches with the top vertex + { + curr_res->emplace(sib->first, sib->second); + bool egalDim = (codimension == length || curr_res->size() == codimension); // dimension of actual simplex == codimension + if (vertices.size() == 1 && curr_res->size() > length && egalDim) + cofaces.push_back(*curr_res); + if (has_children(sib)) { // Rec call + Vertex_handle tmp = vertices[vertices.size() - 1]; + vertices.pop_back(); + rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension); + vertices.push_back(tmp); + } + curr_res->erase(curr_res->end() - 1); + } else // (sib->first < vertices[vertices.size()-1]) + { + if (has_children(sib)) { + curr_res->emplace(sib->first, sib->second); + rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension); + curr_res->erase(curr_res->end() - 1); + } + } + } + } + } + + public: + + /** \brief Compute the cofaces of a n simplex + * \param vertices List of vertices which represent the n simplex. + * \param codimension The function returns the n+codimension-simplices. If codimension = 0, return all cofaces + * \return Vector of Dictionary, empty vector if no cofaces found. + * \warning n+codimension must be lower than Simplex_tree dimension, otherwise an an empty vector is returned. + */ + + template<class RandomAccessVertexRange> + std::vector<Dictionary> coface(const RandomAccessVertexRange &vertices, int codimension) { + RandomAccessVertexRange copy = vertices; + std::vector<Dictionary> cofaces; + std::sort(copy.begin(), copy.end(), std::greater<Vertex_handle>()); // must be sorted in decreasing order + if (root_.members().empty()) { + std::cerr << "Simplex_tree::coface - empty Simplex_tree" << std::endl; + return cofaces; // ----->> + } + if (vertices.empty()) { + std::cerr << "Simplex_tree::coface - empty vertices list" << std::endl; + return cofaces; // ----->> + } + if (codimension < 0) { + std::cerr << "Simplex_tree::coface - codimension is empty" << std::endl; + return cofaces; // ----->> + } + if (codimension + vertices.size() >= (unsigned long) dimension_) { + std::cerr << "Simplex_tree::coface - codimension + vertices list size cannot be greater than Simplex_tree dimension" << std::endl; + return cofaces; // ----->> + } + std::sort(copy.begin(), copy.end(), std::greater<Vertex_handle>()); // must be sorted in decreasing order + Dictionary res; + rec_coface(copy, &root_, &res, cofaces, vertices.size(), codimension + vertices.size()); + return cofaces; + } + + //---------------------------------------------------------------------------------------------- + //---------------------------------------------------------------------------------------------- private: Vertex_handle null_vertex_; @@ -851,6 +950,7 @@ std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) { } return os; } + template<typename T1, typename T2, typename T3> std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) { // assert(st.num_simplices() == 0); diff --git a/src/common/include/gudhi/Off_reader.h b/src/common/include/gudhi/Off_reader.h index e29218d8..618d1b4d 100644 --- a/src/common/include/gudhi/Off_reader.h +++ b/src/common/include/gudhi/Off_reader.h @@ -7,7 +7,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (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 @@ -121,7 +121,7 @@ private: if(!goto_next_uncomment_line(line)) return false; std::istringstream iss(line); - if(is_off_file){ + if((is_off_file) && (!is_noff_file)) { off_info_.dim = 3; if(!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)){ std::cerr << "incorrect number of vertices/faces/edges\n"; |