diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 12 | ||||
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 71 |
2 files changed, 55 insertions, 28 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 04fb185c..d8568be0 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -219,7 +219,7 @@ class Alpha_complex : public Simplex_tree<> { } } // -------------------------------------------------------------------------------------------- - + // -------------------------------------------------------------------------------------------- // 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) { @@ -242,10 +242,6 @@ class Alpha_complex : public Simplex_tree<> { // 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()); - if (!insert_result.second) { - std::cerr << "Alpha_complex::init insert_simplex_and_subfaces failed" << std::endl; - exit(-1); // ----->> - } } // -------------------------------------------------------------------------------------------- @@ -324,8 +320,10 @@ class Alpha_complex : public Simplex_tree<> { pointVector.push_back(get_point(vertex)); } // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function + Point_d point_for_gabriel; for (auto vertex : simplex_vertex_range(f_simplex)) { - if (std::find(pointVector.begin(), pointVector.end(), get_point(vertex)) == pointVector.end()) { + point_for_gabriel = get_point(vertex); + if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) { // vertex is not found in Tau vertexForGabriel = vertex; // No need to continue loop @@ -334,7 +332,7 @@ class Alpha_complex : public Simplex_tree<> { } // is_gabriel function initialization Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); - bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), get_point(vertexForGabriel)) + bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel) != CGAL::ON_BOUNDED_SIDE; #ifdef DEBUG_TRACES std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index fe2faf90..bda7e72a 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -515,33 +515,62 @@ class Simplex_tree { * @param[in] filtration the filtration value assigned to the new N-simplex. */ template<class RandomAccessVertexRange> - std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex, + std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const 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()]); - } - // (N-1)-Simplex recursive call - returned = insert_simplex_and_subfaces(NsimplexMinusOne, filtration); + // Simplex copy + std::vector<Vertex_handle> the_simplex(Nsimplex.begin(), Nsimplex.end()); + // must be sorted in increasing order + std::sort(the_simplex.begin(), the_simplex.end()); + std::vector<std::vector<Vertex_handle>> to_be_inserted; + std::vector<std::vector<Vertex_handle>> to_be_propagated; + return rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration); + } + + private: + std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces(std::vector<Vertex_handle>& the_simplex, + std::vector<std::vector<Vertex_handle>>& to_be_inserted, + std::vector<std::vector<Vertex_handle>>& to_be_propagated, + Filtration_value filtration = 0.0) { + std::pair<Simplex_handle, bool> insert_result; + if (the_simplex.size() > 1) { + // Get and remove last vertex + Vertex_handle last_vertex= the_simplex.back(); + the_simplex.pop_back(); + // Recursive call after last vertex removal + insert_result = rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration); + + // Concatenation of to_be_inserted and to_be_propagated + to_be_inserted.insert(to_be_inserted.end(), to_be_propagated.begin(), to_be_propagated.end()); + to_be_propagated = to_be_inserted; + + // to_be_inserted treatment + for (auto& simplex_tbi : to_be_inserted) { + simplex_tbi.push_back(last_vertex); } - // N-Simplex insert - returned = insert_simplex(Nsimplex, filtration); - if (returned.second == true) { + std::vector<Vertex_handle> last_simplex(1, last_vertex); + to_be_inserted.push_back(last_simplex); + + for (auto& simplex_tbi : to_be_inserted) { + insert_result = insert_simplex(simplex_tbi, filtration); } - } else if (Nsimplex.size() == 1) { - // 1-Simplex insert - End of recursivity - returned = insert_simplex(Nsimplex, filtration); - if (returned.second == true) { + + + } else if (the_simplex.size() == 1) { + // When reaching the end of recursivity, vector of simplices shall be empty and filled on back recursive + if ((to_be_inserted.size() != 0) || (to_be_propagated.size() != 0)){ + std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty" << std::endl; + exit(-1); } + std::vector<Vertex_handle> first_simplex(1, the_simplex.back()); + to_be_inserted.push_back(first_simplex); + + insert_result = insert_simplex(first_simplex, filtration); } else { - // Nothing to insert - empty vector + std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error" << std::endl; + exit(-1); } - return returned; + + return insert_result; } public: |