11 #ifndef PERSISTENCE_GRAPH_H_
12 #define PERSISTENCE_GRAPH_H_
14 #include <gudhi/Internal_point.h>
18 #include <tbb/parallel_sort.h>
28 namespace persistence_diagram {
35 class Persistence_graph {
38 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
39 Persistence_graph(
const Persistence_diagram1& diag1,
const Persistence_diagram2& diag2,
double e);
41 bool on_the_u_diagonal(
int u_point_index)
const;
43 bool on_the_v_diagonal(
int v_point_index)
const;
45 int corresponding_point_in_u(
int v_point_index)
const;
47 int corresponding_point_in_v(
int u_point_index)
const;
49 double distance(
int u_point_index,
int v_point_index)
const;
53 double bottleneck_alive()
const;
55 std::vector<double> sorted_distances()
const;
57 double diameter_bound()
const;
59 Internal_point get_u_point(
int u_point_index)
const;
61 Internal_point get_v_point(
int v_point_index)
const;
64 std::vector<Internal_point> u;
65 std::vector<Internal_point> v;
69 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
70 Persistence_graph::Persistence_graph(
const Persistence_diagram1 &diag1,
71 const Persistence_diagram2 &diag2,
double e)
72 : u(), v(), b_alive(0.) {
73 std::vector<double> u_alive;
74 std::vector<double> v_alive;
75 for (
auto it = std::begin(diag1); it != std::end(diag1); ++it) {
76 if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
77 u_alive.push_back(std::get<0>(*it));
78 else if (std::get<1>(*it) - std::get<0>(*it) > e)
79 u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
81 for (
auto it = std::begin(diag2); it != std::end(diag2); ++it) {
82 if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
83 v_alive.push_back(std::get<0>(*it));
84 else if (std::get<1>(*it) - std::get<0>(*it) > e)
85 v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
87 if (u.size() < v.size())
89 std::sort(u_alive.begin(), u_alive.end());
90 std::sort(v_alive.begin(), v_alive.end());
91 if (u_alive.size() != v_alive.size()) {
92 b_alive = std::numeric_limits<double>::infinity();
94 for (
auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
95 b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
99 inline bool Persistence_graph::on_the_u_diagonal(
int u_point_index)
const {
100 return u_point_index >=
static_cast<int> (u.size());
103 inline bool Persistence_graph::on_the_v_diagonal(
int v_point_index)
const {
104 return v_point_index >=
static_cast<int> (v.size());
107 inline int Persistence_graph::corresponding_point_in_u(
int v_point_index)
const {
108 return on_the_v_diagonal(v_point_index) ?
109 v_point_index -
static_cast<int> (v.size()) : v_point_index +
static_cast<int> (u.size());
112 inline int Persistence_graph::corresponding_point_in_v(
int u_point_index)
const {
113 return on_the_u_diagonal(u_point_index) ?
114 u_point_index -
static_cast<int> (u.size()) : u_point_index +
static_cast<int> (v.size());
117 inline double Persistence_graph::distance(
int u_point_index,
int v_point_index)
const {
118 if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
120 Internal_point p_u = get_u_point(u_point_index);
121 Internal_point p_v = get_v_point(v_point_index);
122 return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
125 inline int Persistence_graph::size()
const {
126 return static_cast<int> (u.size() + v.size());
129 inline double Persistence_graph::bottleneck_alive()
const {
133 inline std::vector<double> Persistence_graph::sorted_distances()
const {
134 std::vector<double> distances;
135 distances.push_back(0.);
136 for (
int u_point_index = 0; u_point_index < size(); ++u_point_index) {
137 distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
138 for (
int v_point_index = 0; v_point_index < size(); ++v_point_index)
139 distances.push_back(distance(u_point_index, v_point_index));
142 tbb::parallel_sort(distances.begin(), distances.end());
144 std::sort(distances.begin(), distances.end());
149 inline Internal_point Persistence_graph::get_u_point(
int u_point_index)
const {
150 if (!on_the_u_diagonal(u_point_index))
151 return u.at(u_point_index);
152 Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
153 double m = (projector.x() + projector.y()) / 2.;
154 return Internal_point(m, m, u_point_index);
157 inline Internal_point Persistence_graph::get_v_point(
int v_point_index)
const {
158 if (!on_the_v_diagonal(v_point_index))
159 return v.at(v_point_index);
160 Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
161 double m = (projector.x() + projector.y()) / 2.;
162 return Internal_point(m, m, v_point_index);
165 inline double Persistence_graph::diameter_bound()
const {
167 for (
auto it = u.cbegin(); it != u.cend(); it++)
168 max = (std::max)(max, it->y());
169 for (
auto it = v.cbegin(); it != v.cend(); it++)
170 max = (std::max)(max, it->y());
178 #endif // PERSISTENCE_GRAPH_H_