Skeleton_blocker_contractor.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): David Salinas
4  *
5  * Copyright (C) 2014 Inria
6  *
7  * Modification(s):
8  * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
9  * - YYYY/MM Author: Description of the modification
10  */
11 
12 #ifndef SKELETON_BLOCKER_CONTRACTOR_H_
13 #define SKELETON_BLOCKER_CONTRACTOR_H_
14 
15 #include <gudhi/Contraction/Edge_profile.h>
16 #include <gudhi/Contraction/policies/Cost_policy.h>
17 #include <gudhi/Contraction/policies/Edge_length_cost.h>
18 #include <gudhi/Contraction/policies/Placement_policy.h>
19 #include <gudhi/Contraction/policies/First_vertex_placement.h>
20 #include <gudhi/Contraction/policies/Valid_contraction_policy.h>
21 #include <gudhi/Contraction/policies/Dummy_valid_contraction.h> // xxx remove
22 #include <gudhi/Contraction/policies/Link_condition_valid_contraction.h> // xxx remove
23 #include <gudhi/Contraction/policies/Contraction_visitor.h>
24 
25 #include <gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h>
26 #include <gudhi/Debug_utils.h>
27 
28 // todo remove the queue to be independent from cgal
29 #include <CGAL/Modifiable_priority_queue.h>
30 #include <CGAL/version.h> // for CGAL_VERSION_NR
31 
32 #include <boost/scoped_array.hpp>
33 #include <boost/scoped_ptr.hpp>
34 
35 #include <memory>
36 #include <cassert>
37 #include <list>
38 #include <utility> // for pair
39 #include <vector>
40 
41 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
42 #if CGAL_VERSION_NR < 1041101000
43 # error Skeleton_blocker_contractor is only available for CGAL >= 4.11
44 #endif
45 
46 namespace Gudhi {
47 
48 namespace contraction {
49 
50 template <class Profile>
51 Placement_policy<Profile>* make_first_vertex_placement() {
52  return new First_vertex_placement<Profile>();
53 }
54 
55 template <class Profile>
56 Valid_contraction_policy<Profile>* make_link_valid_contraction() {
57  return new Link_condition_valid_contraction<Profile>();
58 }
59 
63 template <class Profile>
65  public:
66  typedef typename Profile::Point Point;
67  typedef typename Profile::Complex Complex;
68  typedef typename Complex::Vertex_handle Vertex_handle;
69 
70  void on_contracted(const Profile &profile, boost::optional< Point > placement) override {
71  profile.complex().remove_all_popable_blockers(profile.v0_handle());
72  }
73 };
74 
75 template <class Profile>
76 Contraction_visitor<Profile>* make_remove_popable_blockers_visitor() {
78 }
79 
96 template<class GeometricSimplifiableComplex, class EdgeProfile = Edge_profile<GeometricSimplifiableComplex>>
98 typename GeometricSimplifiableComplex::Vertex_handle> {
99  GeometricSimplifiableComplex& complex_;
100 
101  public:
102  typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex;
103  typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle;
104  typedef typename GeometricSimplifiableComplex::Simplex Simplex;
105  typedef typename GeometricSimplifiableComplex::Root_vertex_handle Root_vertex_handle;
106  typedef typename GeometricSimplifiableComplex::Graph_edge Graph_edge;
107  typedef typename GeometricSimplifiableComplex::Edge_handle Edge_handle;
108  typedef typename GeometricSimplifiableComplex::Point Point;
109 
110  typedef EdgeProfile Profile;
111 
112 
117  typedef Edge_profile_factory<EdgeProfile> Edge_profile_factory_;
118 
119 
120 
121  typedef boost::optional<double> Cost_type;
122  typedef boost::optional<Point> Placement_type;
123 
124  typedef size_t size_type;
125 
127 
128  private:
129  struct Compare_id {
130  Compare_id() : algorithm_(0) { }
131 
132  Compare_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
133 
134  bool operator()(Edge_handle a, Edge_handle b) const {
135  return algorithm_->get_undirected_edge_id(a) < algorithm_->get_undirected_edge_id(b);
136  }
137 
138  Self const* algorithm_;
139  };
140 
141  struct Compare_cost {
142  Compare_cost() : algorithm_(0) { }
143 
144  Compare_cost(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
145 
146  bool operator()(Edge_handle a, Edge_handle b) const {
147  // NOTE: A cost is an optional<> value.
148  // Absent optionals are ordered first; that is, "none < T" and "T > none" for any defined T != none.
149  // In consequence, edges with undefined costs will be promoted to the top of the priority queue and popped out
150  // first.
151  return algorithm_->get_data(a).cost() < algorithm_->get_data(b).cost();
152  }
153 
154  Self const* algorithm_;
155  };
156 
157  struct Undirected_edge_id : boost::put_get_helper<size_type, Undirected_edge_id> {
158  typedef boost::readable_property_map_tag category;
159  typedef size_type value_type;
160  typedef size_type reference;
161  typedef Edge_handle key_type;
162 
163  Undirected_edge_id() : algorithm_(0) { }
164 
165  Undirected_edge_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { }
166 
167  size_type operator[](Edge_handle e) const {
168  return algorithm_->get_undirected_edge_id(e);
169  }
170 
171  Self const* algorithm_;
172  };
173 
174 #if CGAL_VERSION_NR < 1050500000
175  typedef CGAL::Modifiable_priority_queue<Edge_handle, Compare_cost, Undirected_edge_id> PQ;
176 #else
177  typedef CGAL::Modifiable_priority_queue<Edge_handle, Compare_cost, Undirected_edge_id, CGAL::CGAL_BOOST_PAIRING_HEAP> PQ;
178 #endif
179 
180  typedef bool pq_handle;
181 
182 
183  // An Edge_data is associated with EVERY edge in the complex (collapsible or not).
184  // It relates the edge with the PQ-handle needed to update the priority queue
185  // It also relates the edge with a policy-based cache
186 
187  class Edge_data {
188  public:
189  Edge_data() : PQHandle_(), cost_() { }
190 
191  Cost_type const& cost() const {
192  return cost_;
193  }
194 
195  Cost_type & cost() {
196  return cost_;
197  }
198 
199  pq_handle PQ_handle() const {
200  return PQHandle_;
201  }
202 
203  bool is_in_PQ() const {
204  return PQHandle_ != false;
205  }
206 
207  void set_PQ_handle(pq_handle h) {
208  PQHandle_ = h;
209  }
210 
211  void reset_PQ_handle() {
212  PQHandle_ = false;
213  }
214 
215  private:
216  pq_handle PQHandle_;
217  Cost_type cost_;
218  };
219  typedef Edge_data* Edge_data_ptr;
220  typedef boost::scoped_array<Edge_data> Edge_data_array;
221 
222  int get_undirected_edge_id(Edge_handle edge) const {
223  return complex_[edge].index();
224  }
225 
226  const Edge_data& get_data(Edge_handle edge) const {
227  return edge_data_array_[get_undirected_edge_id(edge)];
228  }
229 
230  Edge_data& get_data(Edge_handle edge) {
231  return edge_data_array_[get_undirected_edge_id(edge)];
232  }
233 
234  Cost_type get_cost(const Profile & profile) const {
235  return (*cost_policy_)(profile, get_placement(profile));
236  }
237 
238  Profile create_profile(Edge_handle edge) const {
239  if (edge_profile_factory_)
240  return edge_profile_factory_->make_profile(complex_, edge);
241  else
242  return Profile(complex_, edge);
243  }
244 
245  void insert_in_PQ(Edge_handle edge, Edge_data& data) {
246  heap_PQ_->push(edge);
247  data.set_PQ_handle(true);
248  ++current_num_edges_heap_;
249  }
250 
251  void update_in_PQ(Edge_handle edge, Edge_data& data) {
252 #if CGAL_VERSION_NR < 1050500000
253  data.set_PQ_handle(heap_PQ_->update(edge, data.PQ_handle()));
254 #else
255  heap_PQ_->update(edge);
256 #endif
257  }
258 
259  void remove_from_PQ(Edge_handle edge, Edge_data& data) {
260  heap_PQ_->erase(edge);
261  data.set_PQ_handle(false);
262  --current_num_edges_heap_;
263  }
264 
265  boost::optional<Edge_handle> pop_from_PQ() {
266  boost::optional<Edge_handle> edge = heap_PQ_->extract_top();
267  if (edge)
268  get_data(*edge).reset_PQ_handle();
269  return edge;
270  }
271 
272  private:
280  void collect_edges() {
281  //
282  // Loop over all the edges in the complex in the heap
283  //
284  size_type size = complex_.num_edges();
285  DBG("Collecting edges ...");
286  DBGMSG("num edges :", size);
287 
288  edge_data_array_.reset(new Edge_data[size]);
289 
290  heap_PQ_.reset(new PQ(size, Compare_cost(this), Undirected_edge_id(this)));
291 
292  std::size_t id = 0;
293 
294  // xxx do a parralel for
295  for (auto edge : complex_.edge_range()) {
296  complex_[edge].index() = id++;
297  Profile const& profile = create_profile(edge);
298  Edge_data& data = get_data(edge);
299  data.cost() = get_cost(profile);
300  ++initial_num_edges_heap_;
301  insert_in_PQ(edge, data);
302  if (contraction_visitor_) contraction_visitor_->on_collected(profile, data.cost());
303  }
304 
305  DBG("Edges collected.");
306  }
307 
308  bool should_stop(double lCost, const Profile &profile) const {
309  return false;
310  }
311 
312  boost::optional<Point> get_placement(const Profile& profile) const {
313  return (*placement_policy_)(profile);
314  }
315 
316  bool is_contraction_valid(Profile const& profile, Placement_type placement) const {
317  if (!valid_contraction_policy_) return true;
318  return (*valid_contraction_policy_)(profile, placement);
319  }
320 
321 
322  public:
331  void contract_edges(int num_max_contractions = -1) {
332  DBG("\n\nContract edges");
333  int num_contraction = 0;
334 
335  bool unspecified_num_contractions = (num_max_contractions == -1);
336  //
337  // Pops and processes each edge from the PQ
338  //
339  boost::optional<Edge_handle> edge;
340  while ((edge = pop_from_PQ()) && ((num_contraction < num_max_contractions) || (unspecified_num_contractions))) {
341  Profile const& profile = create_profile(*edge);
342  Cost_type cost(get_data(*edge).cost());
343  if (contraction_visitor_) contraction_visitor_->on_selected(profile, cost, 0, 0);
344 
345  DBGMSG("\n\n---- Pop edge - num vertices :", complex_.num_vertices());
346 
347  if (cost) {
348  DBGMSG("sqrt(cost):", std::sqrt(*cost));
349  if (should_stop(*cost, profile)) {
350  if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached();
351  DBG("should_stop");
352  break;
353  }
354  Placement_type placement = get_placement(profile);
355  if (is_contraction_valid(profile, placement) && placement) {
356  DBG("contraction_valid");
357  contract_edge(profile, placement);
358  ++num_contraction;
359  } else {
360  DBG("contraction not valid");
361  if (contraction_visitor_) contraction_visitor_->on_non_valid(profile);
362  }
363  } else {
364  DBG("uncomputable cost");
365  }
366  }
367  if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached();
368  }
369 
370  bool is_in_heap(Edge_handle edge) const {
371  if (heap_PQ_->empty()) {
372  return false;
373  } else {
374  return edge_data_array_[get_undirected_edge_id(edge)].is_in_PQ();
375  }
376  }
377 
378  bool is_heap_empty() const {
379  return heap_PQ_->empty();
380  }
381 
387  boost::optional<std::pair<Edge_handle, Placement_type > > top_edge() {
388  boost::optional<std::pair<Edge_handle, Placement_type > > res;
389 
390  if (!heap_PQ_->empty()) {
391  auto edge = heap_PQ_->top();
392  Profile const& profile = create_profile(edge);
393  Placement_type placement = get_placement(profile);
394  res = std::make_pair(edge, placement);
395  DBGMSG("top edge:", complex_[edge]);
396  }
397  return res;
398  }
399 
406  Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex)
407  : complex_(complex),
408  cost_policy_(new Edge_length_cost<Profile>),
409  placement_policy_(new First_vertex_placement<Profile>),
410  valid_contraction_policy_(new Link_condition_valid_contraction<Profile>),
411  contraction_visitor_(new Contraction_visitor_()),
412  edge_profile_factory_(0),
413  initial_num_edges_heap_(0),
414  current_num_edges_heap_(0) {
415  complex_.set_visitor(this);
416  if (contraction_visitor_) contraction_visitor_->on_started(complex);
417  collect_edges();
418  }
419 
424  Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex,
425  Cost_policy_ *cost_policy,
426  Placement_policy_ * placement_policy = new First_vertex_placement<Profile>,
427  Valid_contraction_policy_ * valid_contraction_policy =
429  Contraction_visitor_* contraction_visitor = new Contraction_visitor_(),
430  Edge_profile_factory_* edge_profile_factory = NULL) :
431  complex_(complex),
432  cost_policy_(cost_policy),
433  placement_policy_(placement_policy),
434  valid_contraction_policy_(valid_contraction_policy),
435  contraction_visitor_(contraction_visitor),
436  edge_profile_factory_(edge_profile_factory),
437  initial_num_edges_heap_(0),
438  current_num_edges_heap_(0) {
439  complex_.set_visitor(this);
440  if (contraction_visitor) contraction_visitor->on_started(complex);
441  collect_edges();
442  }
443 
445  complex_.set_visitor(0);
446  }
447 
448  private:
449  void contract_edge(const Profile& profile, Placement_type placement) {
450  if (contraction_visitor_) contraction_visitor_->on_contracting(profile, placement);
451 
452  assert(complex_.contains_vertex(profile.v0_handle()));
453  assert(complex_.contains_vertex(profile.v1_handle()));
454  assert(placement);
455 
456  profile.complex().point(profile.v0_handle()) = *placement;
457 
458  // remark : this is not necessary since v1 will be deactivated
459  // profile.complex().point(profile.v1_handle()) = *placement;
460 
461  complex_.contract_edge(profile.v0_handle(), profile.v1_handle());
462 
463  assert(complex_.contains_vertex(profile.v0_handle()));
464  assert(!complex_.contains_vertex(profile.v1_handle()));
465 
466  update_changed_edges();
467 
468  // the visitor could do something as complex_.remove_popable_blockers();
469  if (contraction_visitor_) contraction_visitor_->on_contracted(profile, placement);
470  }
471 
472  private:
473  // every time the visitor's method on_changed_edge is called, it adds an
474  // edge to changed_edges_
475  std::vector< Edge_handle > changed_edges_;
476 
481  inline void on_changed_edge(Vertex_handle a, Vertex_handle b) override {
482  boost::optional<Edge_handle> ab(complex_[std::make_pair(a, b)]);
483  assert(ab);
484  changed_edges_.push_back(*ab);
485  }
486 
487  void update_changed_edges() {
488  // xxx do a parralel for
489  DBG("update edges");
490 
491  // sequential loop
492  for (auto ab : changed_edges_) {
493  // 1-get the Edge_handle corresponding to ab
494  // 2-change the data in mEdgeArray[ab.id()]
495  // 3-update the heap
496  Edge_data& data = get_data(ab);
497  Profile const& profile = create_profile(ab);
498  data.cost() = get_cost(profile);
499  if (data.is_in_PQ()) {
500  update_in_PQ(ab, data);
501  } else {
502  insert_in_PQ(ab, data);
503  }
504  }
505  changed_edges_.clear();
506  }
507 
508 
509  private:
510  void on_remove_edge(Vertex_handle a, Vertex_handle b) override {
511  boost::optional<Edge_handle> ab((complex_[std::make_pair(a, b)]));
512  assert(ab);
513  Edge_data& lData = get_data(*ab);
514  if (lData.is_in_PQ()) {
515  remove_from_PQ(*ab, lData);
516  }
517  }
518 
519  private:
525  void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) override {
526  boost::optional<Edge_handle> ax(complex_[std::make_pair(a, x)]);
527  boost::optional<Edge_handle> bx(complex_[std::make_pair(b, x)]);
528  assert(ax && bx);
529  complex_[*ax].index() = complex_[*bx].index();
530  }
531 
532  private:
538  void on_delete_blocker(const Simplex * blocker) override {
539  // we go for all pairs xy that belongs to the blocker
540  // note that such pairs xy are necessarily edges of the complex
541  // by definition of a blocker
542 
543  // todo uniqument utile pour la link condition
544  // laisser a l'utilisateur ? booleen update_heap_on_removed_blocker?
545  Simplex blocker_copy(*blocker);
546  for (auto x = blocker_copy.begin(); x != blocker_copy.end(); ++x) {
547  for (auto y = x; ++y != blocker_copy.end();) {
548  auto edge_descr(complex_[std::make_pair(*x, *y)]);
549  assert(edge_descr);
550  Edge_data& data = get_data(*edge_descr);
551  Profile const& profile = create_profile(*edge_descr);
552  data.cost() = get_cost(profile);
553 
554  // If the edge is already in the heap
555  // its priority has not changed.
556  // If the edge is not present, we reinsert it
557  // remark : we could also reinsert the edge
558  // only if it is valid
559  if (!data.is_in_PQ()) {
560  insert_in_PQ(*edge_descr, data);
561  }
562  }
563  }
564  }
565 
566 
567  private:
568  std::shared_ptr<Cost_policy_> cost_policy_;
569  std::shared_ptr<Placement_policy_> placement_policy_;
570  std::shared_ptr<Valid_contraction_policy_> valid_contraction_policy_;
571  std::shared_ptr<Contraction_visitor_> contraction_visitor_;
572 
573  // in case the user wants to do something special when the edge profile
574  // are created (for instance add some info)
575  std::shared_ptr<Edge_profile_factory_> edge_profile_factory_;
576  Edge_data_array edge_data_array_;
577 
578  boost::scoped_ptr<PQ> heap_PQ_;
579  int initial_num_edges_heap_;
580  int current_num_edges_heap_;
581 };
582 
583 } // namespace contraction
584 
585 } // namespace Gudhi
586 
587 #endif // SKELETON_BLOCKER_CONTRACTOR_H_
Visitor to remove popable blockers after an edge contraction.
Definition: Skeleton_blocker_contractor.h:64
Interface for a visitor of the edge contraction process.
Definition: Contraction_visitor.h:27
Policy to specify the cost of contracting an edge.
Definition: Cost_policy.h:25
return a cost corresponding to the squared length of the edge
Definition: Edge_length_cost.h:24
Places the contracted point onto the first point of the edge.
Definition: First_vertex_placement.h:24
Policy to specify where the merged point had to be placed after an edge contraction.
Definition: Placement_policy.h:25
Class that allows to contract iteratively edges of a simplicial complex.
Definition: Skeleton_blocker_contractor.h:98
void contract_edges(int num_max_contractions=-1)
Contract edges.
Definition: Skeleton_blocker_contractor.h:331
boost::optional< std::pair< Edge_handle, Placement_type > > top_edge()
Returns an Edge_handle and a Placement_type. This pair consists in the edge with the lowest cost in t...
Definition: Skeleton_blocker_contractor.h:387
Skeleton_blocker_contractor(GeometricSimplifiableComplex &complex, Cost_policy_ *cost_policy, Placement_policy_ *placement_policy=new First_vertex_placement< Profile >, Valid_contraction_policy_ *valid_contraction_policy=new Link_condition_valid_contraction< Profile >, Contraction_visitor_ *contraction_visitor=new Contraction_visitor_(), Edge_profile_factory_ *edge_profile_factory=NULL)
Constructor with customed policies.
Definition: Skeleton_blocker_contractor.h:424
Skeleton_blocker_contractor(GeometricSimplifiableComplex &complex)
Constructor with default policies.
Definition: Skeleton_blocker_contractor.h:406
Policy to specify if an edge contraction is valid or not.
Definition: Valid_contraction_policy.h:23
A dummy visitor of a simplicial complex that does nothing.
Definition: Skeleton_blocker_complex_visitor.h:65
Abstract simplex used in Skeleton blockers data-structure.
Definition: Skeleton_blocker_simplex.h:38
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Aug 16 2022 14:01:50 for GUDHIdev by Doxygen 1.9.1