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