More robust medial axis pruning. #2800

lib/Slic3r/Layer/PerimeterGenerator.pm

@@ -133,12 +133,12 @@ sub process {
                         
                         # the following offset2 ensures almost nothing in @thin_walls is narrower than $min_width
                         # (actually, something larger than that still may exist due to mitering or other causes)
-                        my $min_width = $pwidth / 4;
+                        my $min_width = $ext_pwidth / 4;
                         @thin_walls = @{offset2_ex($diff, -$min_width/2, +$min_width/2)};
-        
+                        
                         # the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop
-                        @thin_walls = grep $_->length > $pwidth*2,
+                        @thin_walls = grep $_->length > $ext_pwidth*2,
-                            map @{$_->medial_axis($pwidth + $pspacing, $min_width)}, @thin_walls;
+                            map @{$_->medial_axis($ext_pwidth + $ext_pspacing, $min_width)}, @thin_walls;
                         Slic3r::debugf "  %d thin walls detected\n", scalar(@thin_walls) if $Slic3r::debug;
         
                         if (0) {

t/thin.t

@@ -1,4 +1,4 @@
-use Test::More tests => 14;
+use Test::More tests => 16;
 use strict;
 use warnings;
 
@@ -133,4 +133,15 @@ if (0) {
     ok sum(map $_->length, @$polylines) > $perimeter/2/4*3, 'medial axis has a reasonable length';
 }
 
+{
+    my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
+        [50, 100],
+        [300, 102],
+        [50, 104],
+    ));
+    my $res = $expolygon->medial_axis(scale 4, scale 0.5);
+    is scalar(@$res), 1, 'medial axis of a narrow triangle is a single line';
+    ok unscale($res->[0]->length) >= (200-100 - (120-100)) - epsilon, 'medial axis has reasonable length';
+}
+
 __END__

xs/src/libslic3r/Geometry.cpp

@@ -449,69 +449,67 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
        "thin" nature of our input, these edges will be very short and not part of
        our wanted output. */
     
+    // retrieve the original line segments which generated the edge we're checking
     const VD::cell_type &cell1 = *edge.cell();
     const VD::cell_type &cell2 = *edge.twin()->cell();
-    if (cell1.contains_segment() && cell2.contains_segment()) {
+    if (!cell1.contains_segment() || !cell2.contains_segment()) return false;
-        Line segment1 = this->retrieve_segment(cell1);
+    const Line &segment1 = this->retrieve_segment(cell1);
-        Line segment2 = this->retrieve_segment(cell2);
+    const Line &segment2 = this->retrieve_segment(cell2);
-        if (segment1.a == segment2.b || segment1.b == segment2.a) return false;
+    
-        
+    // calculate the relative angle between the two boundary segments
-        // calculate relative angle between the two boundary segments
+    double angle = fabs(segment2.orientation() - segment1.orientation());
-        double angle = fabs(segment2.orientation() - segment1.orientation());
+    if (angle > PI) angle -= PI;
-        
+    
-        // fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
+    // fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
-        // we're interested only in segments close to the second case (facing segments)
+    // we're interested only in segments close to the second case (facing segments)
-        // so we allow some tolerance (say, 30°)
+    // so we allow some tolerance (say, 30°)
-        if (angle < PI*2/3 ) {
+    if (angle < PI*2/3) {
-            return false;
+        return false;
-        }
-        
-        // each vertex is equidistant to both cell segments
-        // but such distance might differ between the two vertices;
-        // in this case it means the shape is getting narrow (like a corner)
-        // and we might need to skip the edge since it's not really part of
-        // our skeleton
-        Point v0( edge.vertex0()->x(), edge.vertex0()->y() );
-        Point v1( edge.vertex1()->x(), edge.vertex1()->y() );
-        double dist0 = v0.perp_distance_to(segment1);
-        double dist1 = v1.perp_distance_to(segment1);
-        
-        /*
-        double diff = fabs(dist1 - dist0);
-        double dist_between_segments1 = segment1.a.distance_to(segment2);
-        double dist_between_segments2 = segment1.b.distance_to(segment2);
-        printf("w = %f/%f, dist0 = %f, dist1 = %f, diff = %f, seg1len = %f, seg2len = %f, edgelen = %f, s2s = %f / %f\n",
-            unscale(this->max_width), unscale(this->min_width),
-            unscale(dist0), unscale(dist1), unscale(diff),
-            unscale(segment1.length()), unscale(segment2.length()),
-            unscale(this->edge_to_line(edge).length()),
-            unscale(dist_between_segments1), unscale(dist_between_segments2)
-            );
-        */
-        
-        // if this segment is the centerline for a very thin area, we might want to skip it
-        // in case the area is too thin
-        if (dist0 < this->min_width/2 || dist1 < this->min_width/2) {
-            //printf(" => too thin, skipping\n");
-            return false;
-        }
-        
-        /*
-        // if distance between this edge and the thin area boundary is greater
-        // than half the max width, then it's not a true medial axis segment
-        if (dist1 > this->width*2) {
-            printf(" => too fat, skipping\n");
-            //return false;
-        }
-        */
-        
-        return true;
     }
     
-    return false;
+    // each edge vertex is equidistant to both cell segments
+    // but such distance might differ between the two vertices;
+    // in this case it means the shape is getting narrow (like a corner)
+    // and we might need to skip the edge since it's not really part of
+    // our skeleton
+    
+    // get perpendicular distance of each edge vertex to the segment(s)
+    Line line = this->edge_to_line(edge);
+    double dist0 = line.a.perp_distance_to(segment1);
+    double dist1 = line.b.perp_distance_to(segment1);
+    
+    /*
+    double diff = fabs(dist1 - dist0);
+    double dist_between_segments1 = segment1.a.distance_to(segment2);
+    double dist_between_segments2 = segment1.b.distance_to(segment2);
+    printf("w = %f/%f, dist0 = %f, dist1 = %f, diff = %f, seg1len = %f, seg2len = %f, edgelen = %f, s2s = %f / %f\n",
+        unscale(this->max_width), unscale(this->min_width),
+        unscale(dist0), unscale(dist1), unscale(diff),
+        unscale(segment1.length()), unscale(segment2.length()),
+        unscale(line.length()),
+        unscale(dist_between_segments1), unscale(dist_between_segments2)
+        );
+    */
+
+    // if this edge is the centerline for a very thin area, we might want to skip it
+    // in case the area is too thin
+    if (dist0 < this->min_width/2 && dist1 < this->min_width/2) {
+        //printf(" => too thin, skipping\n");
+        return false;
+    }
+    
+    // if only one of the two edge vertices is the centerline of a very narrow area,
+    // it means the shape is shrinking on that size (dist0 or dist1 might be even 0
+    // when we have a narrow triangle), so in this case we only keep the edge if it's
+    // long enough, otherwise it's an artifact
+    if (dist0 < this->min_width/2 || dist1 < this->min_width/2) {
+        if (line.length() < this->max_width) return false;
+    }
+
+    return true;
 }
 
