Neighbors_finder.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author: Francois Godi
6  *
7  * Copyright (C) 2015 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef NEIGHBORS_FINDER_H_
24 #define NEIGHBORS_FINDER_H_
25 
26 // Inclusion order is important for CGAL patch
27 #include <CGAL/Kd_tree.h>
28 #include <CGAL/Search_traits.h>
29 
30 #include <gudhi/Persistence_graph.h>
31 #include <gudhi/Internal_point.h>
32 
33 #include <unordered_set>
34 #include <vector>
35 #include <algorithm> // for std::max
36 
37 namespace Gudhi {
38 
39 namespace persistence_diagram {
40 
43 struct Square_query {
44  typedef CGAL::Dimension_tag<2> D;
45  typedef Internal_point Point_d;
46  typedef double FT;
47  bool contains(Point_d p) const {
48  return std::max(std::abs(p.x()-c.x()), std::abs(p.y()-c.y())) <= size;
49  }
50  bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT, D> const&r) const {
51  return
52  r.max_coord(0) >= c.x() - size &&
53  r.min_coord(0) <= c.x() + size &&
54  r.max_coord(1) >= c.y() - size &&
55  r.min_coord(1) <= c.y() + size;
56  }
57  bool outer_range_contains(CGAL::Kd_tree_rectangle<FT, D> const&r) const {
58  return
59  r.min_coord(0) >= c.x() - size &&
60  r.max_coord(0) <= c.x() + size &&
61  r.min_coord(1) >= c.y() - size &&
62  r.max_coord(1) <= c.y() + size;
63  }
64  Point_d c;
65  FT size;
66 };
67 
76 class Neighbors_finder {
77  typedef CGAL::Dimension_tag<2> D;
78  typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits;
79  typedef CGAL::Kd_tree<Traits> Kd_tree;
80 
81  public:
83  Neighbors_finder(const Persistence_graph& g, double r);
85  void add(int v_point_index);
88  int pull_near(int u_point_index);
90  std::vector<int> pull_all_near(int u_point_index);
91 
92  private:
93  const Persistence_graph& g;
94  const double r;
95  Kd_tree kd_t;
96  std::unordered_set<int> projections_f;
97 };
98 
107 class Layered_neighbors_finder {
108  public:
110  Layered_neighbors_finder(const Persistence_graph& g, double r);
112  void add(int v_point_index, int vlayer);
115  int pull_near(int u_point_index, int vlayer);
117  int vlayers_number() const;
118 
119  private:
120  const Persistence_graph& g;
121  const double r;
122  std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder;
123 };
124 
125 inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) :
126  g(g), r(r), kd_t(), projections_f() { }
127 
128 inline void Neighbors_finder::add(int v_point_index) {
129  if (g.on_the_v_diagonal(v_point_index))
130  projections_f.emplace(v_point_index);
131  else
132  kd_t.insert(g.get_v_point(v_point_index));
133 }
134 
135 inline int Neighbors_finder::pull_near(int u_point_index) {
136  int tmp;
137  int c = g.corresponding_point_in_v(u_point_index);
138  if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) {
139  // Any pair of projection is at distance 0
140  tmp = *projections_f.cbegin();
141  projections_f.erase(tmp);
142  } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) {
143  // Is the query point near to its projection ?
144  tmp = c;
145  projections_f.erase(tmp);
146  } else {
147  // Is the query point near to a V point in the plane ?
148  Internal_point u_point = g.get_u_point(u_point_index);
149  auto neighbor = kd_t.search_any_point(Square_query{u_point, r});
150  if (!neighbor)
151  return null_point_index();
152  tmp = neighbor->point_index;
153  auto point = g.get_v_point(tmp);
154  int idx = point.point_index;
155  kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;});
156  }
157  return tmp;
158 }
159 
160 inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) {
161  std::vector<int> all_pull;
162  int last_pull = pull_near(u_point_index);
163  while (last_pull != null_point_index()) {
164  all_pull.emplace_back(last_pull);
165  last_pull = pull_near(u_point_index);
166  }
167  return all_pull;
168 }
169 
170 inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) :
171  g(g), r(r), neighbors_finder() { }
172 
173 inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
174  for (int l = neighbors_finder.size(); l <= vlayer; l++)
175  neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r)));
176  neighbors_finder.at(vlayer)->add(v_point_index);
177 }
178 
179 inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
180  if (static_cast<int> (neighbors_finder.size()) <= vlayer)
181  return null_point_index();
182  return neighbors_finder.at(vlayer)->pull_near(u_point_index);
183 }
184 
185 inline int Layered_neighbors_finder::vlayers_number() const {
186  return static_cast<int> (neighbors_finder.size());
187 }
188 
189 } // namespace persistence_diagram
190 
191 } // namespace Gudhi
192 
193 #endif // NEIGHBORS_FINDER_H_
Definition: SimplicialComplexForAlpha.h:26
GUDHI  Version 2.2.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Tue Jul 17 2018 12:56:28 for GUDHI by Doxygen 1.8.13