/* 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 * * 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 . */ #ifndef SKELETON_BLOCKER_CONTRACTOR_H_ #define SKELETON_BLOCKER_CONTRACTOR_H_ // todo remove the queue to be independent from cgald #include #include #include #include #include #include #include #include // xxx remove #include // xxx remove #include #include #include #include #include #include #include #include #include // for pair #include namespace Gudhi { namespace contraction { template Placement_policy* make_first_vertex_placement() { return new First_vertex_placement(); } template Valid_contraction_policy* make_link_valid_contraction() { return new Link_condition_valid_contraction(); } /** *@brief Visitor to remove popable blockers after an edge contraction. */ template class Contraction_visitor_remove_popable : public Contraction_visitor { public: typedef typename Profile::Point Point; typedef typename Profile::Complex Complex; typedef typename Complex::Vertex_handle Vertex_handle; void on_contracted(const Profile &profile, boost::optional< Point > placement) override { profile.complex().remove_all_popable_blockers(profile.v0_handle()); } }; template Contraction_visitor* make_remove_popable_blockers_visitor() { return new Contraction_visitor_remove_popable(); } /** *@class Skeleton_blocker_contractor *@brief Class that allows to contract iteratively edges of a simplicial complex. *@ingroup contr * * @details The simplification algorithm consists in iteratively picking the * edge with lowest cost and performing an edge contraction if the contraction is valid. * This class is policy based (and much inspired from the edge collapse package of CGAL http://doc.cgal.org/latest/Surface_mesh_simplification/index.html). * * Policies that can be changed are : * - the cost policy : how much cost an edge contraction * - the placement policy : where will be placed the contraction point * - the valid contraction policy : is the contraction valid. For instance, it can be * a topological condition (link condition) or a geometrical condition (normals messed up). * */ template> class Skeleton_blocker_contractor : private skeleton_blocker::Dummy_complex_visitor< typename GeometricSimplifiableComplex::Vertex_handle> { GeometricSimplifiableComplex& complex_; public: typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex; typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle; typedef typename GeometricSimplifiableComplex::Simplex Simplex; typedef typename GeometricSimplifiableComplex::Root_vertex_handle Root_vertex_handle; typedef typename GeometricSimplifiableComplex::Graph_edge Graph_edge; typedef typename GeometricSimplifiableComplex::Edge_handle Edge_handle; typedef typename GeometricSimplifiableComplex::Point Point; typedef EdgeProfile Profile; typedef Cost_policy Cost_policy_; typedef Placement_policy Placement_policy_; typedef Valid_contraction_policy Valid_contraction_policy_; typedef Contraction_visitor Contraction_visitor_; typedef Edge_profile_factory Edge_profile_factory_; typedef boost::optional Cost_type; typedef boost::optional Placement_type; typedef size_t size_type; typedef Skeleton_blocker_contractor Self; private: struct Compare_id { Compare_id() : algorithm_(0) { } Compare_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } bool operator()(Edge_handle a, Edge_handle b) const { return algorithm_->get_undirected_edge_id(a) < algorithm_->get_undirected_edge_id(b); } Self const* algorithm_; }; struct Compare_cost { Compare_cost() : algorithm_(0) { } Compare_cost(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } bool operator()(Edge_handle a, Edge_handle b) const { // NOTE: A cost is an optional<> value. // Absent optionals are ordered first; that is, "none < T" and "T > none" for any defined T != none. // In consequence, edges with undefined costs will be promoted to the top of the priority queue and popped out // first. return algorithm_->get_data(a).cost() < algorithm_->get_data(b).cost(); } Self const* algorithm_; }; struct Undirected_edge_id : boost::put_get_helper { typedef boost::readable_property_map_tag category; typedef size_type value_type; typedef size_type reference; typedef Edge_handle key_type; Undirected_edge_id() : algorithm_(0) { } Undirected_edge_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } size_type operator[](Edge_handle e) const { return algorithm_->get_undirected_edge_id(e); } Self const* algorithm_; }; typedef CGAL::Modifiable_priority_queue PQ; typedef typename PQ::handle pq_handle; // An Edge_data is associated with EVERY edge in the complex (collapsible or not). // It relates the edge with the PQ-handle needed to update the priority queue // It also relates the edge with a policy-based cache class Edge_data { public: Edge_data() : PQHandle_(), cost_() { } Cost_type const& cost() const { return cost_; } Cost_type & cost() { return cost_; } pq_handle PQ_handle() const { return PQHandle_; } bool is_in_PQ() const { return PQHandle_ != PQ::null_handle(); } void set_PQ_handle(pq_handle h) { PQHandle_ = h; } void reset_PQ_handle() { PQHandle_ = PQ::null_handle(); } private: pq_handle PQHandle_; Cost_type cost_; }; typedef Edge_data* Edge_data_ptr; typedef boost::scoped_array Edge_data_array; int get_undirected_edge_id(Edge_handle edge) const { return complex_[edge].index(); } const Edge_data& get_data(Edge_handle edge) const { return edge_data_array_[get_undirected_edge_id(edge)]; } Edge_data& get_data(Edge_handle edge) { return edge_data_array_[get_undirected_edge_id(edge)]; } Cost_type get_cost(const Profile & profile) const { return (*cost_policy_)(profile, get_placement(profile)); } Profile create_profile(Edge_handle edge) const { if (edge_profile_factory_) return edge_profile_factory_->make_profile(complex_, edge); else return Profile(complex_, edge); } void insert_in_PQ(Edge_handle edge, Edge_data& data) { data.set_PQ_handle(heap_PQ_->push(edge)); ++current_num_edges_heap_; } void update_in_PQ(Edge_handle edge, Edge_data& data) { data.set_PQ_handle(heap_PQ_->update(edge, data.PQ_handle())); } void remove_from_PQ(Edge_handle edge, Edge_data& data) { data.set_PQ_handle(heap_PQ_->erase(edge, data.PQ_handle())); --current_num_edges_heap_; } boost::optional pop_from_PQ() { boost::optional edge = heap_PQ_->extract_top(); if (edge) get_data(*edge).reset_PQ_handle(); return edge; } private: /** * @brief Collect edges. * * Iterates over all edges of the simplicial complex and * 1) inserts them in the priority queue sorted according to the Cost policy. * 2) set the id() field of each edge */ void collect_edges() { // // Loop over all the edges in the complex in the heap // size_type size = complex_.num_edges(); DBG("Collecting edges ..."); DBGMSG("num edges :", size); edge_data_array_.reset(new Edge_data[size]); heap_PQ_.reset(new PQ(size, Compare_cost(this), Undirected_edge_id(this))); std::size_t id = 0; // xxx do a parralel for for (auto edge : complex_.edge_range()) { complex_[edge].index() = id++; Profile const& profile = create_profile(edge); Edge_data& data = get_data(edge); data.cost() = get_cost(profile); ++initial_num_edges_heap_; insert_in_PQ(edge, data); if (contraction_visitor_) contraction_visitor_->on_collected(profile, data.cost()); } DBG("Edges collected."); } bool should_stop(double lCost, const Profile &profile) const { return false; } boost::optional get_placement(const Profile& profile) const { return (*placement_policy_)(profile); } bool is_contraction_valid(Profile const& profile, Placement_type placement) const { if (!valid_contraction_policy_) return true; return (*valid_contraction_policy_)(profile, placement); } public: /** * \brief Contract edges. * * While the heap is not empty, it extracts the edge with the minimum * cost in the heap then try to contract it. * It stops when the Stop policy says so or when the number of contractions * given by 'num_max_contractions' is reached (if this number is positive). */ void contract_edges(int num_max_contractions = -1) { DBG("\n\nContract edges"); int num_contraction = 0; bool unspecified_num_contractions = (num_max_contractions == -1); // // Pops and processes each edge from the PQ // boost::optional edge; while ((edge = pop_from_PQ()) && ((num_contraction < num_max_contractions) || (unspecified_num_contractions))) { Profile const& profile = create_profile(*edge); Cost_type cost(get_data(*edge).cost()); if (contraction_visitor_) contraction_visitor_->on_selected(profile, cost, 0, 0); DBGMSG("\n\n---- Pop edge - num vertices :", complex_.num_vertices()); if (cost) { DBGMSG("sqrt(cost):", std::sqrt(*cost)); if (should_stop(*cost, profile)) { if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached(); DBG("should_stop"); break; } Placement_type placement = get_placement(profile); if (is_contraction_valid(profile, placement) && placement) { DBG("contraction_valid"); contract_edge(profile, placement); ++num_contraction; } else { DBG("contraction not valid"); if (contraction_visitor_) contraction_visitor_->on_non_valid(profile); } } else { DBG("uncomputable cost"); } } if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached(); } bool is_in_heap(Edge_handle edge) const { if (heap_PQ_->empty()) { return false; } else { return edge_data_array_[get_undirected_edge_id(edge)].is_in_PQ(); } } bool is_heap_empty() const { return heap_PQ_->empty(); } /** * @brief Returns an Edge_handle and a Placement_type. This pair consists in * the edge with the lowest cost in the heap together with its placement. * The returned value is initialized iff the heap is non-empty. */ boost::optional > top_edge() { boost::optional > res; if (!heap_PQ_->empty()) { auto edge = heap_PQ_->top(); Profile const& profile = create_profile(edge); Placement_type placement = get_placement(profile); res = std::make_pair(edge, placement); DBGMSG("top edge:", complex_[edge]); } return res; } /** * @brief Constructor with default policies. * * @details The default cost, placement, valid and visitor policies * are respectively : the edge length, the first point, the link condition */ Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex) : complex_(complex), cost_policy_(new Edge_length_cost), placement_policy_(new First_vertex_placement), valid_contraction_policy_(new Link_condition_valid_contraction), contraction_visitor_(new Contraction_visitor_()), edge_profile_factory_(0), initial_num_edges_heap_(0), current_num_edges_heap_(0) { complex_.set_visitor(this); if (contraction_visitor_) contraction_visitor_->on_started(complex); collect_edges(); } /** * @brief Constructor with customed policies. * @remark Policies destruction is handle by the class with smart pointers. */ Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex, Cost_policy_ *cost_policy, Placement_policy_ * placement_policy = new First_vertex_placement, Valid_contraction_policy_ * valid_contraction_policy = new Link_condition_valid_contraction, Contraction_visitor_* contraction_visitor = new Contraction_visitor_(), Edge_profile_factory_* edge_profile_factory = NULL) : complex_(complex), cost_policy_(cost_policy), placement_policy_(placement_policy), valid_contraction_policy_(valid_contraction_policy), contraction_visitor_(contraction_visitor), edge_profile_factory_(edge_profile_factory), initial_num_edges_heap_(0), current_num_edges_heap_(0) { complex_.set_visitor(this); if (contraction_visitor) contraction_visitor->on_started(complex); collect_edges(); } ~Skeleton_blocker_contractor() { complex_.set_visitor(0); } private: void contract_edge(const Profile& profile, Placement_type placement) { if (contraction_visitor_) contraction_visitor_->on_contracting(profile, placement); assert(complex_.contains_vertex(profile.v0_handle())); assert(complex_.contains_vertex(profile.v1_handle())); assert(placement); profile.complex().point(profile.v0_handle()) = *placement; // remark : this is not necessary since v1 will be deactivated // profile.complex().point(profile.v1_handle()) = *placement; complex_.contract_edge(profile.v0_handle(), profile.v1_handle()); assert(complex_.contains_vertex(profile.v0_handle())); assert(!complex_.contains_vertex(profile.v1_handle())); update_changed_edges(); // the visitor could do something as complex_.remove_popable_blockers(); if (contraction_visitor_) contraction_visitor_->on_contracted(profile, placement); } private: // every time the visitor's method on_changed_edge is called, it adds an // edge to changed_edges_ std::vector< Edge_handle > changed_edges_; /** * @brief we update the cost and the position in the heap of an edge that has * been changed */ inline void on_changed_edge(Vertex_handle a, Vertex_handle b) override { boost::optional ab(complex_[std::make_pair(a, b)]); assert(ab); changed_edges_.push_back(*ab); } void update_changed_edges() { // xxx do a parralel for DBG("update edges"); // sequential loop for (auto ab : changed_edges_) { // 1-get the Edge_handle corresponding to ab // 2-change the data in mEdgeArray[ab.id()] // 3-update the heap Edge_data& data = get_data(ab); Profile const& profile = create_profile(ab); data.cost() = get_cost(profile); if (data.is_in_PQ()) { update_in_PQ(ab, data); } else { insert_in_PQ(ab, data); } } changed_edges_.clear(); } private: void on_remove_edge(Vertex_handle a, Vertex_handle b) override { boost::optional ab((complex_[std::make_pair(a, b)])); assert(ab); Edge_data& lData = get_data(*ab); if (lData.is_in_PQ()) { remove_from_PQ(*ab, lData); } } private: /** * @brief Called when the edge 'ax' has been added while the edge 'bx' * is still there but will be removed on next instruction. * We assign the index of 'bx' to the edge index of 'ax' */ void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) override { boost::optional ax(complex_[std::make_pair(a, x)]); boost::optional bx(complex_[std::make_pair(b, x)]); assert(ax && bx); complex_[*ax].index() = complex_[*bx].index(); } private: /** * @brief Called when a blocker is removed. * All the edges that passes through the blocker may be edge-contractible * again and are thus reinserted in the heap. */ void on_delete_blocker(const Simplex * blocker) override { // we go for all pairs xy that belongs to the blocker // note that such pairs xy are necessarily edges of the complex // by definition of a blocker // todo uniqument utile pour la link condition // laisser a l'utilisateur ? booleen update_heap_on_removed_blocker? Simplex blocker_copy(*blocker); for (auto x = blocker_copy.begin(); x != blocker_copy.end(); ++x) { for (auto y = x; ++y != blocker_copy.end();) { auto edge_descr(complex_[std::make_pair(*x, *y)]); assert(edge_descr); Edge_data& data = get_data(*edge_descr); Profile const& profile = create_profile(*edge_descr); data.cost() = get_cost(profile); // If the edge is already in the heap // its priority has not changed. // If the edge is not present, we reinsert it // remark : we could also reinsert the edge // only if it is valid if (!data.is_in_PQ()) { insert_in_PQ(*edge_descr, data); } } } } private: std::shared_ptr cost_policy_; std::shared_ptr placement_policy_; std::shared_ptr valid_contraction_policy_; std::shared_ptr contraction_visitor_; // in case the user wants to do something special when the edge profile // are created (for instance add some info) std::shared_ptr edge_profile_factory_; Edge_data_array edge_data_array_; boost::scoped_ptr heap_PQ_; int initial_num_edges_heap_; int current_num_edges_heap_; }; } // namespace contraction } // namespace Gudhi #endif // SKELETON_BLOCKER_CONTRACTOR_H_