-Line
+const Line&
 MedialAxis::retrieve_segment(const VD::cell_type& cell) const
 {
     VD::cell_type::source_index_type index = cell.source_index() - this->points.size();

xs/src/libslic3r/Geometry.hpp

@@ -54,7 +54,7 @@ class MedialAxis {
     Line edge_to_line(const VD::edge_type &edge) const;
     void process_edge_neighbors(const voronoi_diagram<double>::edge_type& edge, Points* points);
     bool is_valid_edge(const voronoi_diagram<double>::edge_type& edge) const;
-    Line retrieve_segment(const voronoi_diagram<double>::cell_type& cell) const;
+    const Line& retrieve_segment(const voronoi_diagram<double>::cell_type& cell) const;
 };
 
 } }

xs/t/03_point.t

@@ -4,7 +4,7 @@ use strict;
 use warnings;
 
 use Slic3r::XS;
-use Test::More tests => 22;
+use Test::More tests => 24;
 
 my $point = Slic3r::Point->new(10, 15);
 is_deeply [ @$point ], [10, 15], 'point roundtrip';
@@ -32,12 +32,14 @@ ok !$point->coincides_with($point2), 'coincides_with';
 
 {
     my $line = Slic3r::Line->new([0,0], [100,0]);
-    is +Slic3r::Point->new(0,0)->distance_to_line($line),    0, 'distance_to_line()';
+    is +Slic3r::Point->new(0,0)  ->distance_to_line($line),  0, 'distance_to_line()';
     is +Slic3r::Point->new(100,0)->distance_to_line($line),  0, 'distance_to_line()';
-    is +Slic3r::Point->new(50,0)->distance_to_line($line),   0, 'distance_to_line()';
+    is +Slic3r::Point->new(50,0) ->distance_to_line($line),  0, 'distance_to_line()';
     is +Slic3r::Point->new(150,0)->distance_to_line($line), 50, 'distance_to_line()';
-    is +Slic3r::Point->new(0,50)->distance_to_line($line),  50, 'distance_to_line()';
+    is +Slic3r::Point->new(0,50) ->distance_to_line($line), 50, 'distance_to_line()';
     is +Slic3r::Point->new(50,50)->distance_to_line($line), 50, 'distance_to_line()';
+    is +Slic3r::Point->new(50,50) ->perp_distance_to_line($line), 50, 'perp_distance_to_line()';
+    is +Slic3r::Point->new(150,50)->perp_distance_to_line($line), 50, 'perp_distance_to_line()';
 }
 
 {