diff options
| -rw-r--r-- | src/Alpha_complex/doc/Intro_alpha_complex.h | 5 | ||||
| -rw-r--r-- | src/Alpha_complex/example/Alpha_complex_from_off.cpp | 6 | ||||
| -rw-r--r-- | src/Alpha_complex/example/CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 24 | ||||
| -rw-r--r-- | src/Alpha_complex/test/Alpha_complex_unit_test.cpp | 54 | ||||
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 91 | ||||
| -rw-r--r-- | src/common/include/gudhi/Delaunay_triangulation_off_io.h | 7 | ||||
| -rw-r--r-- | src/common/test/dtoffrw_alphashapedoc_result.off | 12 |
8 files changed, 166 insertions, 35 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index 2cb37578..1fb8fdee 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -20,6 +20,9 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#ifndef INTRO_ALPHA_COMPLEX_H_ +#define INTRO_ALPHA_COMPLEX_H_ + // needs namespace for Doxygen to link on classes namespace Gudhi { // needs namespace for Doxygen to link on classes @@ -117,3 +120,5 @@ namespace alphacomplex { } // namespace alphacomplex } // namespace Gudhi + +#endif // INTRO_ALPHA_COMPLEX_H_ diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp index e140fe3d..cd6f5a4b 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp @@ -25,7 +25,7 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel; Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name, alpha_square_max_value); - + // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- @@ -35,14 +35,14 @@ int main(int argc, char **argv) { 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()) { - if (alpha_complex_from_file.filtration(f_simplex) <= alpha_complex_from_file.filtration()) { + //if (alpha_complex_from_file.filtration(f_simplex) <= alpha_complex_from_file.filtration()) { std::cout << " ( "; for (auto vertex : alpha_complex_from_file.simplex_vertex_range(f_simplex)) { std::cout << vertex << " "; } std::cout << ") -> " << "[" << alpha_complex_from_file.filtration(f_simplex) << "] "; std::cout << std::endl; - } + //} } return 0; } diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index 10b87f04..24f3a9dc 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -1,6 +1,8 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIAlphaShapesExample) +add_executable ( flat flat.cpp ) + # need CGAL 4.7 # cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. if(CGAL_FOUND) diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 562b80c3..10b290b5 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -105,6 +105,7 @@ class Alpha_complex : public Simplex_tree<> { * the Alpha_complex. * * @param[in] off_file_name OFF file [path and] name. + * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$. */ Alpha_complex(const std::string& off_file_name, Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) @@ -115,25 +116,24 @@ class Alpha_complex : public Simplex_tree<> { exit(-1); // ----- >> } triangulation_ = off_reader.get_complex(); - set_filtration(max_alpha_square); - init(); + init(max_alpha_square); } /** \brief Alpha_complex constructor from a Delaunay triangulation. * * @param[in] triangulation_ptr Pointer on a Delaunay triangulation. + * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$. */ Alpha_complex(Delaunay_triangulation* triangulation_ptr, Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) : triangulation_(triangulation_ptr) { - set_filtration(max_alpha_square); - init(); + init(max_alpha_square); } /** \brief Alpha_complex constructor from a list of points. * - * @param[in] dimension Dimension of points to be inserted. * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d + * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$. * * The type InputPointRange must be a range for which std::begin and * std::end return input iterators on a Kernel::Point_d. @@ -155,8 +155,7 @@ class Alpha_complex : public Simplex_tree<> { std::cerr << "Alpha_complex - insertion failed " << inserted << " != " << (last -first) << std::endl; exit(-1); // ----- >> } - set_filtration(max_alpha_square); - init(); + init(max_alpha_square); } /** \brief Alpha_complex destructor from a Delaunay triangulation. @@ -180,12 +179,14 @@ class Alpha_complex : public Simplex_tree<> { private: /** \brief Initialize the Alpha_complex from the Delaunay triangulation. * + * @param[in] max_alpha_square maximum for alpha square value. + * * @warning Delaunay triangulation must be already constructed with at least one vertex and dimension must be more * than 0. * * Initialization can be launched once. */ - void init() { + void init(Filtration_value max_alpha_square) { if (triangulation_ == nullptr) { std::cerr << "Alpha_complex init - Cannot init from a NULL triangulation" << std::endl; return; // ----- >> @@ -287,6 +288,13 @@ class Alpha_complex : public Simplex_tree<> { } } // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- } template<typename Simplex_handle> diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index f64a8ea9..2912019d 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -88,20 +88,9 @@ BOOST_AUTO_TEST_CASE(ALPHA_DOC_OFF_file_filtered) { 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 = 25; + const int NUMBER_OF_SIMPLICES = 23; 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); - - int num_filtered_simplices = 0; - for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { - if (alpha_complex_from_file.filtration(f_simplex) <= alpha_complex_from_file.filtration()) { - num_filtered_simplices++; - } - } - const int NUMBER_OF_FILTERED_SIMPLICES = 23; - std::cout << "num_filtered_simplices=" << num_filtered_simplices << std::endl; - BOOST_CHECK(num_filtered_simplices == NUMBER_OF_FILTERED_SIMPLICES); - } bool are_almost_the_same(float a, float b) { @@ -140,8 +129,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points // ---------------------------------------------------------------------------- - double max_alpha_square_value = 1e10; - Gudhi::alphacomplex::Alpha_complex<Kernel_s> alpha_complex_from_points(points, max_alpha_square_value); + Gudhi::alphacomplex::Alpha_complex<Kernel_s> alpha_complex_from_points(points); std::cout << "========== Alpha_complex_from_points ==========" << std::endl; @@ -210,4 +198,42 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { BOOST_CHECK_THROW (alpha_complex_from_points.get_point(4), std::out_of_range); BOOST_CHECK_THROW (alpha_complex_from_points.get_point(-1), std::out_of_range); BOOST_CHECK_THROW (alpha_complex_from_points.get_point(1234), std::out_of_range); + + // Test after prune_above_filtration + alpha_complex_from_points.prune_above_filtration(0.6); + // Another way to check num_simplices + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + num_simplices = 0; + for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) { + num_simplices++; + 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; + } + BOOST_CHECK(num_simplices == 10); + std::cout << "alpha_complex_from_points.num_simplices()=" << alpha_complex_from_points.num_simplices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_simplices() == 10); + + 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_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == 4); + + for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) { + switch (alpha_complex_from_points.dimension(f_simplex)) { + case 0: + BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 0.0)); + break; + case 1: + BOOST_CHECK(are_almost_the_same(alpha_complex_from_points.filtration(f_simplex), 1.0/2.0)); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } + } diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 35d839e2..8c1beaef 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -39,6 +39,7 @@ #include <utility> #include <vector> #include <functional> // for greater<> +#include <limits> // for numeric_limits infinity namespace Gudhi { /** \defgroup simplex_tree Filtered Complexes @@ -1098,6 +1099,96 @@ class Simplex_tree { os << filtration(sh) << " \n"; } } + + public: + /** \brief Browse the simplex tree to ensure the filtration is not decreasing. + * @return The filtration modification information in order to trigger initialize_filtration. + * \warning initialize_filtration is launched again in case of filtration modification change. + */ + bool make_filtration_non_decreasing() { + bool modified = false; + for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) { + if (has_children(sh)) { + modified = modified || rec_make_filtration_non_decreasing(sh->second.children(), sh->second.filtration()); + } + } + if (modified) { + initialize_filtration(); + } + return modified; + } + + private: + /** \brief Recursively Browse the simplex tree to ensure the filtration is not decreasing. + * @param[in] sib Siblings to be parsed. + * @param[in] upper_filtration Upper level filtration value in the simplex tree. + * @return The filtration modification information in order to trigger initialize_filtration. + */ + bool rec_make_filtration_non_decreasing(Siblings * sib, Filtration_value upper_filtration) { + bool modified = false; + for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) { + if (sh->second.filtration() < upper_filtration) { + // Store the filtration modification information + modified = true; + std::cout << "modified" << std::endl; + sh->second.assign_filtration(upper_filtration); + } + if (has_children(sh)) { + modified = modified || rec_make_filtration_non_decreasing(sh->second.children(), sh->second.filtration()); + } + } + // Make the modified information to be traced by upper call + return modified; + } + + public: + /** \brief Prune above filtration value given as parameter. + * @param[in] filtration Maximum threshold value. + * \warning threshold_ is set from filtration given as parameter. + * \warning The filtration must be valid. If the filtration has not been initialized yet, the method initializes it + * (i.e. order the simplices). If the complex has changed since the last time the filtration was initialized, please + * call `initialize_filtration()` to recompute it. + */ + void prune_above_filtration(Filtration_value filtration) { + threshold_ = filtration; + if (filtration != std::numeric_limits<Filtration_value>::infinity()) { + // Initialize filtration_vect_ if required + if (filtration_vect_.empty()) { + initialize_filtration(); + } + + // Loop in reverse mode until threshold is reached + auto f_simplex = filtration_vect_.rbegin(); + for (; f_simplex != filtration_vect_.rend() && ((*f_simplex)->second.filtration() > threshold_); f_simplex++) { + remove_maximal_simplex(*f_simplex); + } + // Do not forget to update filtration_vect_ - resize is enough + std::size_t new_size = filtration_vect_.size() - (f_simplex - filtration_vect_.rbegin()); + filtration_vect_.resize(new_size); + } + } + + private: + /** \brief Remove a maximal simplex. + * @param[in] sh Simplex handle on the maximal simplex to remove. + * \warning Exception std::invalid_argument is thrown in sh has children. + */ + void remove_maximal_simplex(Simplex_handle sh) { + // Guarantee the simplex is maximal + if (has_children(sh)) { + throw std::invalid_argument ("Simplex_tree::remove_maximal_simplex - argument is not a maximal simplex"); + } + // Simplex is a leaf, it means the child is the Siblings owning the leaf. + Siblings* child = sh->second.children(); + if (child->size() > 1) { + // Not alone, just remove it from members + child->members().erase(sh->first); + } else { + // Sibling is emptied : must be deleted, and its parent must point on his own Sibling + child->oncles()->members().at(child->parent()).assign_children(child->oncles()); + delete child; + } + } private: Vertex_handle null_vertex_; diff --git a/src/common/include/gudhi/Delaunay_triangulation_off_io.h b/src/common/include/gudhi/Delaunay_triangulation_off_io.h index 47066a94..4d26bb71 100644 --- a/src/common/include/gudhi/Delaunay_triangulation_off_io.h +++ b/src/common/include/gudhi/Delaunay_triangulation_off_io.h @@ -19,8 +19,8 @@ * 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_ +#ifndef DELAUNAY_TRIANGULATION_OFF_IO_H_ +#define DELAUNAY_TRIANGULATION_OFF_IO_H_ #include <string> #include <vector> @@ -256,7 +256,6 @@ class Delaunay_triangulation_off_writer { // no endl on next line - don't know why... stream << complex_ptr->current_dimension() << " " << complex_ptr->number_of_vertices() << " " << complex_ptr->number_of_finite_full_cells() << " 0"; - } // bimap to retrieve vertex handles from points and vice versa @@ -305,4 +304,4 @@ class Delaunay_triangulation_off_writer { } // namespace Gudhi -#endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ +#endif // DELAUNAY_TRIANGULATION_OFF_IO_H_ diff --git a/src/common/test/dtoffrw_alphashapedoc_result.off b/src/common/test/dtoffrw_alphashapedoc_result.off index 13c255c6..03b7ca75 100644 --- a/src/common/test/dtoffrw_alphashapedoc_result.off +++ b/src/common/test/dtoffrw_alphashapedoc_result.off @@ -7,9 +7,9 @@ nOFF 0 14 2 19 9 17 -3 1 2 3 -3 4 3 2 -3 5 1 3 -3 5 3 7 -3 7 3 4 -3 6 5 7 +3 0 1 2 +3 3 2 1 +3 4 0 2 +3 4 2 6 +3 6 2 3 +3 5 4 6 |