43 void tss_cleanup_implemented() { }
46 #include <CGAL/Simple_cartesian.h>
47 #include <CGAL/Filtered_kernel.h>
48 #include <CGAL/algorithm.h>
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>
68 typedef double coord_type;
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;
77 typedef CGAL::Alpha_shape_vertex_base_2<K> Avb;
78 typedef CGAL::Triangulation_hierarchy_vertex_base_2<Avb> Av;
80 typedef CGAL::Triangulation_face_base_2<K> Tf;
81 typedef CGAL::Alpha_shape_face_base_2<K,Tf> Af;
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;
88 typedef Alpha_shape_2::Face_circulator Face_circulator;
89 typedef Alpha_shape_2::Vertex_circulator Vertex_circulator;
91 typedef Alpha_shape_2::Alpha_iterator Alpha_iterator;
92 typedef Alpha_shape_2::Alpha_shape_edges_iterator Alpha_shape_edges_iterator;
98 double m_pi(3.14159265358979323846);
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);
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)
114 if(unusedIndexes.empty()
115 || prev_size == unusedIndexes.size())
120 prev_size = unusedIndexes.size();
122 Point start = s.source();
123 Point end = s.target();
124 rings.back().push_back(end);
126 std::vector<int> nextIndexes;
127 for(
unsigned int i = 0;i < segments.size(); i++)
129 if (unusedIndexes.find(i) != unusedIndexes.end())
131 Point source = segments.at(i).source();
134 nextIndexes.push_back(i);
138 if (nextIndexes.size() == 1)
140 int i = nextIndexes.at(0);
141 unusedIndexes.erase(i);
142 find_next_edge(segments.at(i), segments, unusedIndexes, rings);
144 else if (nextIndexes.size() > 1)
146 std::vector< std::pair<double, int> > nextAngles;
147 for (
unsigned int i = 0; i < nextIndexes.size(); i++)
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));
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);
160 if (!unusedIndexes.empty())
162 for (
unsigned int i = 0; i < segments.size(); i++)
164 if (unusedIndexes.find(i) != unusedIndexes.end())
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);
176 template <
class OutputIterator>
178 alpha_edges(
const Alpha_shape_2& A,
182 for(Alpha_shape_edges_iterator it = A.alpha_shape_edges_begin();
183 it != A.alpha_shape_edges_end();
185 *out++ = A.segment(*it);
190 int alpha_shape(
vertex_t *vertices,
size_t count,
double alpha,
191 vertex_t **res,
size_t *res_count,
char **err_msg)
193 std::list<Point> points;
197 for (std::size_t j = 0; j < count; ++j)
199 Point p(vertices[j].x, vertices[j].y);
204 Alpha_shape_2 A(points.begin(), points.end(),
206 Alpha_shape_2::REGULARIZED);
208 std::vector<Segment> segments;
220 alpha = *A.find_optimal_alpha(1);
224 alpha_edges( A, std::back_inserter(segments));
228 if (segments.empty())
235 std::set<int> unusedIndexes;
236 for (
unsigned int i = 0; i < segments.size(); i++)
238 unusedIndexes.insert(i);
241 std::vector<Polygon_2> rings;
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);
248 size_t result_count = 0;
249 for (
unsigned int i = 0; i < rings.size(); i++)
251 Polygon_2 ring = rings.at(i);
252 result_count += ring.size();
254 result_count += rings.size() - 1;
256 *res_count = result_count;
259 for (
unsigned int i = 0; i < rings.size(); i++)
263 (*res)[idx].x = DBL_MAX;
264 (*res)[idx].y = DBL_MAX;
267 Polygon_2 ring = rings.at(i);
268 for(
unsigned int j = 0; j < ring.size(); j++)
271 (*res)[idx].x = point.x();
272 (*res)[idx].y = point.y();