From 7b0decbeb11c53324333b03bc459c657464cb963 Mon Sep 17 00:00:00 2001 From: Alessandro Ranellucci Date: Fri, 10 Jan 2014 16:18:55 +0100 Subject: [PATCH] Finished implementing Boost.Polygon medial axis. Some cleanup needed --- Build.PL | 1 - lib/Slic3r/ExPolygon.pm | 164 ---------------------------------- lib/Slic3r/Layer/Region.pm | 7 +- xs/src/Geometry.cpp | 119 ++++++++++++++++++++---- xs/src/Geometry.hpp | 10 ++- xs/src/PolylineCollection.hpp | 4 +- xs/src/TriangleMesh.cpp | 6 +- 7 files changed, 117 insertions(+), 194 deletions(-) diff --git a/Build.PL b/Build.PL index e5ad62cc..26c6f6f9 100644 --- a/Build.PL +++ b/Build.PL @@ -13,7 +13,6 @@ my %prereqs = qw( File::Basename 0 File::Spec 0 Getopt::Long 0 - Math::Geometry::Voronoi 1.3 Math::PlanePath 53 Module::Build::WithXSpp 0.14 Moo 1.003001 diff --git a/lib/Slic3r/ExPolygon.pm b/lib/Slic3r/ExPolygon.pm index 584eb952..ec3d92f5 100644 --- a/lib/Slic3r/ExPolygon.pm +++ b/lib/Slic3r/ExPolygon.pm @@ -5,7 +5,6 @@ use warnings; # an ExPolygon is a polygon with holes use List::Util qw(first); -use Math::Geometry::Voronoi; use Slic3r::Geometry qw(X Y A B point_in_polygon epsilon scaled_epsilon); use Slic3r::Geometry::Clipper qw(union_ex diff_pl); @@ -43,169 +42,6 @@ sub bounding_box { return $self->contour->bounding_box; } -# this method only works for expolygons having only a contour or -# a contour and a hole, and not being thicker than the supplied -# width. it returns a polyline or a polygon -sub ___medial_axis { - my ($self, $width) = @_; - return $self->_medial_axis_voronoi($width); -} - -sub _medial_axis_clip { - my ($self, $width) = @_; - - my $grow = sub { - my ($line, $distance) = @_; - - my $line_clone = $line->clone; - $line_clone->clip_start(scaled_epsilon); - return () if !$line_clone->is_valid; - $line_clone->clip_end(scaled_epsilon); - return () if !$line_clone->is_valid; - - my ($a, $b) = @$line_clone; - my $dx = $a->x - $b->x; - my $dy = $a->y - $b->y; #- - my $dist = sqrt($dx*$dx + $dy*$dy); - $dx /= $dist; - $dy /= $dist; - return Slic3r::Polygon->new( - Slic3r::Point->new($a->x + $distance*$dy, $a->y - $distance*$dx), #-- - Slic3r::Point->new($b->x + $distance*$dy, $b->y - $distance*$dx), #-- - Slic3r::Point->new($b->x - $distance*$dy, $b->y + $distance*$dx), #++ - Slic3r::Point->new($a->x - $distance*$dy, $a->y + $distance*$dx), #++ - ); - }; - - my @result = (); - my $covered = []; - foreach my $polygon (@$self) { - my @polylines = (); - foreach my $line (@{$polygon->lines}) { - # remove the areas that are already covered from this line - my $clipped = diff_pl([$line->as_polyline], $covered); - - # skip very short segments/dots - @$clipped = grep $_->length > $width/10, @$clipped; - - # grow the remaining lines and add them to the covered areas - push @$covered, map $grow->($_, $width*1.1), @$clipped; - - # if the first remaining segment is connected to the last polyline, append it - # to that -- FIXME: this assumes that diff_pl() - # preserved the orientation of the input linestring but this is not generally true - if (@polylines && @$clipped && $clipped->[0]->first_point->distance_to($polylines[-1]->last_point) <= $width/10) { - $polylines[-1]->append_polyline(shift @$clipped); - } - push @polylines, @$clipped; - } - - foreach my $polyline (@polylines) { - # if this polyline looks like a closed loop, return it as a polygon - if ($polyline->first_point->coincides_with($polyline->last_point)) { - next if @$polyline == 2; - $polyline->pop_back; - push @result, Slic3r::Polygon->new(@$polyline); - } else { - push @result, $polyline; - } - } - } - - return @result; -} - -my $voronoi_lock :shared; - -sub _medial_axis_voronoi { - my ($self, $width) = @_; - - lock($voronoi_lock); - - my $voronoi; - { - my @points = (); - foreach my $polygon (@$self) { - { - my $p = $polygon->pp; - Slic3r::Geometry::polyline_remove_short_segments($p, $width / 2); - $polygon = Slic3r::Polygon->new(@$p); - } - - # subdivide polygon segments so that we don't have anyone of them - # being longer than $width / 2 - $polygon = $polygon->subdivide($width/2); - - push @points, @{$polygon->pp}; - } - $voronoi = Math::Geometry::Voronoi->new(points => \@points); - } - - $voronoi->compute; - my $vertices = $voronoi->vertices; - - my @skeleton_lines = (); - foreach my $edge (@{ $voronoi->edges }) { - # ignore lines going to infinite - next if $edge->[1] == -1 || $edge->[2] == -1; - - my $line = Slic3r::Line->new($vertices->[$edge->[1]], $vertices->[$edge->[2]]); - next if !$self->contains_line($line); - - # contains_point() could be faster, but we need an implementation that - # reliably considers points on boundary - #next if !$self->contains_point(Slic3r::Point->new(@{$vertices->[$edge->[1]]})) - # || !$self->contains_point(Slic3r::Point->new(@{$vertices->[$edge->[2]]})); - - push @skeleton_lines, [$edge->[1], $edge->[2]]; - } - return () if !@skeleton_lines; - - # now walk along the medial axis and build continuos polylines or polygons - my @polylines = (); - { - my @lines = @skeleton_lines; - push @polylines, [ map @$_, shift @lines ]; - CYCLE: while (@lines) { - for my $i (0..$#lines) { - if ($lines[$i][0] == $polylines[-1][-1]) { - push @{$polylines[-1]}, $lines[$i][1]; - } elsif ($lines[$i][1] == $polylines[-1][-1]) { - push @{$polylines[-1]}, $lines[$i][0]; - } elsif ($lines[$i][1] == $polylines[-1][0]) { - unshift @{$polylines[-1]}, $lines[$i][0]; - } elsif ($lines[$i][0] == $polylines[-1][0]) { - unshift @{$polylines[-1]}, $lines[$i][1]; - } else { - next; - } - splice @lines, $i, 1; - next CYCLE; - } - push @polylines, [ map @$_, shift @lines ]; - } - } - - my @result = (); - my $simplify_tolerance = $width / 7; - foreach my $polyline (@polylines) { - next unless @$polyline >= 2; - - # now replace point indexes with coordinates - my @points = map Slic3r::Point->new(@{$vertices->[$_]}), @$polyline; - - if ($points[0]->coincides_with($points[-1])) { - next if @points == 2; - push @result, @{Slic3r::Polygon->new(@points[0..$#points-1])->simplify($simplify_tolerance)}; - } else { - push @result, Slic3r::Polyline->new(@points); - $result[-1]->simplify($simplify_tolerance); - } - } - - return @result; -} - package Slic3r::ExPolygon::Collection; use Slic3r::Geometry qw(X1 Y1); diff --git a/lib/Slic3r/Layer/Region.pm b/lib/Slic3r/Layer/Region.pm index 970dd534..cce81e69 100644 --- a/lib/Slic3r/Layer/Region.pm +++ b/lib/Slic3r/Layer/Region.pm @@ -224,16 +224,17 @@ sub make_perimeters { # process thin walls by collapsing slices to single passes if (@thin_walls) { my @p = map @{$_->medial_axis($pspacing)}, @thin_walls; + if (0) { - use Slic3r::SVG; + require "Slic3r/SVG.pm"; Slic3r::SVG::output( "medial_axis.svg", no_arrows => 1, - expolygons => \@thin_walls, + #expolygons => \@thin_walls, polylines => \@p, ); - exit; } + my @paths = (); for my $p (@p) { next if $p->length <= $pspacing * 2; diff --git a/xs/src/Geometry.cpp b/xs/src/Geometry.cpp index a82680a1..bd17b157 100644 --- a/xs/src/Geometry.cpp +++ b/xs/src/Geometry.cpp @@ -1,4 +1,6 @@ #include "Geometry.hpp" +#include "Line.hpp" +#include "PolylineCollection.hpp" #include "clipper.hpp" #include #include @@ -94,28 +96,107 @@ MedialAxis::build(Polylines* polylines) construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd); - // prepare a cache of twin edges to prevent getting the same edge twice - // (Boost.Polygon returns it duplicated in both directions) - std::set::edge_type*> edge_cache; + // iterate through the diagram by starting from a random edge + this->edge_cache.clear(); + for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) + this->process_edge(*edge, polylines); +} + +void +MedialAxis::process_edge(const VD::edge_type& edge, Polylines* polylines) +{ + // if we already visited this edge or its twin skip it + if (this->edge_cache.count(&edge) > 0) return; - // iterate through the diagram - for (voronoi_diagram::const_edge_iterator it = this->vd.edges().begin(); it != this->vd.edges().end(); ++it) { - (void)edge_cache.insert(it->twin()); - if (edge_cache.count(&*it) > 0) continue; - if (!it->is_primary()) continue; - - Polyline p; - if (!it->is_finite()) { - this->clip_infinite_edge(*it, &p.points); - } else { - p.points.push_back(Point( it->vertex0()->x(), it->vertex0()->y() )); - p.points.push_back(Point( it->vertex1()->x(), it->vertex1()->y() )); - if (it->is_curved()) { - this->sample_curved_edge(*it, &p.points); + // mark this as already visited + (void)this->edge_cache.insert(&edge); + (void)this->edge_cache.insert(edge.twin()); + + if (this->is_valid_edge(edge)) { + Line line = Line( + Point( edge.vertex0()->x(), edge.vertex0()->y() ), + Point( edge.vertex1()->x(), edge.vertex1()->y() ) + ); + bool appended = false; + if (!polylines->empty()) { + Polyline &last_p = polylines->back(); + if (line.a == last_p.points.back()) { + // if this line starts where last polyline ends, just append the other point + last_p.points.push_back(line.b); + appended = true; + } else if (line.b == last_p.points.back()) { + // if this line ends where last polyline ends, just append the other point + last_p.points.push_back(line.a); + appended = true; } } - polylines->push_back(p); + if (polylines->empty() || !appended) { + // start a new polyline + polylines->push_back(Polyline()); + Polyline &p = polylines->back(); + p.points.push_back(line.a); + p.points.push_back(line.b); + } } + + // look for connected edges (on both sides) + this->process_edge_neighbors(edge, polylines); + this->process_edge_neighbors(*edge.twin(), polylines); +} + +void +MedialAxis::process_edge_neighbors(const VD::edge_type& edge, Polylines* polylines) +{ + std::vector neighbors; + for (const VD::edge_type* neighbor = edge.rot_next(); neighbor != &edge; neighbor = neighbor->rot_next()) { + // skip already seen edges + if (this->edge_cache.count(neighbor) > 0) continue; + + // skip edges that we wouldn't include in the MAT anyway + if (!this->is_valid_edge(*neighbor)) continue; + + neighbors.push_back(neighbor); + } + + // process neighbors recursively + if (neighbors.size() == 1) { + this->process_edge(*neighbors.front(), polylines); + } else if (neighbors.size() > 1) { + // close current polyline and start a new one for each branch + for (std::vector::const_iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); ++neighbor) { + Polylines pp; + this->process_edge(**neighbor, &pp); + polylines->insert(polylines->end(), pp.begin(), pp.end()); + } + } +} + +bool +MedialAxis::is_valid_edge(const VD::edge_type& edge) const +{ + // if we only process segments representing closed loops, none if the + // infinite edges (if any) would be part of our MAT anyway + if (edge.is_secondary() || edge.is_infinite()) return false; + + /* If the cells sharing this edge have a common vertex, we're not interested + in this edge. Why? Because it means that the edge lies on the bisector of + two contiguous input lines and it was included in the Voronoi graph because + it's the locus of centers of circles tangent to both vertices. Due to the + "thin" nature of our input, these edges will be very short and not part of + our wanted output. The best way would be to just filter out the edges that + are not the locus of the maximally inscribed disks (requirement of MAT) + 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(); + 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; + } + + return true; } void @@ -203,7 +284,7 @@ MedialAxis::retrieve_point(const voronoi_diagram::cell_type& cell) } Line -MedialAxis::retrieve_segment(const voronoi_diagram::cell_type& cell) +MedialAxis::retrieve_segment(const voronoi_diagram::cell_type& cell) const { voronoi_diagram::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 fba942be..c77fb40d 100644 --- a/xs/src/Geometry.hpp +++ b/xs/src/Geometry.hpp @@ -3,6 +3,7 @@ #include "BoundingBox.hpp" #include "Polygon.hpp" +#include "Polyline.hpp" #include "boost/polygon/voronoi.hpp" using boost::polygon::voronoi_builder; @@ -20,14 +21,19 @@ class MedialAxis { Points points; Lines lines; void build(Polylines* polylines); + void process_edge(const voronoi_diagram::edge_type& edge, Polylines* polylines); + void process_edge_neighbors(const voronoi_diagram::edge_type& edge, Polylines* polylines); + 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); + Line retrieve_segment(const voronoi_diagram::cell_type& cell) const; private: - voronoi_diagram vd; + typedef voronoi_diagram VD; + VD vd; BoundingBox bb; + std::set edge_cache; }; } } diff --git a/xs/src/PolylineCollection.hpp b/xs/src/PolylineCollection.hpp index 2528d0d9..dd2bed71 100644 --- a/xs/src/PolylineCollection.hpp +++ b/xs/src/PolylineCollection.hpp @@ -10,8 +10,8 @@ class PolylineCollection { public: Polylines polylines; - PolylineCollection* chained_path(bool no_reverse) const; - PolylineCollection* chained_path_from(const Point* start_near, bool no_reverse) const; + PolylineCollection* chained_path(bool no_reverse = false) const; + PolylineCollection* chained_path_from(const Point* start_near, bool no_reverse = false) const; Point* leftmost_point() const; }; diff --git a/xs/src/TriangleMesh.cpp b/xs/src/TriangleMesh.cpp index c5f590da..90da0571 100644 --- a/xs/src/TriangleMesh.cpp +++ b/xs/src/TriangleMesh.cpp @@ -580,10 +580,10 @@ TriangleMesh::slice(const std::vector &z, std::vector* layer #ifdef SLIC3R_DEBUG size_t holes_count = 0; for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++e) { - holes_count += e->holes.count(); + holes_count += e->holes.size(); } - printf("Layer %d (slice_z = %.2f): %d surface(s) having %d holes detected from %d polylines\n", - layer_id, z[layer_id], ex_slices.count(), holes_count, loops->count()); + printf("Layer %zu (slice_z = %.2f): %zu surface(s) having %zu holes detected from %zu polylines\n", + layer_id, z[layer_id], ex_slices.size(), holes_count, loops->size()); #endif ExPolygons* layer = &(*layers)[layer_id];