pgRouting
pgRouting extends the PostGIS / PostgreSQL geospatial database to provide geospatial routing functionality.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 #if defined(__MINGW32__) || defined(_MSC_VER)
38 #include <winsock2.h>
39 #include <windows.h>
40 #endif
41 #ifdef __MINGW64__
42 
43 
44 namespace boost {
45  void tss_cleanup_implemented() { }
46 }
47 #endif
48 #include "./alpha_driver.h"
49 
50 #include <CGAL/Simple_cartesian.h>
51 #include <CGAL/Filtered_kernel.h>
52 #include <CGAL/algorithm.h>
53 
54 #include <CGAL/Polygon_2.h>
55 #include <CGAL/Delaunay_triangulation_2.h>
56 #include <CGAL/Triangulation_2.h>
57 #include <CGAL/Triangulation_hierarchy_vertex_base_2.h>
58 #include <CGAL/Triangulation_hierarchy_2.h>
59 #include <CGAL/Triangulation_face_base_2.h>
60 #include <CGAL/Triangulation_euclidean_traits_2.h>
61 #include <CGAL/Alpha_shape_2.h>
62 #include <CGAL/Alpha_shape_face_base_2.h>
63 #include <CGAL/Alpha_shape_vertex_base_2.h>
64 
65 #include <vector>
66 #include <list>
67 #include <cmath>
68 #include <utility>
69 #include <algorithm>
70 #include <set>
71 
72 
73 typedef double coord_type;
74 
75 typedef CGAL::Simple_cartesian<coord_type> SC;
76 typedef CGAL::Filtered_kernel<SC> K;
77 typedef K::Point_2 Point;
78 typedef K::Segment_2 Segment;
79 typedef K::Vector_2 Vector;
80 typedef CGAL::Polygon_2<K> Polygon_2;
81 
82 typedef CGAL::Alpha_shape_vertex_base_2<K> Avb;
83 typedef CGAL::Triangulation_hierarchy_vertex_base_2<Avb> Av;
84 
85 typedef CGAL::Triangulation_face_base_2<K> Tf;
86 typedef CGAL::Alpha_shape_face_base_2<K, Tf> Af;
87 
88 typedef CGAL::Triangulation_default_data_structure_2<K, Av, Af> Tds;
89 typedef CGAL::Delaunay_triangulation_2<K, Tds> Dt;
90 typedef CGAL::Triangulation_hierarchy_2<Dt> Ht;
91 typedef CGAL::Alpha_shape_2<Ht> Alpha_shape_2;
92 
95 
98 
99 // ---------------------------------------------------------------------
100 
101 double get_angle(Point p, Point q, Point r) {
102  double m_pi(3.14159265358979323846);
103  Vector v1(q, p);
104  Vector v2(q, r);
105  double cross = v1.x() * v2.y() - v1.y() * v2.x();
106  double dot = v1.x() * v2.x() + v1.y() * v2.y();
107  double angle = atan2(cross, dot);
108  if (angle < 0.0) {
109  angle += 2 * m_pi;
110  }
111  return angle;
112 }
113 
114 size_t prev_size = 0;
115 void find_next_edge(Segment s, std::vector<Segment>& segments,
116  std::set<int>& unusedIndexes, std::vector<Polygon_2>& rings) {
117  if (unusedIndexes.empty()
118  || prev_size == unusedIndexes.size()) {
119  return;
120  }
121 
122  prev_size = unusedIndexes.size();
123 
124  Point start = s.source();
125  Point end = s.target();
126  rings.back().push_back(end);
127 
128  std::vector<int> nextIndexes;
129  for (unsigned int i = 0; i < segments.size(); i++) {
130  if (unusedIndexes.find(i) != unusedIndexes.end()) {
131  Point source = segments.at(i).source();
132  if (source == end) {
133  nextIndexes.push_back(i);
134  }
135  }
136  }
137  if (nextIndexes.size() == 1) {
138  int i = nextIndexes.at(0);
139  unusedIndexes.erase(i);
140  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
141  } else if (nextIndexes.size() > 1) {
142  std::vector< std::pair<double, int> > nextAngles;
143  for (unsigned int i = 0; i < nextIndexes.size(); i++) {
144  int j = nextIndexes.at(i);
145  Point target = segments.at(j).target();
146  double angle = get_angle(start, end, target);
147  nextAngles.push_back(std::pair<double, int>(angle, j));
148  }
149  std::sort(nextAngles.begin(), nextAngles.end());
150  int i = nextAngles.begin()->second;
151  unusedIndexes.erase(i);
152  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
153  }
154 
155  if (!unusedIndexes.empty()) {
156  for (unsigned int i = 0; i < segments.size(); i++) {
157  if (unusedIndexes.find(i) != unusedIndexes.end()) {
158  Polygon_2 ring;
159  ring.push_back(segments.at(i).source());
160  rings.push_back(ring);
161  unusedIndexes.erase(i);
162  find_next_edge(segments.at(i), segments, unusedIndexes, rings);
163  }
164  }
165  }
166 }
167 
168 template <class OutputIterator>
169 void
171  OutputIterator out) {
172  for (Alpha_shape_edges_iterator it = A.alpha_shape_edges_begin();
173  it != A.alpha_shape_edges_end();
174  ++it) {
175  *out++ = A.segment(*it);
176  }
177 }
178 
179 
180 int alpha_shape(vertex_t *vertices, size_t count, double alpha,
181  vertex_t **res, size_t *res_count, char **err_msg) {
182  try {
183  std::list<Point> points;
184  {
185  std::vector<Point> pv;
186 
187  for (std::size_t j = 0; j < count; ++j) {
188  Point p(vertices[j].x, vertices[j].y);
189  pv.push_back(p);
190  }
191 
192  std::sort(pv.begin(), pv.end(),
193  [](const Point &e1, const Point &e2)->bool {
194  return e2.y() < e1.y();
195  });
196  std::stable_sort(pv.begin(), pv.end(),
197  [](const Point &e1, const Point &e2)->bool {
198  return e2.x() < e1.x();
199  });
200  pv.erase(std::unique(pv.begin(), pv.end()), pv.end());
201  if (pv.size() != count && pv.size() < 3) {
202  *err_msg = strdup("After eliminating duplicated points, less than 3 points remain!!. Alpha shape calculation needs at least 3 vertices.");
203  return -1;
204  }
205  points.insert(points.begin(), pv.begin(), pv.end());
206  }
207 
208  Alpha_shape_2 A(points.begin(), points.end(),
209  coord_type(10000),
210  Alpha_shape_2::REGULARIZED);
211 
212  std::vector<Segment> segments;
213  // std::vector<Segment> result;
214 
215  // Alpha_shape_2::Alpha_shape_vertices_iterator vit;
216  // Alpha_shape_2::Vertex_handle vertex;
217  // Alpha_shape_2::Alpha_shape_edges_iterator eit;
218  // Alpha_shape_2::Edge edge;
219  // Alpha_shape_2::Face_iterator fit;
220  // Alpha_shape_2::Face_handle face;
221 
222  if (alpha <= 0.0) {
223  alpha = *A.find_optimal_alpha(1);
224  }
225  A.set_alpha(alpha);
226 
227  alpha_edges(A, std::back_inserter(segments));
228 
229  // Segment s = segments.at(0);
230  // find_next_edge(s, segments, result);
231  if (segments.empty()) {
232  *res = NULL;
233  *res_count = 0;
234  } else {
235  std::set<int> unusedIndexes;
236  for (unsigned int i = 0; i < segments.size(); i++) {
237  unusedIndexes.insert(i);
238  }
239 
240  std::vector<Polygon_2> rings;
241  Polygon_2 ring;
242  ring.push_back(segments.at(0).source());
243  rings.push_back(ring);
244  unusedIndexes.erase(0);
245  find_next_edge(segments.at(0), segments, unusedIndexes, rings);
246 
247  size_t result_count = 0;
248  for (unsigned int i = 0; i < rings.size(); i++) {
249  Polygon_2 ring = rings.at(i);
250  result_count += ring.size();
251  }
252  result_count += rings.size() - 1;
253  *res = (vertex_t *) malloc(sizeof(vertex_t) * result_count);
254  *res_count = result_count;
255 
256  int idx = 0;
257  for (unsigned int i = 0; i < rings.size(); i++) {
258  if (i > 0) {
259  (*res)[idx].x = DBL_MAX;
260  (*res)[idx].y = DBL_MAX;
261  idx++;
262  }
263  Polygon_2 ring = rings.at(i);
264  for (unsigned int j = 0; j < ring.size(); j++) {
265  Point point = ring.vertex(j);
266  (*res)[idx].x = point.x();
267  (*res)[idx].y = point.y();
268  idx++;
269  }
270  }
271  }
272  *err_msg = NULL;
273 
274  return EXIT_SUCCESS;
275  } catch ( ... ) {
276  *err_msg = strdup("Caught unknown exception!");
277  }
278  return -1;
279 }
float8 x
Definition: tsp.h:40
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
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
CGAL::Alpha_shape_vertex_base_2< K > Avb
double get_angle(Point p, Point q, Point r)
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
Definition: tsp.h:38
CGAL::Delaunay_triangulation_2< K, Tds > Dt
Alpha_shape_2::Vertex_circulator Vertex_circulator
char * err_msg
Definition: BDATester.cpp:50
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
int alpha_shape(vertex_t *vertices, size_t count, double alpha, vertex_t **res, size_t *res_count, char **err_msg)
double x
Definition: alpha_driver.h:34