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