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