diff options
Diffstat (limited to 'src/GudhUI/utils')
-rw-r--r-- | src/GudhUI/utils/Critical_points.h | 203 | ||||
-rw-r--r-- | src/GudhUI/utils/Edge_collapsor.h | 156 | ||||
-rw-r--r-- | src/GudhUI/utils/Edge_contractor.h | 156 | ||||
-rw-r--r-- | src/GudhUI/utils/Furthest_point_epsilon_net.h | 244 | ||||
-rw-r--r-- | src/GudhUI/utils/Is_manifold.h | 117 | ||||
-rw-r--r-- | src/GudhUI/utils/K_nearest_builder.h | 143 | ||||
-rw-r--r-- | src/GudhUI/utils/Lloyd_builder.h | 133 | ||||
-rw-r--r-- | src/GudhUI/utils/MClock.h | 109 | ||||
-rw-r--r-- | src/GudhUI/utils/Persistence_compute.h | 138 | ||||
-rw-r--r-- | src/GudhUI/utils/Rips_builder.h | 99 | ||||
-rw-r--r-- | src/GudhUI/utils/UI_utils.h | 43 | ||||
-rw-r--r-- | src/GudhUI/utils/Vertex_collapsor.h | 138 |
12 files changed, 880 insertions, 799 deletions
diff --git a/src/GudhUI/utils/Critical_points.h b/src/GudhUI/utils/Critical_points.h index 92392d4a..3021a5fe 100644 --- a/src/GudhUI/utils/Critical_points.h +++ b/src/GudhUI/utils/Critical_points.h @@ -1,13 +1,10 @@ -/* - * Critical_points.h - * Created on: Jan 27, 2015 - * This file is part of the Gudhi Library. The Gudhi library +/* 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): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 @@ -24,11 +21,13 @@ * */ - -#ifndef CRITICAL_POINTS_H_ -#define CRITICAL_POINTS_H_ +#ifndef UTILS_CRITICAL_POINTS_H_ +#define UTILS_CRITICAL_POINTS_H_ #include <deque> +#include <utility> // for pair<> +#include <algorithm> // for sort + #include "utils/Edge_contractor.h" /** @@ -37,101 +36,97 @@ * * todo do a sparsification with some parameter eps while growing */ -template<typename SkBlComplex> class Critical_points{ -private: - SkBlComplex filled_complex_; - const SkBlComplex& input_complex_; - double max_length_; - std::ostream& stream_; -public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - typedef typename std::pair<Vertex_handle,Vertex_handle> Edge; - - /** - * @brief check all pair of points with length smaller than max_length - */ - Critical_points(const SkBlComplex& input_complex,std::ostream& stream,double max_length): - input_complex_(input_complex),max_length_(max_length),stream_(stream) - { - - std::deque<Edge> edges; - auto vertices = input_complex.vertex_range(); - for(auto p = vertices.begin(); p!= vertices.end(); ++p){ - filled_complex_.add_vertex(input_complex.point(*p)); - for (auto q = p; ++q != vertices.end(); /**/) - if (squared_eucl_distance(input_complex.point(*p),input_complex.point(*q)) < max_length_*max_length_) - edges.emplace_back(*p,*q); - } - - std::sort(edges.begin(),edges.end(), - [&](Edge e1,Edge e2){ - return squared_edge_length(e1) < squared_edge_length(e2); - }); - - anti_collapse_edges(edges); - - } - -private: - double squared_eucl_distance(const Point& p1,const Point& p2) const{ - return Geometry_trait::Squared_distance_d()(p1,p2); - } - - - void anti_collapse_edges(const std::deque<Edge>& edges){ - unsigned pos = 0; - for(Edge e : edges){ - std::cout<<"edge "<<pos++<<"/"<<edges.size()<<"\n"; - auto eh = filled_complex_.add_edge(e.first,e.second); - int is_contractible(is_link_reducible(eh)); - - switch (is_contractible) { - case 0: - stream_<<"alpha="<<std::sqrt(squared_edge_length(e))<< " topological change"<<std::endl; - break; - case 2: - stream_<<"alpha="<<std::sqrt(squared_edge_length(e))<< " maybe a topological change"<<std::endl; - break; - default: - break; - } - } - - } - - - //0 -> not - //1 -> yes - //2 -> maybe - int is_link_reducible(Edge_handle e){ - auto link = filled_complex_.link(e); - - if(link.empty()) return 0; - - Edge_contractor<Complex> contractor(link,link.num_vertices()-1); - - if(link.num_connected_components()>1) return 0; //one than more CC -> not contractible - - if (link.num_vertices()==1) return 1; //reduced to one point -> contractible - else return 2; //we dont know - } - - - double squared_edge_length(Edge_handle e) const{ - return squared_eucl_distance(input_complex_.point(input_complex_.first_vertex(e)),input_complex_.point(input_complex_.second_vertex(e))); - } - - double squared_edge_length(Edge e) const{ - return squared_eucl_distance(input_complex_.point(e.first),input_complex_.point(e.second)); - } - - - +template<typename SkBlComplex> class Critical_points { + private: + SkBlComplex filled_complex_; + const SkBlComplex& input_complex_; + double max_length_; + std::ostream& stream_; + + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + typedef typename std::pair<Vertex_handle, Vertex_handle> Edge; + + /** + * @brief check all pair of points with length smaller than max_length + */ + Critical_points(const SkBlComplex& input_complex, std::ostream& stream, double max_length) : + input_complex_(input_complex), max_length_(max_length), stream_(stream) { + std::deque<Edge> edges; + auto vertices = input_complex.vertex_range(); + for (auto p = vertices.begin(); p != vertices.end(); ++p) { + filled_complex_.add_vertex(input_complex.point(*p)); + for (auto q = p; ++q != vertices.end(); /**/) + if (squared_eucl_distance(input_complex.point(*p), input_complex.point(*q)) < max_length_ * max_length_) + edges.emplace_back(*p, *q); + } + + std::sort(edges.begin(), edges.end(), + [&](Edge e1, Edge e2) { + return squared_edge_length(e1) < squared_edge_length(e2); + }); + + anti_collapse_edges(edges); + } + + private: + double squared_eucl_distance(const Point& p1, const Point& p2) const { + return Geometry_trait::Squared_distance_d()(p1, p2); + } + + void anti_collapse_edges(const std::deque<Edge>& edges) { + unsigned pos = 0; + for (Edge e : edges) { + std::cout << "edge " << pos++ << "/" << edges.size() << "\n"; + auto eh = filled_complex_.add_edge(e.first, e.second); + int is_contractible(is_link_reducible(eh)); + + switch (is_contractible) { + case 0: + stream_ << "alpha=" << std::sqrt(squared_edge_length(e)) << " topological change" << std::endl; + break; + case 2: + stream_ << "alpha=" << std::sqrt(squared_edge_length(e)) << " maybe a topological change" << std::endl; + break; + default: + break; + } + } + } + + // 0 -> not + // 1 -> yes + // 2 -> maybe + + int is_link_reducible(Edge_handle e) { + auto link = filled_complex_.link(e); + + 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; + + if (link.num_vertices() == 1) + // reduced to one point -> contractible + return 1; + else + // we dont know + return 2; + } + + double squared_edge_length(Edge_handle e) const { + return squared_eucl_distance(input_complex_.point(input_complex_.first_vertex(e)), + input_complex_.point(input_complex_.second_vertex(e))); + } + + double squared_edge_length(Edge e) const { + return squared_eucl_distance(input_complex_.point(e.first), input_complex_.point(e.second)); + } }; - - - - -#endif /* CRITICAL_POINTS_H_ */ +#endif // UTILS_CRITICAL_POINTS_H_ diff --git a/src/GudhUI/utils/Edge_collapsor.h b/src/GudhUI/utils/Edge_collapsor.h index 4dcd18ac..fcf2a248 100644 --- a/src/GudhUI/utils/Edge_collapsor.h +++ b/src/GudhUI/utils/Edge_collapsor.h @@ -1,12 +1,28 @@ -/* - * Collapsor.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 25, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 COLLAPSOR_H_ -#define COLLAPSOR_H_ +#ifndef UTILS_COLLAPSOR_H_ +#define UTILS_COLLAPSOR_H_ #include <list> #include "utils/Edge_contractor.h" @@ -15,71 +31,67 @@ /** * Iteratively puts every vertex at the center of its neighbors */ -template<typename SkBlComplex> class Edge_collapsor{ -private: - SkBlComplex& complex_; - unsigned num_collapses_; -public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - - /** - * @brief Collapse num_collapses edges. If num_collapses<0 then it collapses all possible edges. - * Largest edges are collapsed first. - * - */ - Edge_collapsor(SkBlComplex& complex,unsigned num_collapses): - complex_(complex),num_collapses_(num_collapses) - { - std::list<Edge_handle> edges; - edges.insert(edges.begin(),complex_.edge_range().begin(),complex_.edge_range().end()); - - edges.sort( - [&](Edge_handle e1,Edge_handle e2){ - return squared_edge_length(e1) < squared_edge_length(e2); - }); - - collapse_edges(edges); - - } - -private: - - void collapse_edges(std::list<Edge_handle>& edges){ - while(!edges.empty() && num_collapses_--){ - Edge_handle current_edge = edges.front(); - edges.pop_front(); - if(is_link_reducible(current_edge)) - complex_.remove_edge(current_edge); - } - - } - - bool is_link_reducible(Edge_handle e){ - auto link = complex_.link(e); - - if(link.empty()) return false; - - if(link.is_cone()) return true; - - if(link.num_connected_components()>1) return false; - - Edge_contractor<SkBlComplex> contractor(link,link.num_vertices()-1); - - return (link.num_vertices()==1); - } - - - double squared_edge_length(Edge_handle e) const{ - return squared_eucl_distance(complex_.point(complex_.first_vertex(e)),complex_.point(complex_.second_vertex(e))); - } - - double squared_eucl_distance(const Point& p1,const Point& p2) const{ - return Geometry_trait::Squared_distance_d()(p1,p2); - } - +template<typename SkBlComplex> class Edge_collapsor { + private: + SkBlComplex& complex_; + unsigned num_collapses_; + + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + + /** + * @brief Collapse num_collapses edges. If num_collapses<0 then it collapses all possible edges. + * Largest edges are collapsed first. + * + */ + Edge_collapsor(SkBlComplex& complex, unsigned num_collapses) : + complex_(complex), num_collapses_(num_collapses) { + std::list<Edge_handle> edges; + edges.insert(edges.begin(), complex_.edge_range().begin(), complex_.edge_range().end()); + + edges.sort( + [&](Edge_handle e1, Edge_handle e2) { + return squared_edge_length(e1) < squared_edge_length(e2); + }); + + collapse_edges(edges); + } + + private: + void collapse_edges(std::list<Edge_handle>& edges) { + while (!edges.empty() && num_collapses_--) { + Edge_handle current_edge = edges.front(); + edges.pop_front(); + if (is_link_reducible(current_edge)) + complex_.remove_edge(current_edge); + } + } + + bool is_link_reducible(Edge_handle e) { + auto link = complex_.link(e); + + if (link.empty()) + return false; + + if (link.is_cone()) + return true; + + if (link.num_connected_components() > 1) + return false; + + Edge_contractor<SkBlComplex> contractor(link, link.num_vertices() - 1); + + return (link.num_vertices() == 1); + } + + double squared_edge_length(Edge_handle e) const { + return squared_eucl_distance(complex_.point(complex_.first_vertex(e)), complex_.point(complex_.second_vertex(e))); + } + + double squared_eucl_distance(const Point& p1, const Point& p2) const { + return Geometry_trait::Squared_distance_d()(p1, p2); + } }; - - -#endif /* COLLAPSOR_H_ */ +#endif // UTILS_COLLAPSOR_H_ diff --git a/src/GudhUI/utils/Edge_contractor.h b/src/GudhUI/utils/Edge_contractor.h index d0a1f866..053aac04 100644 --- a/src/GudhUI/utils/Edge_contractor.h +++ b/src/GudhUI/utils/Edge_contractor.h @@ -1,85 +1,97 @@ -/* - * Contractor.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 25, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 EDGE_CONTRACTOR_H_ -#define EDGE_CONTRACTOR_H_ +#ifndef UTILS_EDGE_CONTRACTOR_H_ +#define UTILS_EDGE_CONTRACTOR_H_ -#include "gudhi/Skeleton_blocker_contractor.h" - -#include "gudhi/Contraction/Edge_profile.h" -#include "gudhi/Contraction/policies/Cost_policy.h" +#include <gudhi/Skeleton_blocker_contractor.h> +#include <gudhi/Contraction/Edge_profile.h> +#include <gudhi/Contraction/policies/Cost_policy.h> +#include <vector> /** * Iteratively puts every vertex at the center of its neighbors */ -template<typename SkBlComplex> class Edge_contractor{ -private: - SkBlComplex& complex_; - unsigned num_contractions_; - - /** - * @brief return a cost corresponding to the squared length of the edge - */ - template< typename EdgeProfile> class Length_cost : public contraction::Cost_policy<EdgeProfile>{ - public: - typedef typename contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; - typedef typename EdgeProfile::Point Point; - Cost_type operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const override{ - Cost_type res; - if(!placement) return res; - return Geometry_trait::Squared_distance_d()(profile.p0(),profile.p1()); //not working?? - } - }; - - - /** - * @brief return a cost corresponding to the squared length of the edge - */ - template< typename EdgeProfile> class Middle_placement : public contraction::Placement_policy<EdgeProfile>{ - public: - typedef typename contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; - typedef typename EdgeProfile::Point Point; - Placement_type operator()(const EdgeProfile& profile) const override{ - std::vector<double> mid_coords(profile.p0().dimension(),0); - for (size_t i = 0; i < profile.p0().dimension(); ++i){ - mid_coords[i] = (profile.p0()[i] + profile.p1()[i]) / 2.; - } - return Point(profile.p0().dimension(),mid_coords.begin(), mid_coords.end()); - } - }; - - public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - - /** - * @brief Modify complex to be the expansion of the k-nearest neighbor - * symetric graph. - */ - Edge_contractor(SkBlComplex& complex,unsigned num_contractions): - complex_(complex),num_contractions_(num_contractions) - { - typedef typename contraction::Edge_profile<Complex> Profile; - num_contractions = (std::min)((int)num_contractions,(int)(complex_.num_vertices()-1)); - contraction::Skeleton_blocker_contractor<Complex> contractor( - complex_, - new Length_cost<contraction::Edge_profile<Complex>>(), - new Middle_placement<contraction::Edge_profile<Complex>>(), - contraction::make_link_valid_contraction<Profile>(), - contraction::make_remove_popable_blockers_visitor<Profile>() - ); - contractor.contract_edges(num_contractions); - } - - +template<typename SkBlComplex> class Edge_contractor { + private: + SkBlComplex& complex_; + unsigned num_contractions_; + + /** + * @brief return a cost corresponding to the squared length of the edge + */ + template< typename EdgeProfile> class Length_cost : public contraction::Cost_policy<EdgeProfile> { + public: + typedef typename contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; + typedef typename EdgeProfile::Point Point; + + Cost_type operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const override { + Cost_type res; + if (!placement) + return res; + return Geometry_trait::Squared_distance_d()(profile.p0(), profile.p1()); // not working?? + } + }; + + /** + * @brief return a cost corresponding to the squared length of the edge + */ + template< typename EdgeProfile> class Middle_placement : public contraction::Placement_policy<EdgeProfile> { + public: + typedef typename contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; + typedef typename EdgeProfile::Point Point; + + Placement_type operator()(const EdgeProfile& profile) const override { + std::vector<double> mid_coords(profile.p0().dimension(), 0); + for (size_t i = 0; i < profile.p0().dimension(); ++i) { + mid_coords[i] = (profile.p0()[i] + profile.p1()[i]) / 2.; + } + return Point(profile.p0().dimension(), mid_coords.begin(), mid_coords.end()); + } + }; + + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + + /** + * @brief Modify complex to be the expansion of the k-nearest neighbor + * symetric graph. + */ + Edge_contractor(SkBlComplex& complex, unsigned num_contractions) : + complex_(complex), num_contractions_(num_contractions) { + typedef typename contraction::Edge_profile<Complex> Profile; + num_contractions = (std::min)(static_cast<int>(num_contractions), static_cast<int>(complex_.num_vertices() - 1)); + typedef typename contraction::Skeleton_blocker_contractor<Complex> Contractor; + Contractor contractor(complex_, + new Length_cost<contraction::Edge_profile < Complex >> (), + new Middle_placement<contraction::Edge_profile < Complex >> (), + contraction::make_link_valid_contraction<Profile>(), + contraction::make_remove_popable_blockers_visitor<Profile>()); + contractor.contract_edges(num_contractions); + } }; - - #endif /* EDGE_CONTRACTOR_H_ */ diff --git a/src/GudhUI/utils/Furthest_point_epsilon_net.h b/src/GudhUI/utils/Furthest_point_epsilon_net.h index 590b65c4..f2a216f6 100644 --- a/src/GudhUI/utils/Furthest_point_epsilon_net.h +++ b/src/GudhUI/utils/Furthest_point_epsilon_net.h @@ -1,134 +1,132 @@ -/* - * Furthest_point_epsilon_net.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 26, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 FURTHEST_POINT_EPSILON_NET_H_ -#define FURTHEST_POINT_EPSILON_NET_H_ +#ifndef UTILS_FURTHEST_POINT_EPSILON_NET_H_ +#define UTILS_FURTHEST_POINT_EPSILON_NET_H_ -#include "utils/UI_utils.h" #include <vector> +#include <algorithm> // for sort + +#include "utils/UI_utils.h" /** * Computes an epsilon net with furthest point strategy. */ -template<typename SkBlComplex> class Furthest_point_epsilon_net{ -private: - SkBlComplex& complex_; - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - - /** - * Let V be the set of vertices. - * Initially v0 is one arbitrarly vertex and the set V0 is {v0}. - * Then Vk is computed as follows. - * First we compute the vertex pk that is the furthest from Vk - * then Vk = Vk \cup pk. - * The radius of pk is its distance to Vk and its meeting vertex - * is the vertex of Vk for which this distance is achieved. - */ - struct Net_filtration_vertex{ - Vertex_handle vertex_handle; - Vertex_handle meeting_vertex; - double radius; - - - Net_filtration_vertex( - Vertex_handle vertex_handle_, - Vertex_handle meeting_vertex_, - double radius_): - vertex_handle(vertex_handle_),meeting_vertex(meeting_vertex_),radius(radius_) - {} - - bool operator<(const Net_filtration_vertex& other ) const{ - return radius < other.radius; - } - - }; - -public: - - - std::vector<Net_filtration_vertex> net_filtration_; - - /** - * @brief Modify complex to be the expansion of the k-nearest neighbor - * symetric graph. - */ - Furthest_point_epsilon_net(SkBlComplex& complex): - complex_(complex) - { - if(!complex.empty()){ - init_filtration(); - for(int k = 2; k < net_filtration_.size(); ++k){ - update_radius_value(k); - } - } - } - - //xxx does not work if complex not full - double radius(Vertex_handle v){ - return net_filtration_[v.vertex].radius; - } - - - - -private: - - void init_filtration(){ - Vertex_handle v0 = *(complex_.vertex_range().begin()); - net_filtration_.reserve(complex_.num_vertices()); - for(auto v : complex_.vertex_range()){ - if(v != v0) - net_filtration_.push_back( - Net_filtration_vertex(v, - Vertex_handle(-1), - squared_eucl_distance(v,v0)) - ); - } - net_filtration_.push_back(Net_filtration_vertex(v0,Vertex_handle(-1),1e10)); - auto n = net_filtration_.size(); - sort_filtration(n-1); - } - - void update_radius_value(int k){ - int n = net_filtration_.size(); - int index_to_update = n-k; - for(int i = 0; i< index_to_update; ++i){ - net_filtration_[i].radius = (std::min)(net_filtration_[i].radius , - squared_eucl_distance( - net_filtration_[i].vertex_handle, - net_filtration_[index_to_update].vertex_handle - ) - ); - } - sort_filtration(n-k); - } - - /** - * sort all i first elements. - */ - void sort_filtration(int i){ - std::sort(net_filtration_.begin(),net_filtration_.begin()+ i); - } - - double squared_eucl_distance(Vertex_handle v1,Vertex_handle v2) const{ - return std::sqrt(Geometry_trait::Squared_distance_d()( - complex_.point(v1),complex_.point(v2)) - ); - } - - void print_filtration() const{ - for(auto v : net_filtration_){ - std::cerr <<"v="<<v.vertex_handle<<"-> d="<<v.radius<<std::endl; - } - } - +template<typename SkBlComplex> class Furthest_point_epsilon_net { + private: + SkBlComplex& complex_; + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + + /** + * Let V be the set of vertices. + * Initially v0 is one arbitrarly vertex and the set V0 is {v0}. + * Then Vk is computed as follows. + * First we compute the vertex pk that is the furthest from Vk + * then Vk = Vk \cup pk. + * The radius of pk is its distance to Vk and its meeting vertex + * is the vertex of Vk for which this distance is achieved. + */ + struct Net_filtration_vertex { + Vertex_handle vertex_handle; + Vertex_handle meeting_vertex; + double radius; + + Net_filtration_vertex(Vertex_handle vertex_handle_, + Vertex_handle meeting_vertex_, + double radius_) : + vertex_handle(vertex_handle_), meeting_vertex(meeting_vertex_), radius(radius_) { } + + bool operator<(const Net_filtration_vertex& other) const { + return radius < other.radius; + } + }; + + public: + std::vector<Net_filtration_vertex> net_filtration_; + + /** + * @brief Modify complex to be the expansion of the k-nearest neighbor + * symetric graph. + */ + Furthest_point_epsilon_net(SkBlComplex& complex) : + complex_(complex) { + if (!complex.empty()) { + init_filtration(); + for (int k = 2; k < net_filtration_.size(); ++k) { + update_radius_value(k); + } + } + } + + // xxx does not work if complex not full + + double radius(Vertex_handle v) { + return net_filtration_[v.vertex].radius; + } + + private: + void init_filtration() { + Vertex_handle v0 = *(complex_.vertex_range().begin()); + net_filtration_.reserve(complex_.num_vertices()); + for (auto v : complex_.vertex_range()) { + if (v != v0) + net_filtration_.push_back(Net_filtration_vertex(v, + Vertex_handle(-1), + squared_eucl_distance(v, v0))); + } + net_filtration_.push_back(Net_filtration_vertex(v0, Vertex_handle(-1), 1e10)); + auto n = net_filtration_.size(); + sort_filtration(n - 1); + } + + void update_radius_value(int k) { + int n = net_filtration_.size(); + int index_to_update = n - k; + for (int i = 0; i < index_to_update; ++i) { + net_filtration_[i].radius = (std::min)(net_filtration_[i].radius, + squared_eucl_distance(net_filtration_[i].vertex_handle, + net_filtration_[index_to_update].vertex_handle)); + } + sort_filtration(n - k); + } + + /** + * sort all i first elements. + */ + void sort_filtration(int i) { + std::sort(net_filtration_.begin(), net_filtration_.begin() + i); + } + + double squared_eucl_distance(Vertex_handle v1, Vertex_handle v2) const { + return std::sqrt(Geometry_trait::Squared_distance_d()(complex_.point(v1), complex_.point(v2))); + } + + void print_filtration() const { + for (auto v : net_filtration_) { + std::cerr << "v=" << v.vertex_handle << "-> d=" << v.radius << std::endl; + } + } }; - - -#endif /* FURTHEST_POINT_EPSILON_NET_H_ */ +#endif // UTILS_FURTHEST_POINT_EPSILON_NET_H_ diff --git a/src/GudhUI/utils/Is_manifold.h b/src/GudhUI/utils/Is_manifold.h index e708a6a4..1b25df34 100644 --- a/src/GudhUI/utils/Is_manifold.h +++ b/src/GudhUI/utils/Is_manifold.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-Mediterranee (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 @@ -25,11 +25,10 @@ */ -#ifndef IS_MANIFOLD_H_ -#define IS_MANIFOLD_H_ +#ifndef UTILS_IS_MANIFOLD_H_ +#define UTILS_IS_MANIFOLD_H_ #include "utils/UI_utils.h" - #include "utils/Edge_contractor.h" /** @@ -38,67 +37,67 @@ * * todo do a sparsification with some parameter eps while growing */ -template<typename SkBlComplex> class Is_manifold{ -private: - const SkBlComplex& input_complex_; - typedef typename SkBlComplex::Vertex_handle Vertex_handle; +template<typename SkBlComplex> class Is_manifold { + private: + const SkBlComplex& input_complex_; + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + public: + /* + * return dim the maximum dimension around one simplex and res which is true if the complex is a manifold. + * If the complex has dimension <= 3 then if res is false, the complex is not a manifold. + * For d-manifold with d>=4, res may be false while the complex is a manifold. + */ + Is_manifold(const SkBlComplex& input_complex, unsigned& dim, bool & res) : input_complex_(input_complex) { + res = true; + dim = -1; + if (!input_complex_.empty()) { + for (auto v : input_complex_.vertex_range()) { + dim = local_dimension(v); + break; + } + //check that the link of every vertex is a dim-1 sphere + for (auto v : input_complex_.vertex_range()) { + if (!is_k_sphere(v, dim - 1)) { + res = false; + break; + } + } + } + } -public: - /* - * return dim the maximum dimension around one simplex and res which is true if the complex is a manifold. - * If the complex has dimension <= 3 then if res is false, the complex is not a manifold. - * For d-manifold with d>=4, res may be false while the complex is a manifold. - */ - Is_manifold(const SkBlComplex& input_complex,unsigned& dim,bool & res):input_complex_(input_complex){ - res = true; - dim = -1; - if(!input_complex_.empty()){ - for(auto v : input_complex_.vertex_range()){ - dim = local_dimension(v); - break; - } - //check that the link of every vertex is a dim-1 sphere - for(auto v : input_complex_.vertex_range()){ - if(!is_k_sphere(v,dim-1)) { - res = false; - break; - } - } - } - } + private: + unsigned local_dimension(Vertex_handle v) { + unsigned dim = 0; + for (const auto& s : input_complex_.simplex_range(v)) + dim = (std::max)(dim, (unsigned) s.dimension()); + return dim; + } -private: - unsigned local_dimension(Vertex_handle v){ - unsigned dim = 0; - for(const auto& s: input_complex_.simplex_range(v)) - dim = (std::max)(dim,(unsigned)s.dimension()); - return dim; - } + 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); + } - 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); - } + // A minimal sphere is a complex that contains vertices v1...vn and all faces + // made upon this set except the face {v1,...,vn} + // return -2 if not a minimal sphere + // and d otherwise if complex is a d minimal sphere - // A minimal sphere is a complex that contains vertices v1...vn and all faces - // made upon this set except the face {v1,...,vn} - // return -2 if not a minimal sphere - // and d otherwise if complex is a d minimal sphere - template<typename SubComplex> - int is_sphere_simplex(const SubComplex& complex){ - if(complex.empty()) return -1; - if(complex.num_blockers()!=1) return -2; + template<typename SubComplex> + int is_sphere_simplex(const SubComplex& complex) { + if (complex.empty()) return -1; + if (complex.num_blockers() != 1) return -2; - //necessary and sufficient condition : there exists a unique blocker that passes through all vertices - auto first_blocker = *(complex.const_blocker_range().begin()); + //necessary and sufficient condition : there exists a unique blocker that passes through all vertices + auto first_blocker = *(complex.const_blocker_range().begin()); - if (first_blocker->dimension()+1 != complex.num_vertices()) - return -2; - else - return (first_blocker->dimension()-1); - } + if (first_blocker->dimension() + 1 != complex.num_vertices()) + return -2; + else + return (first_blocker->dimension() - 1); + } }; -#endif /* IS_MANIFOLD_H_ */ +#endif // UTILS_IS_MANIFOLD_H_ diff --git a/src/GudhUI/utils/K_nearest_builder.h b/src/GudhUI/utils/K_nearest_builder.h index b65a79e0..cab24b7c 100644 --- a/src/GudhUI/utils/K_nearest_builder.h +++ b/src/GudhUI/utils/K_nearest_builder.h @@ -1,81 +1,92 @@ -/* - * K_nearest_builder.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 10, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 K_NEAREST_BUILDER_H_ -#define K_NEAREST_BUILDER_H_ +#ifndef UTILS_K_NEAREST_BUILDER_H_ +#define UTILS_K_NEAREST_BUILDER_H_ -#include <unordered_map> -#include <boost/iterator/iterator_facade.hpp> #include <CGAL/Euclidean_distance.h> #include <CGAL/Orthogonal_k_neighbor_search.h> #include <CGAL/Search_traits_d.h> #include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h>
+#include <CGAL/property_map.h> +#include <boost/iterator/iterator_facade.hpp> #include <boost/iterator/zip_iterator.hpp> +#include <unordered_map> +#include <tuple> +#include <list> + #include "utils/UI_utils.h" #include "model/Complex_typedefs.h" - - -template<typename SkBlComplex> class K_nearest_builder{ -private: - - typedef Geometry_trait Kernel; - typedef Point Point_d;
- typedef boost::tuple<Point_d, unsigned> Point_d_with_id;
- typedef CGAL::Search_traits_d<Kernel> Traits_base;
- typedef CGAL::Search_traits_adapter<Point_d_with_id, CGAL::Nth_of_tuple_property_map<0, Point_d_with_id>, Traits_base> Traits;
- typedef CGAL::Orthogonal_k_neighbor_search<Traits> Neighbor_search;
- typedef Neighbor_search::Tree Tree;
- typedef Neighbor_search::Distance Distance;
- - SkBlComplex& complex_; -public: - - /** - * @brief Modify complex to be the expansion of the k-nearest neighbor - * symetric graph. - */ - K_nearest_builder(SkBlComplex& complex,unsigned k):complex_(complex){ - complex.keep_only_vertices(); - compute_edges(k); - } - -private: - - - - double squared_eucl_distance(const Point& p1,const Point& p2) const{ - return Geometry_trait::Squared_distance_d()(p1,p2); - } - - void compute_edges(unsigned k){ - - std::list<Point_d_with_id> points_with_id; - for(auto v: complex_.vertex_range()){ - Point_d_with_id point_v(complex_.point(v),v.vertex); - points_with_id.push_back(point_v); - } - - Tree tree(points_with_id.begin(),points_with_id.end()); - - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - for (auto p : complex_.vertex_range()){ - Neighbor_search search(tree, complex_.point(p),k+1); - for(auto it = ++search.begin(); it != search.end(); ++it){ - Vertex_handle q(boost::get<1>(it->first)); - if (p != q && complex_.contains_vertex(p) && complex_.contains_vertex(q)) - complex_.add_edge(p,q); - } - } - } - - +template<typename SkBlComplex> class K_nearest_builder { + private: + typedef Geometry_trait Kernel; + typedef Point Point_d; + typedef boost::tuple<Point_d, unsigned> Point_d_with_id; + typedef CGAL::Search_traits_d<Kernel> Traits_base; + typedef CGAL::Search_traits_adapter<Point_d_with_id, CGAL::Nth_of_tuple_property_map<0, Point_d_with_id>, + Traits_base> Traits; + typedef CGAL::Orthogonal_k_neighbor_search<Traits> Neighbor_search; + typedef Neighbor_search::Tree Tree; + typedef Neighbor_search::Distance Distance; + + SkBlComplex& complex_; + + public: + /** + * @brief Modify complex to be the expansion of the k-nearest neighbor + * symetric graph. + */ + K_nearest_builder(SkBlComplex& complex, unsigned k) : complex_(complex) { + complex.keep_only_vertices(); + compute_edges(k); + } + + private: + double squared_eucl_distance(const Point& p1, const Point& p2) const { + return Geometry_trait::Squared_distance_d()(p1, p2); + } + + void compute_edges(unsigned k) { + std::list<Point_d_with_id> points_with_id; + for (auto v : complex_.vertex_range()) { + Point_d_with_id point_v(complex_.point(v), v.vertex); + points_with_id.push_back(point_v); + } + + Tree tree(points_with_id.begin(), points_with_id.end()); + + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + for (auto p : complex_.vertex_range()) { + Neighbor_search search(tree, complex_.point(p), k + 1); + for (auto it = ++search.begin(); it != search.end(); ++it) { + Vertex_handle q(boost::get<1>(it->first)); + if (p != q && complex_.contains_vertex(p) && complex_.contains_vertex(q)) + complex_.add_edge(p, q); + } + } + } }; -#endif /* K_NEAREST_BUILDER_H_ */ +#endif // UTILS_K_NEAREST_BUILDER_H_ diff --git a/src/GudhUI/utils/Lloyd_builder.h b/src/GudhUI/utils/Lloyd_builder.h index db2a7973..18ec9fac 100644 --- a/src/GudhUI/utils/Lloyd_builder.h +++ b/src/GudhUI/utils/Lloyd_builder.h @@ -1,78 +1,91 @@ -/* - * Lloyd.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 25, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 LLOYD_H_ -#define LLOYD_H_ - +#ifndef UTILS_LLOYD_BUILDER_H_ +#define UTILS_LLOYD_BUILDER_H_ #include <vector> +#include <list> /** * Iteratively puts every vertex at the center of its neighbors */ -template<typename SkBlComplex> class Lloyd_builder{ -private: - SkBlComplex& complex_; - int dim; -public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - /** - * @brief Modify complex to be the expansion of the k-nearest neighbor - * symetric graph. - */ - Lloyd_builder(SkBlComplex& complex,unsigned num_iterations):complex_(complex),dim(-1){ - if(!complex_.empty()){ - dim = get_dimension(); - while(num_iterations--){ - std::list<Point> new_points; - for(auto v : complex.vertex_range()) - new_points.push_back(barycenter_neighbors(v)); - - auto new_points_it = new_points.begin(); - for(auto v : complex.vertex_range()) - complex_.point(v) = *(new_points_it++); - } - } - } +template<typename SkBlComplex> class Lloyd_builder { + private: + SkBlComplex& complex_; + int dim; -private: + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; - int get_dimension(){ - assert(!complex_.empty()); - for(auto v : complex_.vertex_range()) - return complex_.point(v).dimension(); - return -1; - } + /** + * @brief Modify complex to be the expansion of the k-nearest neighbor + * symetric graph. + */ + Lloyd_builder(SkBlComplex& complex, unsigned num_iterations) : complex_(complex), dim(-1) { + if (!complex_.empty()) { + dim = get_dimension(); + while (num_iterations--) { + std::list<Point> new_points; + for (auto v : complex.vertex_range()) + new_points.push_back(barycenter_neighbors(v)); - Point barycenter_neighbors(Vertex_handle v) const{ - if(complex_.degree(v)==0) - return complex_.point(v); + auto new_points_it = new_points.begin(); + for (auto v : complex.vertex_range()) + complex_.point(v) = *(new_points_it++); + } + } + } - std::vector<double> res(dim,0); - unsigned num_points = 0; - for(auto nv : complex_.vertex_range(v)){ - ++num_points; - const Point& point = complex_.point(nv); - assert(point.dimension() == dim); - for(int i = 0; i < point.dimension() ; ++i) - res[i] += point[i]; - } - for(auto& x : res) - x/= num_points; - return Point(dim,res.begin(),res.end()); - } + private: + int get_dimension() { + assert(!complex_.empty()); + for (auto v : complex_.vertex_range()) + return complex_.point(v).dimension(); + return -1; + } + Point barycenter_neighbors(Vertex_handle v) const { + if (complex_.degree(v) == 0) + return complex_.point(v); + std::vector<double> res(dim, 0); + unsigned num_points = 0; + for (auto nv : complex_.vertex_range(v)) { + ++num_points; + const Point& point = complex_.point(nv); + assert(point.dimension() == dim); + for (int i = 0; i < point.dimension(); ++i) + res[i] += point[i]; + } + for (auto& x : res) + x /= num_points; + return Point(dim, res.begin(), res.end()); + } - double squared_eucl_distance(const Point& p1,const Point& p2) const{ - return Geometry_trait::Squared_distance_d()(p1,p2); - } - + double squared_eucl_distance(const Point& p1, const Point& p2) const { + return Geometry_trait::Squared_distance_d()(p1, p2); + } }; - -#endif /* LLOYD_H_ */ +#endif // UTILS_LLOYD_BUILDER_H_ diff --git a/src/GudhUI/utils/MClock.h b/src/GudhUI/utils/MClock.h index b5c7d191..e8d8918a 100644 --- a/src/GudhUI/utils/MClock.h +++ b/src/GudhUI/utils/MClock.h @@ -1,57 +1,70 @@ -/* - * Clock.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Jun 17, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 CLOCK_H_ -#define CLOCK_H_ +#ifndef UTILS_MCLOCK_H_ +#define UTILS_MCLOCK_H_ #include <sys/time.h> +class MClock { + public: + MClock() { + end_called = false; + begin(); + } -class MClock{ - -public: - MClock(){ - end_called = false; - begin(); - } - - void begin() const{ - end_called = false; - gettimeofday(&startTime, NULL); - } - - void end() const{ - end_called = true; - gettimeofday(&endTime, NULL); - } - - friend std::ostream& operator<< (std::ostream& stream,const MClock& clock){ - // if(!clock.end_called) clock.end(); - if(!clock.end_called) stream << "end not called"; - else{ - long totalTime = (clock.endTime.tv_sec - clock.startTime.tv_sec) * 1000000L; - totalTime += (clock.endTime.tv_usec - clock.startTime.tv_usec); - stream << clock.num_seconds() <<"s"; - } - return stream; - - } - - double num_seconds() const{ - if(!end_called) end(); - long totalTime = (endTime.tv_sec - startTime.tv_sec) * 1000000L; - totalTime += (endTime.tv_usec - startTime.tv_usec); - return (totalTime / 1000L)/(double) 1000; - } - -private: - mutable struct timeval startTime, endTime; - mutable bool end_called; -}; + void begin() const { + end_called = false; + gettimeofday(&startTime, NULL); + } + void end() const { + end_called = true; + gettimeofday(&endTime, NULL); + } + + friend std::ostream& operator<<(std::ostream& stream, const MClock& clock) { + // if(!clock.end_called) clock.end(); + if (!clock.end_called) { + stream << "end not called"; + } else { + long totalTime = (clock.endTime.tv_sec - clock.startTime.tv_sec) * 1000000L; + totalTime += (clock.endTime.tv_usec - clock.startTime.tv_usec); + stream << clock.num_seconds() << "s"; + } + return stream; + } + + double num_seconds() const { + if (!end_called) end(); + long totalTime = (endTime.tv_sec - startTime.tv_sec) * 1000000L; + totalTime += (endTime.tv_usec - startTime.tv_usec); + return (totalTime / 1000L) / static_cast<double>(1000); + } + + private: + mutable struct timeval startTime, endTime; + mutable bool end_called; +}; -#endif /* CLOCK_H_ */ +#endif // UTILS_MCLOCK_H_ diff --git a/src/GudhUI/utils/Persistence_compute.h b/src/GudhUI/utils/Persistence_compute.h index 50d340b8..0b9961d3 100644 --- a/src/GudhUI/utils/Persistence_compute.h +++ b/src/GudhUI/utils/Persistence_compute.h @@ -1,13 +1,10 @@ -/* - * Persistence_compute.h - * Created on: Jan 26, 2015 - * This file is part of the Gudhi Library. The Gudhi library +/* 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): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 @@ -25,83 +22,76 @@ */ -#ifndef PERSISTENCE_COMPUTE_H_ -#define PERSISTENCE_COMPUTE_H_ +#ifndef UTILS_PERSISTENCE_COMPUTE_H_ +#define UTILS_PERSISTENCE_COMPUTE_H_ +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Persistent_cohomology.h> -#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Simplex_tree.h" -#include "gudhi/distance_functions.h" -#include "gudhi/Persistent_cohomology.h" +#include <vector> +struct Persistence_params { + int p; + double threshold; + int max_dim; + double min_pers; -struct Persistence_params{ - int p; - double threshold; - int max_dim; - double min_pers; - - Persistence_params(int p_,double th_,int max_dim_=10,double min_pers_=0) - :p(p_),threshold(th_),max_dim(max_dim_),min_pers(min_pers_){} + Persistence_params(int p_, double th_, int max_dim_ = 10, double min_pers_ = 0) + : p(p_), threshold(th_), max_dim(max_dim_), min_pers(min_pers_) { } }; - /** * Show persistence into output stream */ -template<typename SkBlComplex> class Persistence_compute{ -private: - SkBlComplex& complex_; - std::ostream& stream_; -public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - - /** - * @brief Compute persistence - * parameters : - * unsigned dim_max - * double threshold - * int p for coefficient Z_p - */ - Persistence_compute(SkBlComplex& complex,std::ostream& stream,const Persistence_params& params): -// -// double threshold = 0.5,unsigned dim_max = 8): - complex_(complex),stream_(stream){ - - //for now everything is copied, todo boost adapt iterators to points of SkBlComplex instead of copying to an intial vector - typedef std::vector<double> Point_t; - std::vector< Point_t > points; - points.reserve(complex.num_vertices()); - for(auto v : complex.vertex_range()){ - const auto & pt = complex.point(v); - Point_t pt_to_add(pt.cartesian_begin(),pt.cartesian_end()); - points.emplace_back(std::move(pt_to_add)); - } - - - Graph_t prox_graph = compute_proximity_graph( points, params.threshold, euclidean_distance<Point_t> ); - Gudhi::Simplex_tree<> st; - st.insert_graph(prox_graph); - st.expansion( params.max_dim ); - - Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Simplex_tree<>,Gudhi::persistent_cohomology::Field_Zp > pcoh (st); - pcoh.init_coefficients( params.p ); //initilizes the coefficient field for homology - pcoh.compute_persistent_cohomology( params.min_pers ); //put params.min_pers - stream_<<"persistence: \n"; - stream_<<"p dimension birth death: \n"; - - pcoh.output_diagram(stream_); - } - -private: - - +template<typename SkBlComplex> class Persistence_compute { + private: + SkBlComplex& complex_; + std::ostream& stream_; + + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + + /** + * @brief Compute persistence + * parameters : + * unsigned dim_max + * double threshold + * int p for coefficient Z_p + */ + Persistence_compute(SkBlComplex& complex, std::ostream& stream, const Persistence_params& params) : + // double threshold = 0.5,unsigned dim_max = 8): + complex_(complex), stream_(stream) { + // for now everything is copied, todo boost adapt iterators to points of SkBlComplex instead of copying to an + // initial vector + typedef std::vector<double> Point_t; + std::vector< Point_t > points; + points.reserve(complex.num_vertices()); + for (auto v : complex.vertex_range()) { + const auto & pt = complex.point(v); + Point_t pt_to_add(pt.cartesian_begin(), pt.cartesian_end()); + points.emplace_back(std::move(pt_to_add)); + } + + + Graph_t prox_graph = compute_proximity_graph(points, params.threshold, euclidean_distance<Point_t>); + Gudhi::Simplex_tree<> st; + st.insert_graph(prox_graph); + st.expansion(params.max_dim); + + Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Simplex_tree<>, + Gudhi::persistent_cohomology::Field_Zp > pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(params.p); + // put params.min_pers + pcoh.compute_persistent_cohomology(params.min_pers); + stream_ << "persistence: \n"; + stream_ << "p dimension birth death: \n"; + + pcoh.output_diagram(stream_); + } }; - - - - - -#endif /* PERSISTENCE_COMPUTE_H_ */ +#endif // UTILS_PERSISTENCE_COMPUTE_H_ diff --git a/src/GudhUI/utils/Rips_builder.h b/src/GudhUI/utils/Rips_builder.h index 9484f9ab..b22f4db6 100644 --- a/src/GudhUI/utils/Rips_builder.h +++ b/src/GudhUI/utils/Rips_builder.h @@ -1,56 +1,69 @@ -/* - * Rips_builder.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 10, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 RIPS_BUILDER_H_ -#define RIPS_BUILDER_H_ +#ifndef UTILS_RIPS_BUILDER_H_ +#define UTILS_RIPS_BUILDER_H_ #include <boost/iterator/iterator_facade.hpp> -#include "utils/UI_utils.h" -#include "model/Complex_typedefs.h" #include <CGAL/Euclidean_distance.h> - #include <CGAL/Orthogonal_k_neighbor_search.h> #include <CGAL/Search_traits_d.h> -template<typename SkBlComplex> class Rips_builder{ -private: - SkBlComplex& complex_; -public: - - /** - * @brief Modify complex to be the Rips complex - * of its points with offset alpha. - */ - Rips_builder(SkBlComplex& complex,double alpha):complex_(complex){ - complex.keep_only_vertices(); - if (alpha<=0) return; - compute_edges(alpha); - } - -private: - - - double squared_eucl_distance(const Point& p1,const Point& p2) const{ - return Geometry_trait::Squared_distance_d()(p1,p2); - } - - void compute_edges(double alpha){ - auto vertices = complex_.vertex_range(); - for(auto p = vertices.begin(); p!= vertices.end(); ++p){ - std::cout << *p << " "; std::cout.flush(); - for (auto q = p; ++q != vertices.end(); /**/) - if (squared_eucl_distance(complex_.point(*p),complex_.point(*q)) < 4*alpha*alpha) - complex_.add_edge(*p,*q); - } - std::cout << std::endl; - } +#include "utils/UI_utils.h" +#include "model/Complex_typedefs.h" -}; +template<typename SkBlComplex> class Rips_builder { + private: + SkBlComplex& complex_; + public: + /** + * @brief Modify complex to be the Rips complex + * of its points with offset alpha. + */ + Rips_builder(SkBlComplex& complex, double alpha) : complex_(complex) { + complex.keep_only_vertices(); + if (alpha <= 0) return; + compute_edges(alpha); + } + + private: + double squared_eucl_distance(const Point& p1, const Point& p2) const { + return Geometry_trait::Squared_distance_d()(p1, p2); + } + + void compute_edges(double alpha) { + auto vertices = complex_.vertex_range(); + for (auto p = vertices.begin(); p != vertices.end(); ++p) { + std::cout << *p << " "; + std::cout.flush(); + for (auto q = p; ++q != vertices.end(); /**/) + if (squared_eucl_distance(complex_.point(*p), complex_.point(*q)) < 4 * alpha * alpha) + complex_.add_edge(*p, *q); + } + std::cout << std::endl; + } +}; -#endif /* RIPS_BUILDER_H_ */ +#endif // UTILS_RIPS_BUILDER_H_ diff --git a/src/GudhUI/utils/UI_utils.h b/src/GudhUI/utils/UI_utils.h index 4ade4b98..9cc209d3 100644 --- a/src/GudhUI/utils/UI_utils.h +++ b/src/GudhUI/utils/UI_utils.h @@ -1,30 +1,45 @@ -/* - * UI_utils.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Aug 25, 2014 - * Author: david + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 UI_UTILS_H_ -#define UI_UTILS_H_ +#ifndef UTILS_UI_UTILS_H_ +#define UTILS_UI_UTILS_H_ #define UIDBG_VERBOSE #ifdef UIDBG_VERBOSE -#define UIDBG(a) std::cerr << "UIDBG: " << (a)<<std::endl -#define UIDBGMSG(a,b) std::cerr << "UIDBG: " << a<<b<<std::endl -#define UIDBGVALUE(a) std::cerr << "UIDBG: " << #a << ": " << a<<std::endl -#define UIDBGCONT(a) std::cerr << "UIDBG: container "<< #a<<" -> "; for(auto x:a) std::cerr<< x << ","; std::cerr<<std::endl +#define UIDBG(a) std::cerr << "UIDBG: " << (a) << std::endl +#define UIDBGMSG(a, b) std::cerr << "UIDBG: " << a << b << std::endl +#define UIDBGVALUE(a) std::cerr << "UIDBG: " << #a << ": " << a << std::endl +#define UIDBGCONT(a) std::cerr << "UIDBG: container " << #a << " -> "; for (auto x : a) std::cerr << x << ","; std::cerr << std::endl } #else // #define DBG(a) a -// #define DBGMSG(a,b) b +// #define DBGMSG(a, b) b // #define DBGVALUE(a) a // #define DBGCONT(a) a #define UIDBG(a) -#define UIDBGMSG(a,b) +#define UIDBGMSG(a, b) #define UIDBGVALUE(a) #define UIDBGCONT(a) #endif - -#endif // UI_UTILS_H_ +#endif // UTILS_UI_UTILS_H_ diff --git a/src/GudhUI/utils/Vertex_collapsor.h b/src/GudhUI/utils/Vertex_collapsor.h index d4911a35..be03c765 100644 --- a/src/GudhUI/utils/Vertex_collapsor.h +++ b/src/GudhUI/utils/Vertex_collapsor.h @@ -1,76 +1,86 @@ -/* - * Vertex_collapsor.h +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. * - * Created on: Sep 25, 2014 - * Author: dsalinas + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 VERTEX_COLLAPSOR_H_ -#define VERTEX_COLLAPSOR_H_ +#ifndef UTILS_VERTEX_COLLAPSOR_H_ +#define UTILS_VERTEX_COLLAPSOR_H_ #include "utils/Edge_contractor.h" #include "utils/Furthest_point_epsilon_net.h" #include "utils/UI_utils.h" + /** * Iteratively puts every vertex at the center of its neighbors */ -template<typename SkBlComplex> class Vertex_collapsor{ -private: - SkBlComplex& complex_; - size_t num_collapses_; -public: - typedef typename SkBlComplex::Vertex_handle Vertex_handle; - typedef typename SkBlComplex::Edge_handle Edge_handle; - - /** - * @brief Modify complex to be the expansion of the k-nearest neighbor - * symetric graph. - */ - Vertex_collapsor(SkBlComplex& complex, size_t num_collapses) : - complex_(complex),num_collapses_(num_collapses) - { -// std::list<Vertex_handle> vertices; -// vertices.insert(vertices.begin(),complex_.vertex_range().begin(),complex_.vertex_range().end()); -// UIDBG("Collapse vertices"); -// collapse_vertices(vertices); - - std::list<Vertex_handle> vertices; - - UIDBG("Compute eps net"); - Furthest_point_epsilon_net<Complex> eps_net(complex_); - - for(auto vh : eps_net.net_filtration_) - vertices.push_back(vh.vertex_handle); - - UIDBG("Collapse vertices"); - collapse_vertices(vertices); - - - - } - -private: - - - void collapse_vertices(std::list<Vertex_handle>& vertices){ - while(!vertices.empty() && num_collapses_--){ - Vertex_handle current_vertex = vertices.front(); - vertices.pop_front(); - if(is_link_reducible(current_vertex)) - complex_.remove_vertex(current_vertex); - } - } - - bool is_link_reducible(Vertex_handle v){ - auto link = complex_.link(v); - 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); - } - +template<typename SkBlComplex> class Vertex_collapsor { + private: + SkBlComplex& complex_; + size_t num_collapses_; + + public: + typedef typename SkBlComplex::Vertex_handle Vertex_handle; + typedef typename SkBlComplex::Edge_handle Edge_handle; + + /** + * @brief Modify complex to be the expansion of the k-nearest neighbor + * symetric graph. + */ + Vertex_collapsor(SkBlComplex& complex, size_t num_collapses) : + complex_(complex), num_collapses_(num_collapses) { + // std::list<Vertex_handle> vertices; + // vertices.insert(vertices.begin(),complex_.vertex_range().begin(),complex_.vertex_range().end()); + // UIDBG("Collapse vertices"); + // collapse_vertices(vertices); + + std::list<Vertex_handle> vertices; + + UIDBG("Compute eps net"); + Furthest_point_epsilon_net<Complex> eps_net(complex_); + + for (auto vh : eps_net.net_filtration_) + vertices.push_back(vh.vertex_handle); + + UIDBG("Collapse vertices"); + collapse_vertices(vertices); + } + + private: + void collapse_vertices(std::list<Vertex_handle>& vertices) { + while (!vertices.empty() && num_collapses_--) { + Vertex_handle current_vertex = vertices.front(); + vertices.pop_front(); + if (is_link_reducible(current_vertex)) + complex_.remove_vertex(current_vertex); + } + } + + bool is_link_reducible(Vertex_handle v) { + auto link = complex_.link(v); + 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); + } }; - -#endif /* VERTEX_COLLAPSOR_H_ */ +#endif // UTILS_VERTEX_COLLAPSOR_H_ |