Better pruning

t/geometry.t

@@ -2,7 +2,7 @@ use Test::More;
 use strict;
 use warnings;
 
-plan tests => 23;
+plan tests => 28;
 
 BEGIN {
     use FindBin;
@@ -167,3 +167,14 @@ my $polygons = [
 }
 
 #==========================================================
+
+{
+    my $line = Slic3r::Line->new([0, 0], [20, 0]);
+    is +Slic3r::Point->new(10, 10)->distance_to_line($line), 10, 'distance_to';
+    is +Slic3r::Point->new(50, 10)->distance_to_line($line), 10, 'distance_to';
+    is +Slic3r::Point->new(0, 0)->distance_to_line($line), 0, 'distance_to';
+    is +Slic3r::Point->new(20, 0)->distance_to_line($line), 0, 'distance_to';
+    is +Slic3r::Point->new(10, 0)->distance_to_line($line), 0, 'distance_to';
+}
+
+#==========================================================

t/thin.t

@@ -11,7 +11,7 @@ use Slic3r;
 use List::Util qw(first);
 use Slic3r::Geometry qw(epsilon scale unscale);
 use Slic3r::Test;
-
+goto TTT;
 {
     my $config = Slic3r::Config->new_from_defaults;
     $config->set('layer_height', 0.2);
@@ -65,7 +65,7 @@ use Slic3r::Test;
         'medial axis loop has reasonable length';
 }
 
-{
+TTT: {
     my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
         [100, 100],
         [120, 100],
@@ -75,6 +75,21 @@ use Slic3r::Test;
     my $res = $expolygon->medial_axis(scale 10);
     is scalar(@$res), 1, 'medial axis of a narrow rectangle is a single line';
     ok unscale($res->[0]->length) >= (200-100 - (120-100)) - epsilon, 'medial axis has reasonable length';
+    
+    $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
+        [100, 100],
+        [120, 100],
+        [120, 200],
+        [105, 200],  # extra point in the short side
+        [100, 200],
+    ));
+    my $res2 = $expolygon->medial_axis(scale 10);
+    use Slic3r::SVG;
+    Slic3r::SVG::output(
+        "thin.svg",
+        expolygons => [$expolygon],
+        polylines => $res2,
+    );
 }
 
 {

xs/src/ExPolygon.cpp

@@ -141,9 +141,13 @@ ExPolygon::medial_axis(double width, Polylines* polylines) const
     Slic3r::Geometry::MedialAxis ma(width);
     
     // populate list of segments for the Voronoi diagram
-    this->contour.lines(&ma.lines);
+    ExPolygons expp;
-    for (Polygons::const_iterator hole = this->holes.begin(); hole != this->holes.end(); ++hole)
+    this->simplify(scale_(0.01), expp);
-        hole->lines(&ma.lines);
+    for (ExPolygons::const_iterator expolygon = expp.begin(); expolygon != expp.end(); ++expolygon) {
+        expolygon->contour.lines(&ma.lines);
+        for (Polygons::const_iterator hole = expolygon->holes.begin(); hole != expolygon->holes.end(); ++hole)
+            hole->lines(&ma.lines);
+    }
     
     // compute the Voronoi diagram
     ma.build(polylines);

xs/src/Geometry.cpp

@@ -245,10 +245,7 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
        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
+       our wanted output. */
-       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 VD::cell_type &cell1 = *edge.cell();
     const VD::cell_type &cell2 = *edge.twin()->cell();
@@ -257,12 +254,17 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
         Line segment2 = this->retrieve_segment(cell2);
         if (segment1.a == segment2.b || segment1.b == segment2.a) return false;
         
-        /*
+        // calculate relative angle between the two boundary segments
         Vector vec1 = segment1.vector();
         Vector vec2 = segment2.vector();
         double angle = atan2(vec1.x*vec2.y - vec1.y*vec2.x, vec1.x*vec2.x + vec1.y*vec2.y);
-        //if (angle > PI/2) return false;
         
+        // 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)
+        // so we allow some tolerance (say, 30°)
+        if (fabs(angle) < PI - PI/3) 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)
@@ -273,12 +275,26 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
         double dist0 = v0.distance_to(segment1);
         double dist1 = v1.distance_to(segment1);
         double diff = fabs(dist1 - dist0);
-        //if (diff > this->edge_to_line(edge).length()/2 && diff > this->width/5) return false;
+        double dist_between_segments1 = segment1.a.distance_to(segment2);
-        
+        double dist_between_segments2 = segment1.b.distance_to(segment2);
+        printf("w = %f, dist0 = %f, dist1 = %f, diff = %f, seglength = %f, edgelen = %f, s2s = %f / %f\n",
+            unscale(this->width),
+            unscale(dist0), unscale(dist1), unscale(diff), unscale(segment1.length()),
+            unscale(this->edge_to_line(edge).length()),
+            unscale(dist_between_segments1), unscale(dist_between_segments2)
+            );
+        if (dist0 < SCALED_EPSILON && dist1 < SCALED_EPSILON) {
+            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 (dist0 > this->width/2) return false;
+        if (dist1 > this->width*2) {
+            printf(" => too fat, skipping\n");
+            //return false;
+        }
         */
+        
     }
     
     return true;