1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
* See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
* Author(s): David Salinas
*
* Copyright (C) 2014 Inria
*
* Modification(s):
* - YYYY/MM Author: Description of the modification
*/
#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"
/**
* Iteratively tries to anticollapse smallest edge non added so far.
* If its link is contractible then no topological change and else possible topological change.
*
* 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_without_blockers(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);
(void)contractor;
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 // UTILS_CRITICAL_POINTS_H_
|