Finished implementing Boost.Polygon medial axis. Some cleanup needed

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

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);
 

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;

xs/src/Geometry.cpp

@@ -1,4 +1,6 @@
 #include "Geometry.hpp"
+#include "Line.hpp"
+#include "PolylineCollection.hpp"
 #include "clipper.hpp"
 #include <algorithm>
 #include <map>
@@ -94,30 +96,109 @@ 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
+    // iterate through the diagram by starting from a random edge
-    // (Boost.Polygon returns it duplicated in both directions)
+    this->edge_cache.clear();
-    std::set<const voronoi_diagram<double>::edge_type*> edge_cache;
+    for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge)
+        this->process_edge(*edge, polylines);
+}
 
-    // iterate through the diagram
+void
-    for (voronoi_diagram<double>::const_edge_iterator it = this->vd.edges().begin(); it != this->vd.edges().end(); ++it) {
+MedialAxis::process_edge(const VD::edge_type& edge, Polylines* polylines)
-        (void)edge_cache.insert(it->twin());
+{
-        if (edge_cache.count(&*it) > 0) continue;
+    // if we already visited this edge or its twin skip it
-        if (!it->is_primary()) continue;
+    if (this->edge_cache.count(&edge) > 0) return;
     
-        Polyline p;
+    // mark this as already visited
-        if (!it->is_finite()) {
+    (void)this->edge_cache.insert(&edge);
-            this->clip_infinite_edge(*it, &p.points);
+    (void)this->edge_cache.insert(edge.twin());
-        } else {
+    
-            p.points.push_back(Point( it->vertex0()->x(), it->vertex0()->y() ));
+    if (this->is_valid_edge(edge)) {
-            p.points.push_back(Point( it->vertex1()->x(), it->vertex1()->y() ));
+        Line line = Line(
-            if (it->is_curved()) {
+            Point( edge.vertex0()->x(), edge.vertex0()->y() ),
-                this->sample_curved_edge(*it, &p.points);
+            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<const VD::edge_type*> 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 VD::edge_type*>::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<double>::cell_type &cell1 = *edge.cell();
+    const voronoi_diagram<double>::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
 MedialAxis::clip_infinite_edge(const voronoi_diagram<double>::edge_type& edge, Points* clipped_edge)
 {
@@ -203,7 +284,7 @@ MedialAxis::retrieve_point(const voronoi_diagram<double>::cell_type& cell)
 }
 
 Line
-MedialAxis::retrieve_segment(const voronoi_diagram<double>::cell_type& cell)
+MedialAxis::retrieve_segment(const voronoi_diagram<double>::cell_type& cell) const
 {
     voronoi_diagram<double>::cell_type::source_index_type index = cell.source_index() - this->points.size();
     return this->lines[index];

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<double>::edge_type& edge, Polylines* polylines);
+    void process_edge_neighbors(const voronoi_diagram<double>::edge_type& edge, Polylines* polylines);
+    bool is_valid_edge(const voronoi_diagram<double>::edge_type& edge) const;
     void clip_infinite_edge(const voronoi_diagram<double>::edge_type& edge, Points* clipped_edge);
     void sample_curved_edge(const voronoi_diagram<double>::edge_type& edge, Points* sampled_edge);
     Point retrieve_point(const voronoi_diagram<double>::cell_type& cell);
-    Line retrieve_segment(const voronoi_diagram<double>::cell_type& cell);
+    Line retrieve_segment(const voronoi_diagram<double>::cell_type& cell) const;
     
     private:
-    voronoi_diagram<double> vd;
+    typedef voronoi_diagram<double> VD;
+    VD vd;
     BoundingBox bb;
+    std::set<const VD::edge_type*> edge_cache;
 };
 
 } }

xs/src/PolylineCollection.hpp

@@ -10,8 +10,8 @@ class PolylineCollection
 {
     public:
     Polylines polylines;
-    PolylineCollection* chained_path(bool no_reverse) const;
+    PolylineCollection* chained_path(bool no_reverse = false) const;
-    PolylineCollection* chained_path_from(const Point* start_near, bool no_reverse) const;
+    PolylineCollection* chained_path_from(const Point* start_near, bool no_reverse = false) const;
     Point* leftmost_point() const;
 };
 

xs/src/TriangleMesh.cpp

@@ -580,10 +580,10 @@ TriangleMesh::slice(const std::vector<double> &z, std::vector<ExPolygons>* 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",
+        printf("Layer %zu (slice_z = %.2f): %zu surface(s) having %zu holes detected from %zu polylines\n",
-            layer_id, z[layer_id], ex_slices.count(), holes_count, loops->count());
+            layer_id, z[layer_id], ex_slices.size(), holes_count, loops->size());
         #endif
         
         ExPolygons* layer = &(*layers)[layer_id];