PGROUTING  2.6-dev
alpha_driver.cpp
Go to the documentation of this file.
1 /*PGR-GNU*****************************************************************
2 File: alpha_drivedist.cpp
3 
4 Copyright (c) 2006 Anton A. Patrushev, Orkney, Inc.
5 Mail: project@pgrouting.org
6 
7 ------
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 2 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, write to the Free Software
21 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
22 
23 ********************************************************************PGR-GNU*/
24 /*
25  * As a special exception, you have permission to link this program
26  * with the CGAL library and distribute executables, as long as you
27  * follow the requirements of the GNU GPL in regard to all of the
28  * software in the executable aside from CGAL.
29  *
30  */
31 
32 
33 /***********************************************************************
34 Takes a list of points and returns a list of segments
35 corresponding to the Alpha shape.
36 ************************************************************************/
37 
39 
40 #include <vector>
41 #include <list>
42 #include <cmath>
43 #include <utility>
44 #include <algorithm>
45 #include <set>
46 
47 #include "cpp_common/pgr_assert.h"
48 
49 #include "CGAL/Simple_cartesian.h"
50 #include "CGAL/Filtered_kernel.h"
51 #include "CGAL/algorithm.h"
52 
53 #include "CGAL/Polygon_2.h"
54 #include "CGAL/Delaunay_triangulation_2.h"
55 #include "CGAL/Triangulation_2.h"
56 #include "CGAL/Triangulation_hierarchy_vertex_base_2.h"
57 #include "CGAL/Triangulation_hierarchy_2.h"
58 #include "CGAL/Triangulation_face_base_2.h"
59 #include "CGAL/Triangulation_euclidean_traits_2.h"
60 #include "CGAL/Alpha_shape_2.h"
61 #include "CGAL/Alpha_shape_face_base_2.h"
62 #include "CGAL/Alpha_shape_vertex_base_2.h"
63 
64 
65 #include "cpp_common/pgr_alloc.hpp"
66 
67 typedef double coord_type;
68 
69 typedef CGAL::Simple_cartesian<coord_type> SC;
70 typedef CGAL::Filtered_kernel<SC> K;
71 typedef K::Point_2 Point;
72 typedef K::Segment_2 Segment;
73 typedef K::Vector_2 Vector;
74 typedef CGAL::Polygon_2<K> Polygon_2;
75 
76 typedef CGAL::Alpha_shape_vertex_base_2<K> Avb;
77 typedef CGAL::Triangulation_hierarchy_vertex_base_2<Avb> Av;
78 
79 typedef CGAL::Triangulation_face_base_2<K> Tf;
80 typedef CGAL::Alpha_shape_face_base_2<K, Tf> Af;
81 
82 typedef CGAL::Triangulation_default_data_structure_2<K, Av, Af> Tds;
83 typedef CGAL::Delaunay_triangulation_2<K, Tds> Dt;
84 typedef CGAL::Triangulation_hierarchy_2<Dt> Ht;
85 typedef CGAL::Alpha_shape_2<Ht> Alpha_shape_2;
86 
89 
92 
93 // ---------------------------------------------------------------------
94 
95 double get_angle(Point p, Point q, Point r) {
96  double m_pi(3.14159265358979323846);
97  Vector v1(q, p);
98  Vector v2(q, r);
99  double cross = v1.x() * v2.y() - v1.y() * v2.x();
100  double dot = v1.x() * v2.x() + v1.y() * v2.y();
101  double angle = atan2(cross, dot);
102  if (angle < 0.0) {
103  angle += 2 * m_pi;
104  }
105  return angle;
106 }
107 
108 size_t prev_size = 0;
109 void find_next_edge(Segment s, std::vector<Segment>& segments,
110  std::set<int>& unusedIndexes, std::vector<Polygon_2>& rings) {
111  if (unusedIndexes.empty()
112  || prev_size == unusedIndexes.size()) {
113  return;
114  }
115 
116  prev_size = unusedIndexes.size();
117 
118  Point start = s.source();
119  Point end = s.target();
120  rings.back().push_back(end);
121 
122  std::vector<int> nextIndexes;
123  for (unsigned int i = 0; i < segments.size(); i++) {
124  if (unusedIndexes.find(i) != unusedIndexes.end()) {
125  Point source = segments.at(i).source();
126  if (source == end) {
127  nextIndexes.push_back(i);
128  }
129  }
130  }
131  if (nextIndexes.size() == 1) {
132  int i = nextIndexes.at(0);
133  unusedIndexes.erase(i);
134  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
135  } else if (nextIndexes.size() > 1) {
136  std::vector< std::pair<double, int> > nextAngles;
137  for (unsigned int i = 0; i < nextIndexes.size(); i++) {
138  int j = nextIndexes.at(i);
139  Point target = segments.at(j).target();
140  double angle = get_angle(start, end, target);
141  nextAngles.push_back(std::pair<double, int>(angle, j));
142  }
143  std::sort(nextAngles.begin(), nextAngles.end());
144  int i = nextAngles.begin()->second;
145  unusedIndexes.erase(i);
146  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
147  }
148 
149  if (!unusedIndexes.empty()) {
150  for (unsigned int i = 0; i < segments.size(); i++) {
151  if (unusedIndexes.find(i) != unusedIndexes.end()) {
152  Polygon_2 ring;
153  ring.push_back(segments.at(i).source());
154  rings.push_back(ring);
155  unusedIndexes.erase(i);
156  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
157  }
158  }
159  }
160 }
161 
162 template <class OutputIterator>
163 void
165  OutputIterator out) {
166  for (Alpha_shape_edges_iterator it = A.alpha_shape_edges_begin();
167  it != A.alpha_shape_edges_end();
168  ++it) {
169  *out++ = A.segment(*it);
170  }
171 }
172 
173 
174 int alpha_shape(vertex_t *vertices, size_t count, double alpha,
175  vertex_t **return_tuples, size_t *res_count, char **err_msg) {
176  std::ostringstream err;
177 
178  try {
179  std::list<Point> points;
180  {
181  std::vector<Point> pv;
182 
183  for (std::size_t j = 0; j < count; ++j) {
184  Point p(vertices[j].x, vertices[j].y);
185  pv.push_back(p);
186  }
187 
188  std::sort(pv.begin(), pv.end(),
189  [](const Point &e1, const Point &e2)->bool {
190  return e2.y() < e1.y();
191  });
192  std::stable_sort(pv.begin(), pv.end(),
193  [](const Point &e1, const Point &e2)->bool {
194  return e2.x() < e1.x();
195  });
196  pv.erase(std::unique(pv.begin(), pv.end()), pv.end());
197  if (pv.size() != count && pv.size() < 3) {
198  err << "After eliminating duplicated points,"
199  " less than 3 points remain!!."
200  " Alpha shape calculation needs at least 3 vertices.";
201  *err_msg = pgr_msg(err.str().c_str());
202  return -1;
203  }
204  points.insert(points.begin(), pv.begin(), pv.end());
205  }
206 
207  Alpha_shape_2 A(points.begin(), points.end(),
208  coord_type(10000),
209  Alpha_shape_2::REGULARIZED);
210 
211  std::vector<Segment> segments;
212  // std::vector<Segment> result;
213 
214  // Alpha_shape_2::Alpha_shape_vertices_iterator vit;
215  // Alpha_shape_2::Vertex_handle vertex;
216  // Alpha_shape_2::Alpha_shape_edges_iterator eit;
217  // Alpha_shape_2::Edge edge;
218  // Alpha_shape_2::Face_iterator fit;
219  // Alpha_shape_2::Face_handle face;
220 
221  if (alpha <= 0.0) {
222  alpha = *A.find_optimal_alpha(1);
223  }
224  A.set_alpha(alpha);
225 
226  alpha_edges(A, std::back_inserter(segments));
227 
228  // Segment s = segments.at(0);
229  // find_next_edge(s, segments, result);
230  if (segments.empty()) {
231  *return_tuples = NULL;
232  *res_count = 0;
233  } else {
234  std::set<int> unusedIndexes;
235  for (unsigned int i = 0; i < segments.size(); i++) {
236  unusedIndexes.insert(i);
237  }
238 
239  std::vector<Polygon_2> rings;
240  Polygon_2 ring;
241  ring.push_back(segments.at(0).source());
242  rings.push_back(ring);
243  unusedIndexes.erase(0);
244  find_next_edge(segments.at(0), segments, unusedIndexes, rings);
245 
246  size_t result_count = 0;
247  for (unsigned int i = 0; i < rings.size(); i++) {
248  Polygon_2 ring = rings.at(i);
249  result_count += ring.size();
250  }
251  result_count += rings.size() - 1;
252  *return_tuples = pgr_alloc(result_count, (*return_tuples));
253  *res_count = result_count;
254 
255  int idx = 0;
256  for (unsigned int i = 0; i < rings.size(); i++) {
257  if (i > 0) {
258  (*return_tuples)[idx].x = DBL_MAX;
259  (*return_tuples)[idx].y = DBL_MAX;
260  idx++;
261  }
262  Polygon_2 ring = rings.at(i);
263  for (unsigned int j = 0; j < ring.size(); j++) {
264  Point point = ring.vertex(j);
265  (*return_tuples)[idx].x = point.x();
266  (*return_tuples)[idx].y = point.y();
267  idx++;
268  }
269  }
270  }
271  *err_msg = NULL;
272 
273  return EXIT_SUCCESS;
274  } catch (AssertFailedException &except) {
275  (*return_tuples) = pgr_free(*return_tuples);
276  (*res_count) = 0;
277  err << except.what();
278  *err_msg = pgr_msg(err.str().c_str());
279  } catch (std::exception &except) {
280  (*return_tuples) = pgr_free(*return_tuples);
281  (*res_count) = 0;
282  err << except.what();
283  *err_msg = pgr_msg(err.str().c_str());
284  } catch(...) {
285  (*return_tuples) = pgr_free(*return_tuples);
286  (*res_count) = 0;
287  err << "Caught unknown exception!";
288  *err_msg = pgr_msg(err.str().c_str());
289  }
290  return -1;
291 }
CGAL::Triangulation_face_base_2< K > Tf
size_t prev_size
K::Point_2 Point
CGAL::Triangulation_default_data_structure_2< K, Av, Af > Tds
Extends std::exception and is the exception that we throw if an assert fails.
Definition: pgr_assert.h:126
void alpha_edges(const Alpha_shape_2 &A, OutputIterator out)
K::Segment_2 Segment
CGAL::Triangulation_hierarchy_vertex_base_2< Avb > Av
CGAL::Simple_cartesian< coord_type > SC
int alpha_shape(vertex_t *vertices, size_t count, double alpha, vertex_t **return_tuples, size_t *res_count, char **err_msg)
CGAL::Alpha_shape_vertex_base_2< K > Avb
double get_angle(Point p, Point q, Point r)
T * pgr_free(T *ptr)
Definition: pgr_alloc.hpp:77
CGAL::Filtered_kernel< SC > K
CGAL::Alpha_shape_2< Ht > Alpha_shape_2
CGAL::Triangulation_hierarchy_2< Dt > Ht
double coord_type
Alpha_shape_2::Face_circulator Face_circulator
CGAL::Polygon_2< K > Polygon_2
char * pgr_msg(const std::string &msg)
Definition: pgr_alloc.cpp:30
T * pgr_alloc(std::size_t size, T *ptr)
allocates memory
Definition: pgr_alloc.hpp:66
Assertions Handling.
CGAL::Delaunay_triangulation_2< K, Tds > Dt
virtual const char * what() const
Definition: pgr_assert.cpp:53
Alpha_shape_2::Vertex_circulator Vertex_circulator
Alpha_shape_2::Alpha_iterator Alpha_iterator
K::Vector_2 Vector
void find_next_edge(Segment s, std::vector< Segment > &segments, std::set< int > &unusedIndexes, std::vector< Polygon_2 > &rings)
Alpha_shape_2::Alpha_shape_edges_iterator Alpha_shape_edges_iterator
CGAL::Alpha_shape_face_base_2< K, Tf > Af