diff --git a/lib/Slic3r/Geometry.pm b/lib/Slic3r/Geometry.pm index e48e56bb..d647f2a8 100644 --- a/lib/Slic3r/Geometry.pm +++ b/lib/Slic3r/Geometry.pm @@ -5,7 +5,7 @@ use warnings; require Exporter; our @ISA = qw(Exporter); our @EXPORT_OK = qw( - PI X Y Z A B X1 Y1 X2 Y2 Z1 Z2 MIN MAX epsilon slope line_atan lines_parallel + PI X Y Z A B X1 Y1 X2 Y2 Z1 Z2 MIN MAX epsilon slope line_point_belongs_to_segment points_coincide distance_between_points normalize tan move_points_3D point_in_polygon point_in_segment segment_in_segment @@ -15,7 +15,7 @@ our @EXPORT_OK = qw( rotate_points move_points dot perp polygon_points_visibility line_intersection bounding_box bounding_box_intersect - angle3points three_points_aligned line_direction + angle3points chained_path chained_path_from collinear scale unscale rad2deg_dir bounding_box_center line_intersects_any douglas_peucker polyline_remove_short_segments normal triangle_normal polygon_is_convex @@ -57,30 +57,6 @@ sub slope { return ($line->[B][Y] - $line->[A][Y]) / ($line->[B][X] - $line->[A][X]); } -sub line_atan { - my ($line) = @_; - return atan2($line->[B][Y] - $line->[A][Y], $line->[B][X] - $line->[A][X]); -} - -sub line_direction { - my ($line) = @_; - my $atan2 = line_atan($line); - return ($atan2 == PI) ? 0 - : ($atan2 < 0) ? ($atan2 + PI) - : $atan2; -} - -sub lines_parallel { - my ($line1, $line2) = @_; - - return abs(line_direction($line1) - line_direction($line2)) < $parallel_degrees_limit; -} - -sub three_points_aligned { - my ($p1, $p2, $p3) = @_; - return lines_parallel([$p1, $p2], [$p2, $p3]); -} - # this subroutine checks whether a given point may belong to a given # segment given the hypothesis that it belongs to the line containing # the segment diff --git a/lib/Slic3r/Layer/Region.pm b/lib/Slic3r/Layer/Region.pm index 37cfbdf7..f4c425a1 100644 --- a/lib/Slic3r/Layer/Region.pm +++ b/lib/Slic3r/Layer/Region.pm @@ -224,7 +224,7 @@ sub make_perimeters { # process thin walls by collapsing slices to single passes my $min_thin_wall_width = $pwidth/3; my $min_thin_wall_length = 2*$pwidth; - @thin_walls = @{offset2_ex([ map @$_, @thin_walls ], -0.5*$min_thin_wall_width, +0.5*$min_thin_wall_width)}; + #@thin_walls = @{offset2_ex([ map @$_, @thin_walls ], -0.5*$min_thin_wall_width, +0.5*$min_thin_wall_width)}; if (@thin_walls) { my @p = map @{$_->medial_axis($pspacing)}, @thin_walls; diff --git a/lib/Slic3r/Line.pm b/lib/Slic3r/Line.pm index 034770c7..a50e7064 100644 --- a/lib/Slic3r/Line.pm +++ b/lib/Slic3r/Line.pm @@ -7,16 +7,6 @@ use parent 'Slic3r::Polyline'; use Slic3r::Geometry qw(A B X Y); -sub atan { - my $self = shift; - return Slic3r::Geometry::line_atan($self); -} - -sub direction { - my $self = shift; - return Slic3r::Geometry::line_direction($self); -} - sub intersection { my $self = shift; my ($line, $require_crossing) = @_; diff --git a/t/angles.t b/t/angles.t index dc1a671b..3b5d64e2 100644 --- a/t/angles.t +++ b/t/angles.t @@ -10,7 +10,7 @@ BEGIN { } use Slic3r; -use Slic3r::Geometry qw(line_atan line_direction rad2deg_dir angle3points PI); +use Slic3r::Geometry qw(rad2deg_dir angle3points PI); #========================================================== @@ -61,3 +61,13 @@ use Slic3r::Geometry qw(line_atan line_direction rad2deg_dir angle3points PI); } #========================================================== + +sub line_atan { + my ($l) = @_; + return Slic3r::Line->new(@$l)->atan2_; +} + +sub line_direction { + my ($l) = @_; + return Slic3r::Line->new(@$l)->direction; +} \ No newline at end of file diff --git a/xs/src/ExPolygon.cpp b/xs/src/ExPolygon.cpp index 3cb54583..0eb1c8d5 100644 --- a/xs/src/ExPolygon.cpp +++ b/xs/src/ExPolygon.cpp @@ -135,10 +135,10 @@ ExPolygon::simplify(double tolerance, ExPolygons &expolygons) const } void -ExPolygon::medial_axis(Polylines* polylines) const +ExPolygon::medial_axis(double width, Polylines* polylines) const { // init helper object - Slic3r::Geometry::MedialAxis ma; + Slic3r::Geometry::MedialAxis ma(width); // populate list of segments for the Voronoi diagram this->contour.lines(&ma.lines); diff --git a/xs/src/ExPolygon.hpp b/xs/src/ExPolygon.hpp index 0108f566..e7db3e62 100644 --- a/xs/src/ExPolygon.hpp +++ b/xs/src/ExPolygon.hpp @@ -26,7 +26,7 @@ class ExPolygon Polygons simplify_p(double tolerance) const; ExPolygons simplify(double tolerance) const; void simplify(double tolerance, ExPolygons &expolygons) const; - void medial_axis(Polylines* polylines) const; + void medial_axis(double width, Polylines* polylines) const; #ifdef SLIC3RXS void from_SV(SV* poly_sv); diff --git a/xs/src/Geometry.cpp b/xs/src/Geometry.cpp index bc491bee..baf3803b 100644 --- a/xs/src/Geometry.cpp +++ b/xs/src/Geometry.cpp @@ -3,10 +3,10 @@ #include "PolylineCollection.hpp" #include "clipper.hpp" #include +#include #include #include #include -//#include "voronoi_visual_utils.hpp" #ifdef SLIC3R_DEBUG #include "SVG.hpp" @@ -92,6 +92,16 @@ chained_path_items(Points &points, T &items, T &retval) } template void chained_path_items(Points &points, ClipperLib::PolyNodes &items, ClipperLib::PolyNodes &retval); +Line +MedialAxis::edge_to_line(const VD::edge_type &edge) { + Line line; + line.a.x = edge.vertex0()->x(); + line.a.y = edge.vertex0()->y(); + line.b.x = edge.vertex1()->x(); + line.b.y = edge.vertex1()->y(); + return line; +} + void MedialAxis::build(Polylines* polylines) { @@ -104,11 +114,73 @@ MedialAxis::build(Polylines* polylines) construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd); // collect valid edges (i.e. prune those not belonging to MAT) + // note: this keeps twins, so it contains twice the number of the valid edges this->edges.clear(); for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) { if (this->is_valid_edge(*edge)) this->edges.insert(&*edge); } + // count valid segments for each vertex + std::map< const VD::vertex_type*,std::list > vertex_edges; + std::list entry_nodes; + for (VD::const_vertex_iterator vertex = this->vd.vertices().begin(); vertex != this->vd.vertices().end(); ++vertex) { + // get a reference to the list of valid edges originating from this vertex + std::list& edges = vertex_edges[&*vertex]; + + // get one random edge originating from this vertex + const VD::edge_type* edge = vertex->incident_edge(); + do { + if (this->edges.count(edge) > 0) // only count valid edges + edges.push_back(edge); + edge = edge->rot_next(); // next edge originating from this vertex + } while (edge != vertex->incident_edge()); + + // if there's only one edge starting at this vertex then it's a leaf + if (edges.size() == 1) entry_nodes.push_back(&*vertex); + } + + // iterate through the leafs to prune short branches + for (std::list::const_iterator vertex = entry_nodes.begin(); vertex != entry_nodes.end(); ++vertex) { + const VD::vertex_type* v = *vertex; + + // start a polyline from this vertex + Polyline polyline; + polyline.points.push_back(Point(v->x(), v->y())); + + // keep track of visited edges to prevent infinite loops + std::set visited_edges; + + do { + // get edge starting from v + const VD::edge_type* edge = vertex_edges[v].front(); + + // if we picked the edge going backwards (thus the twin of the previous edge) + if (visited_edges.count(edge->twin()) > 0) { + edge = vertex_edges[v].back(); + } + + // avoid getting twice on the same edge + if (visited_edges.count(edge) > 0) break; + visited_edges.insert(edge); + + // get ending vertex for this edge and append it to the polyline + v = edge->vertex1(); + polyline.points.push_back(Point( v->x(), v->y() )); + + // if two edges start at this vertex (one forward one backwards) then + // it's not branching and we can go on + } while (vertex_edges[v].size() == 2); + + // if this branch is too short, invalidate all of its edges so that + // they will be ignored when building actual polylines in the loop below + if (polyline.length() < this->width) { + for (std::set::const_iterator edge = visited_edges.begin(); edge != visited_edges.end(); ++edge) { + (void)this->edges.erase(*edge); + (void)this->edges.erase((*edge)->twin()); + } + } + } + // iterate through the valid edges to build polylines while (!this->edges.empty()) { const VD::edge_type& edge = **this->edges.begin(); @@ -176,107 +248,26 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const but I don't know how to do it. Maybe we could check the relative angle of the two segments (we are only interested in facing segments). */ - const voronoi_diagram::cell_type &cell1 = *edge.cell(); - const voronoi_diagram::cell_type &cell2 = *edge.twin()->cell(); + const VD::cell_type &cell1 = *edge.cell(); + const VD::cell_type &cell2 = *edge.twin()->cell(); if (cell1.contains_segment() && cell2.contains_segment()) { Line segment1 = this->retrieve_segment(cell1); Line segment2 = this->retrieve_segment(cell2); if (segment1.a == segment2.b || segment1.b == segment2.a) return false; + if (fabs(segment1.atan2_() - segment2.atan2_()) < PI/3) { + //printf("segment1 atan2 = %f, segment2 atan2 = %f\n", segment1.atan2_(), segment2.atan2_()); + //printf(" => SAME ATAN2\n"); + return false; + } } return true; } -/* -void -MedialAxis::clip_infinite_edge(const voronoi_diagram::edge_type& edge, Points* clipped_edge) -{ - const voronoi_diagram::cell_type& cell1 = *edge.cell(); - const voronoi_diagram::cell_type& cell2 = *edge.twin()->cell(); - Point origin, direction; - // Infinite edges could not be created by two segment sites. - if (cell1.contains_point() && cell2.contains_point()) { - Point p1 = retrieve_point(cell1); - Point p2 = retrieve_point(cell2); - origin.x = (p1.x + p2.x) * 0.5; - origin.y = (p1.y + p2.y) * 0.5; - direction.x = p1.y - p2.y; - direction.y = p2.x - p1.x; - } else { - origin = cell1.contains_segment() - ? retrieve_point(cell2) - : retrieve_point(cell1); - Line segment = cell1.contains_segment() - ? retrieve_segment(cell1) - : retrieve_segment(cell2); - coord_t dx = high(segment).x - low(segment).x; - coord_t dy = high(segment).y - low(segment).y; - if ((low(segment) == origin) ^ cell1.contains_point()) { - direction.x = dy; - direction.y = -dx; - } else { - direction.x = -dy; - direction.y = dx; - } - } - coord_t side = this->bb.size().x; - coord_t koef = side / (std::max)(fabs(direction.x), fabs(direction.y)); - if (edge.vertex0() == NULL) { - clipped_edge->push_back(Point( - origin.x - direction.x * koef, - origin.y - direction.y * koef - )); - } else { - clipped_edge->push_back( - Point(edge.vertex0()->x(), edge.vertex0()->y())); - } - if (edge.vertex1() == NULL) { - clipped_edge->push_back(Point( - origin.x + direction.x * koef, - origin.y + direction.y * koef - )); - } else { - clipped_edge->push_back( - Point(edge.vertex1()->x(), edge.vertex1()->y())); - } -} - -void -MedialAxis::sample_curved_edge(const voronoi_diagram::edge_type& edge, Points* sampled_edge) -{ - Point point = edge.cell()->contains_point() - ? retrieve_point(*edge.cell()) - : retrieve_point(*edge.twin()->cell()); - - Line segment = edge.cell()->contains_point() - ? retrieve_segment(*edge.twin()->cell()) - : retrieve_segment(*edge.cell()); - - double max_dist = 1E-3 * this->bb.size().x; - voronoi_visual_utils::discretize(point, segment, max_dist, sampled_edge); -} -*/ - -Point -MedialAxis::retrieve_point(const voronoi_diagram::cell_type& cell) -{ - voronoi_diagram::cell_type::source_index_type index = cell.source_index(); - voronoi_diagram::cell_type::source_category_type category = cell.source_category(); - if (category == SOURCE_CATEGORY_SINGLE_POINT) { - return this->points[index]; - } - index -= this->points.size(); - if (category == SOURCE_CATEGORY_SEGMENT_START_POINT) { - return low(this->lines[index]); - } else { - return high(this->lines[index]); - } -} - Line -MedialAxis::retrieve_segment(const voronoi_diagram::cell_type& cell) const +MedialAxis::retrieve_segment(const VD::cell_type& cell) const { - voronoi_diagram::cell_type::source_index_type index = cell.source_index() - this->points.size(); + VD::cell_type::source_index_type index = cell.source_index() - this->points.size(); return this->lines[index]; } diff --git a/xs/src/Geometry.hpp b/xs/src/Geometry.hpp index 6094da9b..b194e6ab 100644 --- a/xs/src/Geometry.hpp +++ b/xs/src/Geometry.hpp @@ -20,19 +20,18 @@ class MedialAxis { public: Points points; Lines lines; + double width; + MedialAxis(double _width) : width(_width) {}; void build(Polylines* polylines); - void process_edge_neighbors(const voronoi_diagram::edge_type& edge, Points* points); - bool is_valid_edge(const voronoi_diagram::edge_type& edge) const; - //void clip_infinite_edge(const voronoi_diagram::edge_type& edge, Points* clipped_edge); - //void sample_curved_edge(const voronoi_diagram::edge_type& edge, Points* sampled_edge); - Point retrieve_point(const voronoi_diagram::cell_type& cell); - Line retrieve_segment(const voronoi_diagram::cell_type& cell) const; private: typedef voronoi_diagram VD; VD vd; - //BoundingBox bb; std::set edges; + Line edge_to_line(const VD::edge_type &edge); + void process_edge_neighbors(const voronoi_diagram::edge_type& edge, Points* points); + bool is_valid_edge(const voronoi_diagram::edge_type& edge) const; + Line retrieve_segment(const voronoi_diagram::cell_type& cell) const; }; } } diff --git a/xs/src/Line.cpp b/xs/src/Line.cpp index 16a3241a..5daceb44 100644 --- a/xs/src/Line.cpp +++ b/xs/src/Line.cpp @@ -1,6 +1,7 @@ #include "Line.hpp" #include "Polyline.hpp" #include +#include #include namespace Slic3r { @@ -85,6 +86,21 @@ Line::distance_to(const Point* point) const return point->distance_to(this); } +double +Line::atan2_() const +{ + return atan2(this->b.y - this->a.y, this->b.x - this->a.x); +} + +double +Line::direction() const +{ + double atan2 = this->atan2_(); + return (atan2 == PI) ? 0 + : (atan2 < 0) ? (atan2 + PI) + : atan2; +} + #ifdef SLIC3RXS void Line::from_SV(SV* line_sv) diff --git a/xs/src/Line.hpp b/xs/src/Line.hpp index e0298216..a3983cdb 100644 --- a/xs/src/Line.hpp +++ b/xs/src/Line.hpp @@ -28,6 +28,8 @@ class Line Point* point_at(double distance) const; bool coincides_with(const Line* line) const; double distance_to(const Point* point) const; + double atan2_() const; + double direction() const; #ifdef SLIC3RXS void from_SV(SV* line_sv); diff --git a/xs/src/voronoi_visual_utils.hpp b/xs/src/voronoi_visual_utils.hpp deleted file mode 100644 index 6ae93e32..00000000 --- a/xs/src/voronoi_visual_utils.hpp +++ /dev/null @@ -1,186 +0,0 @@ -// Boost.Polygon library voronoi_graphic_utils.hpp header file - -// Copyright Andrii Sydorchuk 2010-2012. -// Distributed under the Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -// See http://www.boost.org for updates, documentation, and revision history. - -#ifndef BOOST_POLYGON_VORONOI_VISUAL_UTILS -#define BOOST_POLYGON_VORONOI_VISUAL_UTILS - -#include -#include - -#include -#include -#include -#include - -namespace boost { -namespace polygon { -// Utilities class, that contains set of routines handful for visualization. -template -class voronoi_visual_utils { - public: - // Discretize parabolic Voronoi edge. - // Parabolic Voronoi edges are always formed by one point and one segment - // from the initial input set. - // - // Args: - // point: input point. - // segment: input segment. - // max_dist: maximum discretization distance. - // discretization: point discretization of the given Voronoi edge. - // - // Template arguments: - // InCT: coordinate type of the input geometries (usually integer). - // Point: point type, should model point concept. - // Segment: segment type, should model segment concept. - // - // Important: - // discretization should contain both edge endpoints initially. - template - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - void - >::type discretize( - const Point& point, - const Segment& segment, - const CT max_dist, - std::vector< Point >* discretization) { - // Apply the linear transformation to move start point of the segment to - // the point with coordinates (0, 0) and the direction of the segment to - // coincide the positive direction of the x-axis. - CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y; - - // Compute x-coordinates of the endpoints of the edge - // in the transformed space. - CT projection_start = sqr_segment_length * - get_point_projection((*discretization)[0], segment); - CT projection_end = sqr_segment_length * - get_point_projection((*discretization)[1], segment); - - // Compute parabola parameters in the transformed space. - // Parabola has next representation: - // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y). - CT point_vec_x = cast(x(point)) - cast(x(low(segment))); - CT point_vec_y = cast(y(point)) - cast(y(low(segment))); - CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y; - CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x; - - // Save the last point. - Point last_point = (*discretization)[1]; - discretization->pop_back(); - - // Use stack to avoid recursion. - std::stack point_stack; - point_stack.push(projection_end); - CT cur_x = projection_start; - CT cur_y = parabola_y(cur_x, rot_x, rot_y); - - // Adjust max_dist parameter in the transformed space. - const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length; - while (!point_stack.empty()) { - CT new_x = point_stack.top(); - CT new_y = parabola_y(new_x, rot_x, rot_y); - - // Compute coordinates of the point of the parabola that is - // furthest from the current line segment. - CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x; - CT mid_y = parabola_y(mid_x, rot_x, rot_y); - - // Compute maximum distance between the given parabolic arc - // and line segment that discretize it. - CT dist = (new_y - cur_y) * (mid_x - cur_x) - - (new_x - cur_x) * (mid_y - cur_y); - dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) + - (new_x - cur_x) * (new_x - cur_x)); - if (dist <= max_dist_transformed) { - // Distance between parabola and line segment is less than max_dist. - point_stack.pop(); - CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) / - sqr_segment_length + cast(x(low(segment))); - CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) / - sqr_segment_length + cast(y(low(segment))); - discretization->push_back(Point(inter_x, inter_y)); - cur_x = new_x; - cur_y = new_y; - } else { - point_stack.push(mid_x); - } - } - - // Update last point. - discretization->back() = last_point; - } - - private: - // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b). - static CT parabola_y(CT x, CT a, CT b) { - return ((x - a) * (x - a) + b * b) / (b + b); - } - - // Get normalized length of the distance between: - // 1) point projection onto the segment - // 2) start point of the segment - // Return this length divided by the segment length. This is made to avoid - // sqrt computation during transformation from the initial space to the - // transformed one and vice versa. The assumption is made that projection of - // the point lies between the start-point and endpoint of the segment. - template - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - CT - >::type get_point_projection( - const Point& point, const Segment& segment) { - CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT point_vec_x = x(point) - cast(x(low(segment))); - CT point_vec_y = y(point) - cast(y(low(segment))); - CT sqr_segment_length = - segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y; - CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y; - return vec_dot / sqr_segment_length; - } - - template - static CT cast(const InCT& value) { - return static_cast(value); - } -}; -} -} - -#endif // BOOST_POLYGON_VORONOI_VISUAL_UTILS diff --git a/xs/xsp/ExPolygon.xsp b/xs/xsp/ExPolygon.xsp index b9075f3d..caadb9e5 100644 --- a/xs/xsp/ExPolygon.xsp +++ b/xs/xsp/ExPolygon.xsp @@ -26,7 +26,7 @@ ExPolygons simplify(double tolerance); Polygons simplify_p(double tolerance); Polylines medial_axis(double width) - %code{% THIS->medial_axis(&RETVAL); %}; + %code{% THIS->medial_axis(width, &RETVAL); %}; %{ ExPolygon* diff --git a/xs/xsp/Line.xsp b/xs/xsp/Line.xsp index 8e8ff49c..c82f458d 100644 --- a/xs/xsp/Line.xsp +++ b/xs/xsp/Line.xsp @@ -21,6 +21,8 @@ void scale(double factor); void translate(double x, double y); double length(); + double atan2_(); + double direction(); Point* midpoint() %code{% const char* CLASS = "Slic3r::Point"; RETVAL = THIS->midpoint(); %}; Point* point_at(double distance)