1
#include "Geometry.hpp"
3
#include "PolylineCollection.hpp"
16
using namespace boost::polygon; // provides also high() and low()
18
namespace Slic3r { namespace Geometry {
21
sort_points (Point a, Point b)
23
return (a.x < b.x) || (a.x == b.x && a.y < b.y);
26
/* This implementation is based on Andrew's monotone chain 2D convex hull algorithm */
28
convex_hull(Points &points, Polygon* hull)
30
assert(points.size() >= 3);
32
std::sort(points.begin(), points.end(), sort_points);
34
int n = points.size(), k = 0;
35
hull->points.resize(2*n);
38
for (int i = 0; i < n; i++) {
39
while (k >= 2 && points[i].ccw(hull->points[k-2], hull->points[k-1]) <= 0) k--;
40
hull->points[k++] = points[i];
44
for (int i = n-2, t = k+1; i >= 0; i--) {
45
while (k >= t && points[i].ccw(hull->points[k-2], hull->points[k-1]) <= 0) k--;
46
hull->points[k++] = points[i];
49
hull->points.resize(k);
51
assert( hull->points.front().coincides_with(hull->points.back()) );
52
hull->points.pop_back();
55
/* accepts an arrayref of points and returns a list of indices
56
according to a nearest-neighbor walk */
58
chained_path(Points &points, std::vector<Points::size_type> &retval, Point start_near)
61
std::map<Point*,Points::size_type> indices;
62
my_points.reserve(points.size());
63
for (Points::iterator it = points.begin(); it != points.end(); ++it) {
64
my_points.push_back(&*it);
65
indices[&*it] = it - points.begin();
68
retval.reserve(points.size());
69
while (!my_points.empty()) {
70
Points::size_type idx = start_near.nearest_point_index(my_points);
71
start_near = *my_points[idx];
72
retval.push_back(indices[ my_points[idx] ]);
73
my_points.erase(my_points.begin() + idx);
78
chained_path(Points &points, std::vector<Points::size_type> &retval)
80
if (points.empty()) return; // can't call front() on empty vector
81
chained_path(points, retval, points.front());
84
/* retval and items must be different containers */
87
chained_path_items(Points &points, T &items, T &retval)
89
std::vector<Points::size_type> indices;
90
chained_path(points, indices);
91
for (std::vector<Points::size_type>::const_iterator it = indices.begin(); it != indices.end(); ++it)
92
retval.push_back(items[*it]);
94
template void chained_path_items(Points &points, ClipperLib::PolyNodes &items, ClipperLib::PolyNodes &retval);
97
directions_parallel(double angle1, double angle2, double max_diff)
99
double diff = fabs(angle1 - angle2);
101
return diff < max_diff || fabs(diff - PI) < max_diff;
105
MedialAxis::edge_to_line(const VD::edge_type &edge) const
108
line.a.x = edge.vertex0()->x();
109
line.a.y = edge.vertex0()->y();
110
line.b.x = edge.vertex1()->x();
111
line.b.y = edge.vertex1()->y();
116
MedialAxis::build(Polylines* polylines)
119
// build bounding box (we use it for clipping infinite segments)
120
// --> we have no infinite segments
121
this->bb = BoundingBox(this->lines);
124
construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd);
127
// DEBUG: dump all Voronoi edges
129
for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
130
if (edge->is_infinite()) continue;
133
polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
134
polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
135
polylines->push_back(polyline);
141
// collect valid edges (i.e. prune those not belonging to MAT)
142
// note: this keeps twins, so it contains twice the number of the valid edges
144
for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
145
// if we only process segments representing closed loops, none if the
146
// infinite edges (if any) would be part of our MAT anyway
147
if (edge->is_secondary() || edge->is_infinite()) continue;
148
this->edges.insert(&*edge);
151
// count valid segments for each vertex
152
std::map< const VD::vertex_type*,std::set<const VD::edge_type*> > vertex_edges;
153
std::set<const VD::vertex_type*> entry_nodes;
154
for (VD::const_vertex_iterator vertex = this->vd.vertices().begin(); vertex != this->vd.vertices().end(); ++vertex) {
155
// get a reference to the list of valid edges originating from this vertex
156
std::set<const VD::edge_type*>& edges = vertex_edges[&*vertex];
158
// get one random edge originating from this vertex
159
const VD::edge_type* edge = vertex->incident_edge();
161
if (this->edges.count(edge) > 0) // only count valid edges
163
edge = edge->rot_next(); // next edge originating from this vertex
164
} while (edge != vertex->incident_edge());
166
// if there's only one edge starting at this vertex then it's a leaf
167
size_t edge_count = edges.size();
168
if (edge_count == 1) {
169
entry_nodes.insert(&*vertex);
174
while (!entry_nodes.empty()) {
175
// get a random entry node
176
const VD::vertex_type* v = *entry_nodes.begin();
178
// get edge starting from v
179
assert(!vertex_edges[v].empty());
180
const VD::edge_type* edge = *vertex_edges[v].begin();
182
if (!this->is_valid_edge(*edge)) {
183
// if edge is not valid, erase it from edge list
184
(void)this->edges.erase(edge);
185
(void)this->edges.erase(edge->twin());
187
// decrement edge counters for the affected nodes
188
const VD::vertex_type* v1 = edge->vertex1();
189
(void)vertex_edges[v].erase(edge);
190
(void)vertex_edges[v1].erase(edge->twin());
192
// also, check whether the end vertex is a new leaf
193
if (vertex_edges[v1].size() == 1) {
194
entry_nodes.insert(v1);
195
} else if (vertex_edges[v1].empty()) {
196
entry_nodes.erase(v1);
200
// remove node from the set to prevent it from being visited again
201
entry_nodes.erase(v);
204
// iterate through the valid edges to build polylines
205
while (!this->edges.empty()) {
206
const VD::edge_type& edge = **this->edges.begin();
210
polyline.points.push_back(Point( edge.vertex0()->x(), edge.vertex0()->y() ));
211
polyline.points.push_back(Point( edge.vertex1()->x(), edge.vertex1()->y() ));
213
// remove this edge and its twin from the available edges
214
(void)this->edges.erase(&edge);
215
(void)this->edges.erase(edge.twin());
218
this->process_edge_neighbors(edge, &polyline.points);
220
// get previous points
222
this->process_edge_neighbors(*edge.twin(), &pp);
223
polyline.points.insert(polyline.points.begin(), pp.rbegin(), pp.rend());
225
// append polyline to result if it's not too small
226
if (polyline.length() > this->max_width)
227
polylines->push_back(polyline);
232
MedialAxis::process_edge_neighbors(const VD::edge_type& edge, Points* points)
234
// Since rot_next() works on the edge starting point but we want
235
// to find neighbors on the ending point, we just swap edge with
237
const VD::edge_type& twin = *edge.twin();
239
// count neighbors for this edge
240
std::vector<const VD::edge_type*> neighbors;
241
for (const VD::edge_type* neighbor = twin.rot_next(); neighbor != &twin; neighbor = neighbor->rot_next()) {
242
if (this->edges.count(neighbor) > 0) neighbors.push_back(neighbor);
245
// if we have a single neighbor then we can continue recursively
246
if (neighbors.size() == 1) {
247
const VD::edge_type& neighbor = *neighbors.front();
248
points->push_back(Point( neighbor.vertex1()->x(), neighbor.vertex1()->y() ));
249
(void)this->edges.erase(&neighbor);
250
(void)this->edges.erase(neighbor.twin());
251
this->process_edge_neighbors(neighbor, points);
256
MedialAxis::is_valid_edge(const VD::edge_type& edge) const
258
/* If the cells sharing this edge have a common vertex, we're not interested
259
in this edge. Why? Because it means that the edge lies on the bisector of
260
two contiguous input lines and it was included in the Voronoi graph because
261
it's the locus of centers of circles tangent to both vertices. Due to the
262
"thin" nature of our input, these edges will be very short and not part of
263
our wanted output. */
265
const VD::cell_type &cell1 = *edge.cell();
266
const VD::cell_type &cell2 = *edge.twin()->cell();
267
if (cell1.contains_segment() && cell2.contains_segment()) {
268
Line segment1 = this->retrieve_segment(cell1);
269
Line segment2 = this->retrieve_segment(cell2);
270
if (segment1.a == segment2.b || segment1.b == segment2.a) return false;
272
// calculate relative angle between the two boundary segments
273
Vector vec1 = segment1.vector();
274
Vector vec2 = segment2.vector();
275
double angle = atan2(vec1.x*vec2.y - vec1.y*vec2.x, vec1.x*vec2.x + vec1.y*vec2.y);
277
// fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
278
// we're interested only in segments close to the second case (facing segments)
279
// so we allow some tolerance (say, 30°)
280
if (fabs(angle) < PI - PI/3) return false;
283
// each vertex is equidistant to both cell segments
284
// but such distance might differ between the two vertices;
285
// in this case it means the shape is getting narrow (like a corner)
286
// and we might need to skip the edge since it's not really part of
288
Point v0( edge.vertex0()->x(), edge.vertex0()->y() );
289
Point v1( edge.vertex1()->x(), edge.vertex1()->y() );
290
double dist0 = v0.distance_to(segment1);
291
double dist1 = v1.distance_to(segment1);
294
double diff = fabs(dist1 - dist0);
295
double dist_between_segments1 = segment1.a.distance_to(segment2);
296
double dist_between_segments2 = segment1.b.distance_to(segment2);
297
printf("w = %f/%f, dist0 = %f, dist1 = %f, diff = %f, seg1len = %f, seg2len = %f, edgelen = %f, s2s = %f / %f\n",
298
unscale(this->max_width), unscale(this->min_width),
299
unscale(dist0), unscale(dist1), unscale(diff),
300
unscale(segment1.length()), unscale(segment2.length()),
301
unscale(this->edge_to_line(edge).length()),
302
unscale(dist_between_segments1), unscale(dist_between_segments2)
306
// if this segment is the centerline for a very thin area, we might want to skip it
307
// in case the area is too thin
308
if (dist0 < this->min_width/2 || dist1 < this->min_width/2) {
309
//printf(" => too thin, skipping\n");
314
// if distance between this edge and the thin area boundary is greater
315
// than half the max width, then it's not a true medial axis segment
316
if (dist1 > this->width*2) {
317
printf(" => too fat, skipping\n");
329
MedialAxis::retrieve_segment(const VD::cell_type& cell) const
331
VD::cell_type::source_index_type index = cell.source_index() - this->points.size();
332
return this->lines[index];
added
removed
xs/src/Geometry.cpp
