PGROUTING  3.2
pgr_alphaShape.cpp
Go to the documentation of this file.
1 /*PGR-GNU*****************************************************************
2 file: pgr_alphaShape.cpp
3 
4 Copyright (c) 2018 pgRouting developers
5 Mail: project@pgrouting.org
6 
7 Copyright (c) 2018 Celia Virginia Vergara Castillo
8 vicky_vergara@hotmail.com
9 
10 ------
11 
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version.
16 
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU General Public License for more details.
21 
22 You should have received a copy of the GNU General Public License
23 along with this program; if not, write to the Free Software
24 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 
26  ********************************************************************PGR-GNU*/
27 
29 
31 #include <boost/graph/dijkstra_shortest_paths.hpp>
32 #include <boost/graph/filtered_graph.hpp>
33 
34 #include <boost/geometry/algorithms/num_points.hpp>
35 #include <boost/geometry/algorithms/append.hpp>
36 #include <boost/geometry/algorithms/area.hpp>
37 #if Boost_VERSION_MACRO >= 107500
38 # include <boost/geometry/strategies/strategies.hpp>
39 #else
40 # include <boost/geometry/strategies/agnostic/point_in_point.hpp>
41 #endif
42 #include <boost/geometry/algorithms/correct.hpp>
43 #include <boost/geometry/algorithms/distance.hpp>
44 
45 
46 #include <sstream>
47 #include <set>
48 #include <vector>
49 #include <algorithm>
50 #include <utility>
51 #include <map>
53 
54 namespace bg = boost::geometry;
55 
56 namespace pgrouting {
57 namespace alphashape {
58 
59 namespace {
60 
61 
62 /*
63  * Determinant of this matrix
64  * r00, r01
65  * r10, r11
66  */
67 double
68 det(double r00, double r01, double r10, double r11) {
69  return r00 * r11 - r01 * r10;
70 }
71 
72 Bpoint
73 circumcenter(const Bpoint a, const Bpoint b, const Bpoint c) {
74  auto cx = c.x();
75  auto cy = c.y();
76  auto ax = a.x() - cx;
77  auto ay = a.y() - cy;
78  auto bx = b.x() - cx;
79  auto by = b.y() - cy;
80 
81  auto denom = 2 * det(ax, ay, bx, by);
82  /*
83  * If denom == 0 then points are colinear
84  */
85  pgassert(denom != 0);
86 
87  auto numx = det(ay, ax * ax + ay * ay, by, bx * bx + by * by);
88  auto numy = det(ax, ax * ax + ay * ay, bx, bx * bx + by * by);
89 
90  return Bpoint {cx - numx / denom, cy + numy / denom};
91 }
92 
93 template <typename B_G, typename V>
94 std::vector<V>
95 get_predecessors(V source, V target, const B_G &subg) {
96  pgassert(boost::num_vertices(subg));
97  std::vector<V> predecessors(boost::num_vertices(subg));
98  std::vector<double> distances(num_vertices(subg));
99  pgassert(predecessors.size() == boost::num_vertices(subg));
100  pgassert(distances.size() == boost::num_vertices(subg));
101 
102  /* abort in case of an interruption occurs (e.g. the query is being cancelled) */
103  CHECK_FOR_INTERRUPTS();
104 
105  try {
106  boost::dijkstra_shortest_paths(subg, source,
107  boost::predecessor_map(&predecessors[0])
108  .weight_map(get(&Basic_edge::cost, subg))
109  .distance_map(&distances[0])
110  .visitor(visitors::dijkstra_one_goal_visitor<V>(target)));
111  } catch(found_goals &) {
112  } catch (boost::exception const& ex) {
113  (void)ex;
114  throw;
115  } catch (std::exception &e) {
116  (void)e;
117  throw;
118  } catch (...) {
119  throw;
120  }
121  return predecessors;
122 }
123 
124 
125 template <typename B_G, typename V>
126 Bpoly
127 get_polygon(V source, V target, const std::vector<V> & predecessors, const B_G &graph, std::set<E> &edges_used) {
128  Bpoly polygon;
129  /*
130  * There is no path -> returning empty polygon
131  */
132  if (target == predecessors[target]) {
133  pgassert(bg::num_points(polygon) == 0);
134  return polygon;
135  }
136 
137  /*
138  * the last stop is the target
139  */
140  bg::append(polygon.outer(), graph[source].point);
141 
142  /*
143  * get the path
144  */
145  while (target != source || target != predecessors[target]) {
146  bg::append(polygon.outer(), graph[target].point);
147  edges_used.insert(boost::edge(predecessors[target], target, graph).first);
148  target = predecessors[target];
149  }
150  bg::correct(polygon);
151  return polygon;
152 }
153 
154 typedef std::pair<Triangle, double> MyPairType;
156  bool operator()(const MyPairType& lhs, const MyPairType& rhs) const {
157  return lhs.second > rhs.second;
158  }
159 };
160 
161 } // namespace
162 
163 
164 /*
165  * Constructor
166  */
167 Pgr_alphaShape::Pgr_alphaShape(const std::vector<Pgr_edge_xy_t> &edges) :
168 graph(UNDIRECTED) {
169  graph.insert_edges(edges);
170  make_triangles();
171 }
172 
173 void
175  /*
176  * triangle sides:
177  * a_r, b, c edge descriptor
178  */
179  BGL_FORALL_EDGES(c, graph.graph, BG) {
180  /*
181  * triangle vertices:
182  * u, v, w vertex descriptor
183  */
184  auto u = graph.source(c);
185  auto v = graph.target(c);
186 
187  std::vector<Triangle> adjacent_to_side;
188 
189  size_t i = 0;
190  BGL_FORALL_OUTEDGES(u, b, graph.graph, BG) {
191  ++i;
192  auto w = graph.adjacent(u, b);
193  if (w == v) {
194  pgassert(b == c);
195  continue;
196  }
197 
198  auto a_r = boost::edge(v, w, graph.graph);
199  if (!a_r.second) continue;
200 
201  Triangle face{{a_r.first, b, c}};
202  adjacent_to_side.push_back(face);
203  }
204 
205  /*
206  * All vertices must have at least 2 edges
207  * So cycle above must have passed at least twice
208  */
209  pgassert(i > 1);
210 
211  if (adjacent_to_side.size() == 2) {
212  m_adjacent_triangles[adjacent_to_side[0]].insert(adjacent_to_side[1]);
213  m_adjacent_triangles[adjacent_to_side[1]].insert(adjacent_to_side[0]);
214  } else {
215  if (m_adjacent_triangles.find(adjacent_to_side[0]) == m_adjacent_triangles.end()) {
216  m_adjacent_triangles[adjacent_to_side[0]].clear();
217  }
218  }
219  }
220 }
221 
222 
223 /*
224  * Radius of triangle's circumcenter
225  */
226 double
228  std::vector<E> edges(t.begin(), t.end());
229  auto a = graph.source(edges[0]);
230  auto b = graph.target(edges[0]);
231  auto c = graph.source(edges[1]);
232  c = (c == a || c == b)? graph.target(edges[1]) : c;
233  pgassert(a != b && a != c && b!= c);
234 
235  auto center = circumcenter(graph[a].point, graph[b].point, graph[c].point);
236  return bg::distance(center, graph[a].point);
237 }
238 
239 /*
240  * The whole traingle face belongs to the shape?
241  */
242 bool
243 Pgr_alphaShape::faceBelongs(const Triangle t, double alpha) const {
244  return radius(t) <= alpha;
245 }
246 
247 std::vector<Bpoly>
249  std::map<Triangle, double> border_triangles;
250  std::map<Triangle, double> inner_triangles;
251 
252  size_t i(0);
253  for (const auto &t : m_adjacent_triangles) {
254  if (t.second.size() == 2) {
255  border_triangles[t.first] = radius(t.first);
256  } else {
257  inner_triangles[t.first] = radius(t.first);
258  }
259  ++i;
260  }
261  pgassert(border_triangles.size() + inner_triangles.size() == i);
262 
263  auto max_border_triangle = *std::min_element(border_triangles.begin(), border_triangles.end(), CompareRadius());
264  auto max_inner_triangle = *std::min_element(inner_triangles.begin(), inner_triangles.end(), CompareRadius());
265 
266  double max_border_radius = max_border_triangle.second;
267  double max_inner_radius = max_inner_triangle.second;
268 
269 
270  auto count = border_triangles.size() + inner_triangles.size();
271  while (max_border_radius >= max_inner_radius) {
272  auto max_border_triangle = *std::min_element(border_triangles.begin(), border_triangles.end(), CompareRadius());
273  auto max_inner_triangle = *std::min_element(inner_triangles.begin(), inner_triangles.end(), CompareRadius());
274  /*
275  * Removing largest border triangle
276  */
277  border_triangles.erase(max_border_triangle.first);
278 
279  /*
280  * Adjacent triangles of a border triangle
281  * - are no longer inner triangles
282  * - are now border triangles
283  */
284  for (const auto &t : m_adjacent_triangles.at(max_border_triangle.first)) {
285  if (inner_triangles.find(t) != inner_triangles.end()) {
286  inner_triangles.erase(t);
287  border_triangles[t] = radius(t);
288  }
289  }
290 
291  auto new_max_border_triangle = *std::min_element(
292  border_triangles.begin(), border_triangles.end(), CompareRadius());
293  auto new_max_inner_triangle = *std::min_element(
294  inner_triangles.begin(), inner_triangles.end(), CompareRadius());
295 
296  if (new_max_border_triangle.second < new_max_inner_triangle.second) {
297  /*
298  * Roll back and exit loop
299  */
300  for (const auto &t : m_adjacent_triangles.at(max_border_triangle.first)) {
301  border_triangles.erase(t);
302  inner_triangles[t] = radius(t);
303  }
304  border_triangles[max_border_triangle.first] = max_border_triangle.second;
305  break;
306  }
307  max_border_radius = new_max_border_triangle.second;
308  max_inner_radius = new_max_inner_triangle.second;
309  pgassert(count > (border_triangles.size() + inner_triangles.size()));
310  count = border_triangles.size() + inner_triangles.size();;
311  }
312 
313  pgassert(max_border_radius > 0);
314  log << "Using Alpha:" << max_border_radius * max_border_radius;
315  return this->operator()(max_border_radius);
316 }
317 
318 void
320  const Triangle face,
321  std::set<Triangle> &used,
322  std::set<E> &border_edges,
323  double alpha) const {
324  /*
325  * Do nothing when the face does not belong to the polygon of the alphashape
326  */
327  if (!faceBelongs(face, alpha)) return;
328 
329  /*
330  * Do nothing when the face has being processed before
331  */
332  auto original = used.size();
333  used.insert(face);
334  if (original == used.size()) return;
335 
336  std::set<E> common_sides;
337  for (const auto &adj_t : m_adjacent_triangles.at(face)) {
338  if (!faceBelongs(adj_t, alpha)) {
339  /*
340  * Adjacent face is not on shape
341  * then side is on the border
342  */
343  std::set_intersection(
344  face.begin(), face.end(),
345  adj_t.begin(), adj_t.end(),
346  std::inserter(border_edges, border_edges.end()));
347  }
348  std::set_intersection(
349  face.begin(), face.end(),
350  adj_t.begin(), adj_t.end(),
351  std::inserter(common_sides, common_sides.end()));
352  recursive_build(adj_t, used, border_edges, alpha);
353  }
354 
355  if (m_adjacent_triangles.at(face).size() == 2) {
356  /*
357  * Side without adjacent triangle (in convex hull)
358  * is part of the border
359  */
360  std::set_difference(
361  face.begin(), face.end(),
362  common_sides.begin(), common_sides.end(),
363  std::inserter(border_edges, border_edges.end()));
364  }
365 }
366 
367 std::vector<Bpoly>
368 Pgr_alphaShape::operator() (double alpha) const {
369  std::vector<Bpoly> shape;
370 
371  if (alpha <= 0) return build_best_alpha();
372 
373  std::set<Triangle> used;
374  using Subgraph = boost::filtered_graph<BG, EdgesFilter, boost::keep_all>;
375 
376  for (const auto &t : m_adjacent_triangles) {
377  EdgesFilter border_edges;
378  /*
379  * Recurse through the triangles to get the border sides
380  */
381  recursive_build(t.first, used, border_edges.edges, alpha);
382  /*
383  * triangle was already processed
384  * or is not part of the shape
385  */
386  if (border_edges.edges.empty()) continue;
387 
388  std::vector<Bpoly> polys;
389  Bpoly polygon;
390  double area = 0;
391  while (!border_edges.edges.empty()) {
392  auto first_edge = *border_edges.edges.begin();
393  border_edges.edges.erase(first_edge);
394 
395  Subgraph subg(graph.graph, border_edges, {});
396 
397  auto source = boost::source(first_edge, subg);
398  auto target = boost::target(first_edge, subg);
399 
400 
401  auto predecessors = get_predecessors(source, target, subg);
402  std::set<E> edges_used;
403  auto poly = get_polygon(source, target, predecessors, subg, edges_used);
404 
405  for (const auto &e : edges_used) {
406  border_edges.edges.erase(e);
407  }
408 
409  pgassert(bg::num_points(poly) >= 3 || bg::num_points(poly) == 0);
410 
411  if (bg::num_points(poly) == 0) continue;
412 
413  /*
414  * A polygon must have at least 3 vertices
415  */
416  if (bg::num_points(poly) >= 3) {
417  if (area == 0) {
418  /*
419  * consider the first polygon as the outer polygon
420  */
421  polygon = poly;
422  area = bg::area(poly);
423  } else {
424  auto new_area = bg::area(poly);
425  if (new_area > area) {
426  /*
427  * if new poly is larger than the current outer polygon
428  * - save current outer polygon as inner ring
429  * - save new poly as outer polygon
430  */
431  polys.push_back(polygon);
432  area = new_area;
433  polygon = poly;
434  } else {
435  /*
436  * New poly is inner ring
437  */
438  polys.push_back(poly);
439  }
440  }
441  }
442  }
443 
444  /*
445  * Add inner rings to the polygon
446  */
447  if (!polys.empty()) {
448  polygon.inners().resize(polys.size());
449  size_t i(0);
450  for (const auto &inner_ring : polys) {
451  bg::assign(polygon.inners()[i++], inner_ring);
452  }
453  }
454 
455  /*
456  * Add the polygon to the alpha shape
457  */
458  shape.push_back(polygon);
459  }
460 
461  return shape;
462 }
463 
464 
465 
466 
467 std::ostream&
468 operator<<(std::ostream& os, const Pgr_alphaShape &d) {
469  os << d.graph;
470 
471  return os;
472 }
473 
474 } // namespace alphashape
475 } // namespace pgrouting
interruption.h
pgrouting::graph::Pgr_base_graph::graph
G graph
The graph.
Definition: pgr_base_graph.hpp:260
pgrouting::Bpoly
bg::model::polygon< Bpoint > Bpoly
Definition: bline.hpp:46
pgrouting::Bpoint
bg::model::d2::point_xy< double > Bpoint
Definition: bline.hpp:40
pgrouting::alphashape::Pgr_alphaShape::make_triangles
void make_triangles()
Definition: pgr_alphaShape.cpp:174
pgrouting::alphashape::Pgr_alphaShape::faceBelongs
bool faceBelongs(const Triangle face, double alpha) const
Definition: pgr_alphaShape.cpp:243
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::MyPairType
std::pair< Triangle, double > MyPairType
Definition: pgr_alphaShape.cpp:154
dijkstra_one_goal_visitor.hpp
pgrouting::graph::Pgr_base_graph::target
V target(E e_idx) const
Definition: pgr_base_graph.hpp:561
pgrouting::alphashape::Triangle
std::set< E > Triangle
Definition: pgr_alphaShape.h:60
pgrouting::Pgr_messages::log
std::ostringstream log
Stores the hint information.
Definition: pgr_messages.h:81
pgrouting::alphashape::Pgr_alphaShape::Pgr_alphaShape
Pgr_alphaShape()=delete
pgrouting::visitors::dijkstra_one_goal_visitor
Definition: dijkstra_one_goal_visitor.hpp:39
pgassert
#define pgassert(expr)
Uses the standard assert syntax.
Definition: pgr_assert.h:94
pgrouting::alphashape::operator<<
std::ostream & operator<<(std::ostream &os, const Pgr_alphaShape &d)
Definition: pgr_alphaShape.cpp:468
UNDIRECTED
@ UNDIRECTED
Definition: graph_enum.h:30
pgrouting::graph::Pgr_base_graph::insert_edges
void insert_edges(const T *edges, size_t count)
Inserts count edges of type T into the graph.
Definition: pgr_base_graph.hpp:357
pgrouting::alphashape::Pgr_alphaShape::graph
G graph
sides graph
Definition: pgr_alphaShape.h:93
pgrouting::alphashape::Pgr_alphaShape::build_best_alpha
std::vector< Bpoly > build_best_alpha() const
Definition: pgr_alphaShape.cpp:248
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::CompareRadius
Definition: pgr_alphaShape.cpp:155
pgrouting::graph::Pgr_base_graph::source
V source(E e_idx) const
Definition: pgr_base_graph.hpp:560
pgrouting::alphashape::Pgr_alphaShape::m_adjacent_triangles
std::map< Triangle, std::set< Triangle > > m_adjacent_triangles
t -> {adj_t}
Definition: pgr_alphaShape.h:95
pgrouting::alphashape::Pgr_alphaShape::EdgesFilter::edges
std::set< E > edges
Definition: pgr_alphaShape.h:86
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::get_polygon
Bpoly get_polygon(V source, V target, const std::vector< V > &predecessors, const B_G &graph, std::set< E > &edges_used)
Definition: pgr_alphaShape.cpp:127
pgrouting::alphashape::BG
boost::adjacency_list< boost::setS, boost::vecS, boost::undirectedS, XY_vertex, Basic_edge > BG
Definition: pgr_alphaShape.h:55
pgrouting::alphashape::Pgr_alphaShape::EdgesFilter
Definition: pgr_alphaShape.h:85
pgrouting::Basic_edge::cost
double cost
Definition: basic_edge.h:44
pgr_alphaShape.h
pgrouting::alphashape::Pgr_alphaShape::radius
double radius(const Triangle t) const
Definition: pgr_alphaShape.cpp:227
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::get_predecessors
std::vector< V > get_predecessors(V source, V target, const B_G &subg)
Definition: pgr_alphaShape.cpp:95
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::det
double det(double r00, double r01, double r10, double r11)
Definition: pgr_alphaShape.cpp:68
pgrouting::found_goals
exception for visitor termination
Definition: found_goals.hpp:33
pgrouting::alphashape::V
boost::graph_traits< BG >::vertex_descriptor V
Definition: pgr_alphaShape.h:58
pgrouting::alphashape::Pgr_alphaShape::operator()
std::vector< Bpoly > operator()(double alpha) const
Definition: pgr_alphaShape.cpp:368
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::CompareRadius::operator()
bool operator()(const MyPairType &lhs, const MyPairType &rhs) const
Definition: pgr_alphaShape.cpp:156
pgrouting::alphashape::Pgr_alphaShape
Definition: pgr_alphaShape.h:63
pgrouting
Book keeping class for swapping orders between vehicles.
Definition: pgr_alphaShape.cpp:56
pgrouting::alphashape::Pgr_alphaShape::recursive_build
void recursive_build(const Triangle face, std::set< Triangle > &used, std::set< E > &border_edges, double alpha) const
Definition: pgr_alphaShape.cpp:319
pgrouting::alphashape::anonymous_namespace{pgr_alphaShape.cpp}::circumcenter
Bpoint circumcenter(const Bpoint a, const Bpoint b, const Bpoint c)
Definition: pgr_alphaShape.cpp:73
pgrouting::graph::Pgr_base_graph::adjacent
V adjacent(V v_idx, E e_idx) const
Definition: pgr_base_graph.hpp:562