/* 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-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 .
*/
#ifndef GUDHI_SKELETON_BLOCKER_CONTRACTOR_H_
#define GUDHI_SKELETON_BLOCKER_CONTRACTOR_H_
#include
#include
// todo remove the queue to be independent from cgald
#include "gudhi/Contraction/CGAL_queue/Modifiable_priority_queue.h"
//#include
#include
#include
#include
#include "gudhi/Contraction/Edge_profile.h"
#include "gudhi/Contraction/policies/Cost_policy.h"
#include "gudhi/Contraction/policies/Edge_length_cost.h"
#include "gudhi/Contraction/policies/Placement_policy.h"
#include "gudhi/Contraction/policies/First_vertex_placement.h"
#include "gudhi/Contraction/policies/Valid_contraction_policy.h"
#include "gudhi/Contraction/policies/Dummy_valid_contraction.h" //xxx remove
#include "gudhi/Contraction/policies/Link_condition_valid_contraction.h" //xxx remove
#include "gudhi/Contraction/policies/Contraction_visitor.h"
#include "gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h"
#include "gudhi/Utils.h"
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();
}
// visitor that 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;
//todo explore only neighborhood with stack
void on_contracted(const Profile &profile, boost::optional< Point > placement) override{
std::list vertex_to_check;
vertex_to_check.push_front(profile.v0_handle());
while(!vertex_to_check.empty()){
Vertex_handle v = vertex_to_check.front();
vertex_to_check.pop_front();
bool blocker_popable_found=true;
while (blocker_popable_found){
blocker_popable_found = false;
for(auto block : profile.complex().blocker_range(v)){
if (profile.complex().is_popable_blocker(block)) {
for(Vertex_handle nv : *block)
if(nv!=v) vertex_to_check.push_back(nv);
profile.complex().delete_blocker(block);
blocker_popable_found = true;
break;
}
}
}
}
}
};
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.
*
* @details Basically, 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).
*
* 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).
*
* TODO expliquer la pile
*/
template<
class GeometricSimplifiableComplex,
class EdgeProfile = Edge_profile
>
class Skeleton_blocker_contractor : private skbl::Dummy_complex_visitor{
GeometricSimplifiableComplex& complex_;
public:
typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex;
typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle;
typedef typename GeometricSimplifiableComplex::Simplex_handle Simplex_handle;
typedef typename GeometricSimplifiableComplex::Simplex_handle_iterator Simplex_handle_iterator;
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 poped 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()){
// Edge_handle edge = *edge_it;
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{
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_contractionon_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(profile);
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");
}
}
}
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 = 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_handle * 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 à l'utilisateur? booleen update_heap_on_removed_blocker?
Simplex_handle 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 /* GUDHI_SKELETON_BLOCKER_CONTRACTOR_H_ */