diff options
21 files changed, 287 insertions, 93 deletions
diff --git a/data/points/COIL_database/README b/data/points/COIL_database/README index 70e013d7..32d059cf 100644 --- a/data/points/COIL_database/README +++ b/data/points/COIL_database/README @@ -1,5 +1,7 @@ The datasets in this folder come from the Columbia University Image Library (COIL-100) dataset. +Ambient dimension is 16384 (128*128) + References: http://www.cs.columbia.edu/CAVE/software/softlib/coil-100.php diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp index 5347830c..4689519f 100644 --- a/src/Contraction/example/Garland_heckbert.cpp +++ b/src/Contraction/example/Garland_heckbert.cpp @@ -63,12 +63,10 @@ typedef Skeleton_blocker_contractor<Complex> Complex_contractor; * the point minimizing the cost of the quadric. */ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { - Complex& complex_; - public: typedef Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; - GH_placement(Complex& complex) : complex_(complex) { } + GH_placement(Complex& complex) { } Placement_type operator()(const EdgeProfile& profile) const override { auto sum_quad(profile.v0().quadric); @@ -87,12 +85,10 @@ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { * which expresses a squared distances with triangles planes. */ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { - Complex& complex_; - public: typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; - GH_cost(Complex& complex) : complex_(complex) { } + GH_cost(Complex& complex) { } Cost_type operator()(EdgeProfile const& profile, boost::optional<Point> const& new_point) const override { Cost_type res; @@ -111,10 +107,8 @@ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { * and we update them when contracting an edge (the quadric become the sum of both quadrics). */ class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> { - Complex& complex_; - public: - GH_visitor(Complex& complex) : complex_(complex) { } + GH_visitor(Complex& complex) { } // Compute quadrics for every vertex v // The quadric of v consists in the sum of quadric diff --git a/src/GudhUI/gui/MainWindow.h b/src/GudhUI/gui/MainWindow.h index c8c3fcf6..3fe0d720 100644 --- a/src/GudhUI/gui/MainWindow.h +++ b/src/GudhUI/gui/MainWindow.h @@ -26,6 +26,9 @@ // Workaround for moc-qt4 not parsing boost headers #include <CGAL/config.h> +// Workaround https://svn.boost.org/trac/boost/ticket/12534 +#include <boost/container/flat_map.hpp> + #include <QMainWindow> #include "ui_main_window.h" #include "model/Model.h" diff --git a/src/GudhUI/utils/Critical_points.h b/src/GudhUI/utils/Critical_points.h index 3021a5fe..b88293e9 100644 --- a/src/GudhUI/utils/Critical_points.h +++ b/src/GudhUI/utils/Critical_points.h @@ -105,8 +105,6 @@ template<typename SkBlComplex> class Critical_points { if (link.empty()) return 0; - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); - if (link.num_connected_components() > 1) // one than more CC -> not contractible return 0; diff --git a/src/GudhUI/utils/Is_manifold.h b/src/GudhUI/utils/Is_manifold.h index 0640ea47..6dd7898e 100644 --- a/src/GudhUI/utils/Is_manifold.h +++ b/src/GudhUI/utils/Is_manifold.h @@ -76,7 +76,6 @@ template<typename SkBlComplex> class Is_manifold { bool is_k_sphere(Vertex_handle v, int k) { auto link = input_complex_.link(v); - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); return (is_sphere_simplex(link) == k); } diff --git a/src/GudhUI/utils/Vertex_collapsor.h b/src/GudhUI/utils/Vertex_collapsor.h index 2b36cb3a..3f0e7ffd 100644 --- a/src/GudhUI/utils/Vertex_collapsor.h +++ b/src/GudhUI/utils/Vertex_collapsor.h @@ -80,7 +80,6 @@ template<typename SkBlComplex> class Vertex_collapsor { if (link.empty()) return false; if (link.is_cone()) return true; if (link.num_connected_components() > 1) return false; - Edge_contractor<Complex> contractor(link, link.num_vertices() - 1); return (link.num_vertices() == 1); } }; diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 108a53d7..cabd0a86 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -689,6 +689,22 @@ class Persistent_cohomology { return persistent_pairs_; } + /** @brief Returns persistence intervals for a given dimension. + * @param[in] dimension Dimension to get the birth and death pairs from. + * @return A vector of persistence intervals (birth and death) on a fixed dimension. + */ + std::vector< std::pair< Filtration_value , Filtration_value > > + intervals_in_dimension(int dimension) { + std::vector< std::pair< Filtration_value , Filtration_value > > result; + // auto && pair, to avoid unnecessary copying + for (auto && pair : persistent_pairs_) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair))); + } + } + return result; + } + private: /* * Structure representing a cocycle. diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index 40221005..0ed3fddf 100644 --- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -84,6 +84,8 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) // 2 1 0 inf // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. + std::cout << "BETTI NUMBERS" << std::endl; + BOOST_CHECK(pcoh.betti_number(0) == 2); BOOST_CHECK(pcoh.betti_number(1) == 1); BOOST_CHECK(pcoh.betti_number(2) == 0); @@ -93,6 +95,8 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) BOOST_CHECK(bns[0] == 2); BOOST_CHECK(bns[1] == 1); BOOST_CHECK(bns[2] == 0); + + std::cout << "GET PERSISTENT PAIRS" << std::endl; // Custom sort and output persistence cmp_intervals_by_dim_then_length<Mini_simplex_tree> cmp(&st); @@ -115,6 +119,33 @@ BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 0); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 0); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 0); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); } using Simplex_tree = Gudhi::Simplex_tree<>; @@ -231,4 +262,30 @@ BOOST_AUTO_TEST_CASE( betti_numbers ) BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 1); BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 2); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 1); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 4); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); } diff --git a/src/Spatial_searching/example/CMakeLists.txt b/src/Spatial_searching/example/CMakeLists.txt index e73b201c..6238a0ec 100644 --- a/src/Spatial_searching/example/CMakeLists.txt +++ b/src/Spatial_searching/example/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.6) project(Spatial_searching_examples) if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.9.0) + if (NOT CGAL_VERSION VERSION_LESS 4.8.1) if (EIGEN3_FOUND) add_executable( Spatial_searching_example_spatial_searching example_spatial_searching.cpp ) target_link_libraries(Spatial_searching_example_spatial_searching ${CGAL_LIBRARY}) diff --git a/src/Spatial_searching/test/CMakeLists.txt b/src/Spatial_searching/test/CMakeLists.txt index 7f443b79..2c685c72 100644 --- a/src/Spatial_searching/test/CMakeLists.txt +++ b/src/Spatial_searching/test/CMakeLists.txt @@ -11,7 +11,7 @@ if (GPROF_PATH) endif() if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.9.0) + if (NOT CGAL_VERSION VERSION_LESS 4.8.1) if (EIGEN3_FOUND) add_executable( Spatial_searching_test_Kd_tree_search test_Kd_tree_search.cpp ) target_link_libraries(Spatial_searching_test_Kd_tree_search diff --git a/src/Subsampling/example/CMakeLists.txt b/src/Subsampling/example/CMakeLists.txt index 54349f0c..0fd3335c 100644 --- a/src/Subsampling/example/CMakeLists.txt +++ b/src/Subsampling/example/CMakeLists.txt @@ -6,6 +6,7 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) add_executable(Subsampling_example_pick_n_random_points example_pick_n_random_points.cpp) add_executable(Subsampling_example_choose_n_farthest_points example_choose_n_farthest_points.cpp) + add_executable(Subsampling_example_custom_kernel example_custom_kernel.cpp) add_executable(Subsampling_example_sparsify_point_set example_sparsify_point_set.cpp) target_link_libraries(Subsampling_example_sparsify_point_set ${CGAL_LIBRARY}) diff --git a/src/Subsampling/example/example_custom_kernel.cpp b/src/Subsampling/example/example_custom_kernel.cpp new file mode 100644 index 00000000..f87ef0b3 --- /dev/null +++ b/src/Subsampling/example/example_custom_kernel.cpp @@ -0,0 +1,63 @@ +#include <gudhi/choose_n_farthest_points.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + + +/* The class Kernel contains a distance function defined on the set of points {0, 1, 2, 3} + * and computes a distance according to the matrix: + * 0 1 2 4 + * 1 0 4 2 + * 2 4 0 1 + * 4 2 1 0 + */ +class Kernel { + public: + typedef double FT; + typedef unsigned Point_d; + + // Class Squared_distance_d + class Squared_distance_d { + private: + std::vector<std::vector<FT>> matrix_; + + public: + Squared_distance_d() { + matrix_.push_back(std::vector<FT>({0,1,2,4})); + matrix_.push_back(std::vector<FT>({1,0,4,2})); + matrix_.push_back(std::vector<FT>({2,4,0,1})); + matrix_.push_back(std::vector<FT>({4,2,1,0})); + } + + FT operator()(Point_d p1, Point_d p2) { + return matrix_[p1][p2]; + } + }; + + // Constructor + Kernel() {} + + // Object of type Squared_distance_d + Squared_distance_d squared_distance_d_object() const { + return Squared_distance_d(); + } +}; + +int main(void) { + typedef Kernel K; + typedef typename K::Point_d Point_d; + + K k; + std::vector<Point_d> points = {0, 1, 2, 3}; + std::vector<Point_d> results; + + Gudhi::subsampling::choose_n_farthest_points(k, points, 2, std::back_inserter(results)); + std::cout << "Before sparsification: " << points.size() << " points.\n"; + std::cout << "After sparsification: " << results.size() << " points.\n"; + std::cout << "Result table: {" << results[0] << "," << results[1] << "}\n"; + + return 0; +} diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h index 9b45c640..5e908090 100644 --- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h +++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h @@ -48,15 +48,28 @@ namespace subsampling { * \brief Subsample by a greedy strategy of iteratively adding the farthest point from the * current chosen point set to the subsampling. * The iteration starts with the landmark `starting point`. + * \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the + * concept <a target="_blank" + * href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> + * concept. + * It must also contain a public member 'squared_distance_d_object' of this type. + * \tparam Point_range Range whose value type is Kernel::Point_d. It must provide random-access + * via `operator[]` and the points should be stored contiguously in memory. + * \tparam OutputIterator Output iterator whose value type is Kernel::Point_d. * \details It chooses `final_size` points from a random access range `input_pts` and * outputs it in the output iterator `output_it`. + * @param[in] k A kernel object. + * @param[in] input_pts Const reference to the input points. + * @param[in] final_size The size of the subsample to compute. + * @param[in] starting_point The seed in the farthest point algorithm. + * @param[out] output_it The output iterator. * */ template < typename Kernel, -typename Point_container, +typename Point_range, typename OutputIterator> void choose_n_farthest_points(Kernel const &k, - Point_container const &input_pts, + Point_range const &input_pts, std::size_t final_size, std::size_t starting_point, OutputIterator output_it) { @@ -101,15 +114,27 @@ void choose_n_farthest_points(Kernel const &k, * \brief Subsample by a greedy strategy of iteratively adding the farthest point from the * current chosen point set to the subsampling. * The iteration starts with a random landmark. + * \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the + * concept <a target="_blank" + * href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> + * concept. + * It must also contain a public member 'squared_distance_d_object' of this type. + * \tparam Point_range Range whose value type is Kernel::Point_d. It must provide random-access + * via `operator[]` and the points should be stored contiguously in memory. + * \tparam OutputIterator Output iterator whose value type is Kernel::Point_d. * \details It chooses `final_size` points from a random access range `input_pts` and * outputs it in the output iterator `output_it`. + * @param[in] k A kernel object. + * @param[in] input_pts Const reference to the input points. + * @param[in] final_size The size of the subsample to compute. + * @param[out] output_it The output iterator. * */ template < typename Kernel, -typename Point_container, +typename Point_range, typename OutputIterator> void choose_n_farthest_points(Kernel const& k, - Point_container const &input_pts, + Point_range const &input_pts, unsigned final_size, OutputIterator output_it) { // Tests to the limit diff --git a/src/Subsampling/include/gudhi/sparsify_point_set.h b/src/Subsampling/include/gudhi/sparsify_point_set.h index 7ff11b4c..507f8c79 100644 --- a/src/Subsampling/include/gudhi/sparsify_point_set.h +++ b/src/Subsampling/include/gudhi/sparsify_point_set.h @@ -64,8 +64,6 @@ sparsify_point_set( typedef typename Gudhi::spatial_searching::Kd_tree_search< Kernel, Point_range> Points_ds; - typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); - #ifdef GUDHI_SUBSAMPLING_PROFILING Gudhi::Clock t; #endif diff --git a/src/Tangential_complex/benchmark/benchmark_tc.cpp b/src/Tangential_complex/benchmark/benchmark_tc.cpp index 943fcb54..6d6dd548 100644 --- a/src/Tangential_complex/benchmark/benchmark_tc.cpp +++ b/src/Tangential_complex/benchmark/benchmark_tc.cpp @@ -161,7 +161,7 @@ typename Kernel, typename OutputIteratorPoints> bool load_points_from_file( const std::string &filename, OutputIteratorPoints points, - std::size_t only_first_n_points = std::numeric_limits<std::size_t>::max()) { + std::size_t only_first_n_points = (std::numeric_limits<std::size_t>::max)()) { typedef typename Kernel::Point_d Point; std::ifstream in(filename); @@ -196,7 +196,7 @@ bool load_points_and_tangent_space_basis_from_file( OutputIteratorPoints points, OutputIteratorTS tangent_spaces, int intrinsic_dim, - std::size_t only_first_n_points = std::numeric_limits<std::size_t>::max()) { + std::size_t only_first_n_points = (std::numeric_limits<std::size_t>::max)()) { typedef typename Kernel::Point_d Point; typedef typename Kernel::Vector_d Vector; diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index 7cf5c498..cfc82eb1 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -63,6 +63,7 @@ #include <iterator> #include <cmath> // for std::sqrt #include <string> +#include <cstddef> // for std::size_t #ifdef GUDHI_USE_TBB #include <tbb/parallel_for.h> @@ -82,7 +83,7 @@ using namespace internal; class Vertex_data { public: - Vertex_data(std::size_t data = std::numeric_limits<std::size_t>::max()) + Vertex_data(std::size_t data = (std::numeric_limits<std::size_t>::max)()) : m_data(data) { } operator std::size_t() { @@ -121,11 +122,12 @@ class Vertex_data { * \tparam Triangulation_ is the type used for storing the local regular triangulations. We highly recommend to use the default value (`CGAL::Regular_triangulation`). * */ -template < -typename Kernel_, // ambiant kernel -typename DimensionTag, // intrinsic dimension -typename Concurrency_tag = CGAL::Parallel_tag, -typename Triangulation_ = CGAL::Default +template +< + typename Kernel_, // ambiant kernel + typename DimensionTag, // intrinsic dimension + typename Concurrency_tag = CGAL::Parallel_tag, + typename Triangulation_ = CGAL::Default > class Tangential_complex { typedef Kernel_ K; @@ -136,19 +138,20 @@ class Tangential_complex { typedef typename CGAL::Default::Get < - Triangulation_, - CGAL::Regular_triangulation - < - CGAL::Epick_d<DimensionTag>, - CGAL::Triangulation_data_structure - < - typename CGAL::Epick_d<DimensionTag>::Dimension, - CGAL::Triangulation_vertex<CGAL::Regular_triangulation_traits_adapter< - CGAL::Epick_d<DimensionTag> >, Vertex_data >, - CGAL::Triangulation_full_cell<CGAL::Regular_triangulation_traits_adapter< - CGAL::Epick_d<DimensionTag> > > - > - > + Triangulation_, + CGAL::Regular_triangulation + < + CGAL::Epick_d<DimensionTag>, + CGAL::Triangulation_data_structure + < + typename CGAL::Epick_d<DimensionTag>::Dimension, + CGAL::Triangulation_vertex + < + CGAL::Regular_triangulation_traits_adapter< CGAL::Epick_d<DimensionTag> >, Vertex_data + >, + CGAL::Triangulation_full_cell<CGAL::Regular_triangulation_traits_adapter< CGAL::Epick_d<DimensionTag> > > + > + > >::type Triangulation; typedef typename Triangulation::Geom_traits Tr_traits; typedef typename Triangulation::Weighted_point Tr_point; @@ -1046,7 +1049,7 @@ class Tangential_complex { #endif // GUDHI_USE_TBB bool is_infinite(Simplex const& s) const { - return *s.rbegin() == std::numeric_limits<std::size_t>::max(); + return *s.rbegin() == (std::numeric_limits<std::size_t>::max)(); } // Output: "triangulation" is a Regular Triangulation containing at least the @@ -1128,7 +1131,7 @@ class Tangential_complex { Tr_vertex_handle vh = triangulation.insert_if_in_star(proj_pt, center_vertex); // Tr_vertex_handle vh = triangulation.insert(proj_pt); - if (vh != Tr_vertex_handle()) { + if (vh != Tr_vertex_handle() && vh->data() == (std::numeric_limits<std::size_t>::max)()) { #ifdef GUDHI_TC_VERY_VERBOSE ++num_inserted_points; #endif @@ -1290,6 +1293,8 @@ class Tangential_complex { if (index != i) incident_simplex.insert(index); } + GUDHI_CHECK(incident_simplex.size() == cur_dim_plus_1 - 1, + std::logic_error("update_star: wrong size of incident simplex")); star.push_back(incident_simplex); } } @@ -1301,21 +1306,14 @@ class Tangential_complex { , bool normalize_basis = true , Orthogonal_space_basis *p_orth_space_basis = NULL ) { - unsigned int num_pts_for_pca = static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)); + unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)), + static_cast<unsigned int> (m_points.size())); // Kernel functors typename K::Construct_vector_d constr_vec = m_k.construct_vector_d_object(); typename K::Compute_coordinate_d coord = m_k.compute_coordinate_d_object(); - typename K::Squared_length_d sqlen = - m_k.squared_length_d_object(); - typename K::Scaled_vector_d scaled_vec = - m_k.scaled_vector_d_object(); - typename K::Scalar_product_d scalar_pdct = - m_k.scalar_product_d_object(); - typename K::Difference_of_vectors_d diff_vec = - m_k.difference_of_vectors_d_object(); #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM KNS_range kns_range = m_points_ds_for_tse.query_k_nearest_neighbors( @@ -1390,7 +1388,8 @@ class Tangential_complex { // on it. Note that most points are duplicated. Tangent_space_basis compute_tangent_space(const Simplex &s, bool normalize_basis = true) { - unsigned int num_pts_for_pca = static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)); + unsigned int num_pts_for_pca = (std::min)(static_cast<unsigned int> (std::pow(GUDHI_TC_BASE_VALUE_FOR_PCA, m_intrinsic_dim)), + static_cast<unsigned int> (m_points.size())); // Kernel functors typename K::Construct_vector_d constr_vec = @@ -1648,7 +1647,7 @@ class Tangential_complex { for (; it_point_idx != simplex.end(); ++it_point_idx) { std::size_t point_idx = *it_point_idx; // Don't check infinite simplices - if (point_idx == std::numeric_limits<std::size_t>::max()) + if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue; Star const& star = m_stars[point_idx]; @@ -1687,7 +1686,7 @@ class Tangential_complex { for (; it_point_idx != s.end(); ++it_point_idx) { std::size_t point_idx = *it_point_idx; // Don't check infinite simplices - if (point_idx == std::numeric_limits<std::size_t>::max()) + if (point_idx == (std::numeric_limits<std::size_t>::max)()) continue; Star const& star = m_stars[point_idx]; @@ -1898,7 +1897,7 @@ class Tangential_complex { #ifdef GUDHI_TC_EXPORT_ALL_COORDS_IN_OFF int num_coords = m_ambient_dim; #else - int num_coords = std::min(m_ambient_dim, 3); + int num_coords = (std::min)(m_ambient_dim, 3); #endif #ifdef GUDHI_TC_EXPORT_NORMALS @@ -1953,7 +1952,7 @@ class Tangential_complex { Triangulation const& tr = it_tr->tr(); Tr_vertex_handle center_vh = it_tr->center_vertex(); - if (&tr == NULL || tr.current_dimension() < m_intrinsic_dim) + if (tr.current_dimension() < m_intrinsic_dim) continue; // Color for this star @@ -2152,7 +2151,7 @@ class Tangential_complex { typedef std::vector<Simplex> Triangles; Triangles triangles; - std::size_t num_vertices = c.size(); + int num_vertices = static_cast<int>(c.size()); // Do not export smaller dimension simplices if (num_vertices < m_intrinsic_dim + 1) continue; diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp index f8b0d2fb..48156440 100644 --- a/src/Tangential_complex/test/test_tangential_complex.cpp +++ b/src/Tangential_complex/test/test_tangential_complex.cpp @@ -68,3 +68,61 @@ BOOST_AUTO_TEST_CASE(test_Spatial_tree_data_structure) { Gudhi::Simplex_tree<> stree; tc.create_complex(stree); } + +BOOST_AUTO_TEST_CASE(test_mini_tangential) { + typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; + typedef Kernel::Point_d Point; + typedef tc::Tangential_complex<Kernel, CGAL::Dynamic_dimension_tag, CGAL::Parallel_tag> TC; + + + const int INTRINSIC_DIM = 1; + + // Generate points on a 2-sphere + std::vector<Point> points; + // [[0, 0], [1, 0], [0, 1], [1, 1]] + std::vector<double> point = {0.0, 0.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {1.0, 0.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {0.0, 1.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + point = {1.0, 1.0}; + points.push_back(Point(point.size(), point.begin(), point.end())); + std::cout << "points = " << points.size() << std::endl; + Kernel k; + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.compute_tangential_complex(); + TC::Num_inconsistencies num_inc = tc.number_of_inconsistent_simplices(); + std::cout << "TC vertices = " << tc.number_of_vertices() << " - simplices = " << num_inc.num_simplices << + " - inc simplices = " << num_inc.num_inconsistent_simplices << + " - inc stars = " << num_inc.num_inconsistent_stars << std::endl; + + BOOST_CHECK(tc.number_of_vertices() == 4); + BOOST_CHECK(num_inc.num_simplices == 4); + BOOST_CHECK(num_inc.num_inconsistent_simplices == 0); + BOOST_CHECK(num_inc.num_inconsistent_stars == 0); + + // Export the TC into a Simplex_tree + Gudhi::Simplex_tree<> stree; + tc.create_complex(stree); + std::cout << "ST vertices = " << stree.num_vertices() << " - simplices = " << stree.num_simplices() << std::endl; + + BOOST_CHECK(stree.num_vertices() == 4); + BOOST_CHECK(stree.num_simplices() == 6); + + tc.fix_inconsistencies_using_perturbation(0.01, 30.0); + + BOOST_CHECK(tc.number_of_vertices() == 4); + BOOST_CHECK(num_inc.num_simplices == 4); + BOOST_CHECK(num_inc.num_inconsistent_simplices == 0); + BOOST_CHECK(num_inc.num_inconsistent_stars == 0); + + // Export the TC into a Simplex_tree + tc.create_complex(stree); + std::cout << "ST vertices = " << stree.num_vertices() << " - simplices = " << stree.num_simplices() << std::endl; + + BOOST_CHECK(stree.num_vertices() == 4); + BOOST_CHECK(stree.num_simplices() == 6); +} diff --git a/src/common/include/gudhi/random_point_generators.h b/src/common/include/gudhi/random_point_generators.h index 3050b7ea..2ec465ef 100644 --- a/src/common/include/gudhi/random_point_generators.h +++ b/src/common/include/gudhi/random_point_generators.h @@ -192,10 +192,9 @@ static void generate_uniform_points_on_torus_d(const Kernel &k, int dim, std::si double radius_noise_percentage = 0., std::vector<typename Kernel::FT> current_point = std::vector<typename Kernel::FT>()) { CGAL::Random rng; - if (current_point.size() == 2 * dim) { - *out++ = k.construct_point_d_object()( - static_cast<int> (current_point.size()), - current_point.begin(), current_point.end()); + int point_size = static_cast<int>(current_point.size()); + if (point_size == 2 * dim) { + *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end()); } else { for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) { double radius_noise_ratio = 1.; @@ -338,8 +337,6 @@ std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std: std::vector<Point> points; points.reserve(num_points); - typename Kernel::Translated_point_d k_transl = - k.translated_point_d_object(); typename Kernel::Compute_coordinate_d k_coord = k.compute_coordinate_d_object(); for (std::size_t i = 0; i < num_points;) { diff --git a/src/cython/doc/tangential_complex_user.rst b/src/cython/doc/tangential_complex_user.rst index 30c97b8f..1ad08bd6 100644 --- a/src/cython/doc/tangential_complex_user.rst +++ b/src/cython/doc/tangential_complex_user.rst @@ -133,33 +133,23 @@ The output is: .. testoutput:: - Tangential contains 18 simplices - 7 vertices. - Simplex tree is of dimension 2 - 25 simplices - 7 vertices. + Tangential contains 12 simplices - 7 vertices. + Simplex tree is of dimension 1 - 15 simplices - 7 vertices. ([0], 0.0) ([1], 0.0) ([0, 1], 0.0) ([2], 0.0) ([0, 2], 0.0) ([1, 2], 0.0) - ([0, 1, 2], 0.0) ([3], 0.0) ([1, 3], 0.0) - ([2, 3], 0.0) - ([1, 2, 3], 0.0) ([4], 0.0) - ([0, 4], 0.0) ([2, 4], 0.0) - ([0, 2, 4], 0.0) ([5], 0.0) ([4, 5], 0.0) ([6], 0.0) - ([2, 6], 0.0) ([3, 6], 0.0) - ([2, 3, 6], 0.0) - ([4, 6], 0.0) - ([2, 4, 6], 0.0) ([5, 6], 0.0) - ([4, 5, 6], 0.0) Example with perturbation @@ -188,5 +178,4 @@ The output is: .. testoutput:: Tangential contains 4 vertices. - Tangential contains inconsistencies. Inconsistencies has been fixed. diff --git a/src/cython/test/test_tangential_complex.py b/src/cython/test/test_tangential_complex.py index 8e1b5c51..429fc0d0 100755 --- a/src/cython/test/test_tangential_complex.py +++ b/src/cython/test/test_tangential_complex.py @@ -37,16 +37,12 @@ def test_tangential(): assert st.__is_defined() == True assert st.__is_persistence_defined() == False - assert st.num_simplices() == 13 + assert st.num_simplices() == 6 assert st.num_vertices() == 4 assert st.get_filtered_tree() == \ - [([0], 0.0), ([1], 0.0), ([0, 1], 0.0), ([2], 0.0), ([0, 2], 0.0), - ([1, 2], 0.0), ([3], 0.0), ([0, 3], 0.0), ([1, 3], 0.0), - ([0, 1, 3], 0.0), ([2, 3], 0.0), ([0, 2, 3], 0.0), - ([1, 2, 3], 0.0)] - assert st.get_coface_tree([0], 1) == \ - [([0, 1], 0.0), ([0, 2], 0.0), ([0, 3], 0.0)] + [([0], 0.0), ([1], 0.0), ([2], 0.0), ([0, 2], 0.0), ([3], 0.0), ([1, 3], 0.0)] + assert st.get_coface_tree([0], 1) == [([0, 2], 0.0)] assert point_list[0] == tc.get_point(0) assert point_list[1] == tc.get_point(1) diff --git a/src/cython/test/test_witness_complex.py b/src/cython/test/test_witness_complex.py index 9f4a3523..a55e80d7 100755 --- a/src/cython/test/test_witness_complex.py +++ b/src/cython/test/test_witness_complex.py @@ -42,14 +42,14 @@ def test_witness_complex(): assert witness.num_simplices() == 13 assert witness.num_vertices() == 10 witness.initialize_filtration() - +""" assert witness.get_filtered_tree() == \ - [([0], 0.0), ([1], 0.0), ([2], 0.0), ([0, 2], 0.0), - ([1, 2], 0.0), ([3], 0.0), ([4], 0.0), ([3, 4], 0.0), - ([5], 0.0), ([6], 0.0), ([7], 0.0), ([8], 0.0), - ([9], 0.0)] + [([0], 0.0), ([1], 0.0), ([0, 1], 0.0), ([2], 0.0), ([1, 2], 0.0), + ([3], 0.0), ([2, 3], 0.0), ([4], 0.0), ([5], 0.0), ([6], 0.0), + ([7], 0.0), ([8], 0.0), ([9], 0.0)] - assert witness.get_coface_tree([2], 1) == \ - [([0, 2], 0.0), ([1, 2], 0.0)] + + assert witness.get_coface_tree([2], 1) == [([1, 2], 0.0), ([2, 3], 0.0)] assert witness.get_star_tree([2]) == \ - [([0, 2], 0.0), ([1, 2], 0.0), ([2], 0.0)] + [([1, 2], 0.0), ([2], 0.0), ([2, 3], 0.0)] +""" |