Bugfix: medial axis missed some segments. #2144

t/angles.t

@@ -2,7 +2,7 @@ use Test::More;
 use strict;
 use warnings;
 
-plan tests => 26;
+plan tests => 34;
 
 BEGIN {
     use FindBin;
@@ -28,6 +28,20 @@ use Slic3r::Geometry qw(rad2deg_dir angle3points PI);
 
 #==========================================================
 
+{
+    is line_orientation([ [0, 0],  [10, 0] ]),  (0),      'E orientation';
+    is line_orientation([ [0, 0],  [0, 10] ]),  (PI/2),   'N orientation';
+    is line_orientation([ [10, 0], [0, 0]  ]),  (PI),     'W orientation';
+    is line_orientation([ [0, 10], [0, 0]  ]),  (PI*3/2), 'S orientation';
+    
+    is line_orientation([ [0, 0], [10, 10] ]),  (PI*1/4), 'NE orientation';
+    is line_orientation([ [10, 0], [0, 10] ]),  (PI*3/4), 'NW orientation';
+    is line_orientation([ [10, 10], [0, 0] ]),  (PI*5/4), 'SW orientation';
+    is line_orientation([ [0, 10], [10, 0] ]),  (PI*7/4), 'SE orientation';
+}
+
+#==========================================================
+
 {
     is line_direction([ [0, 0],  [10, 0] ]), (0),      'E direction';
     is line_direction([ [10, 0], [0, 0]  ]), (0),      'W direction';
@@ -67,6 +81,11 @@ sub line_atan {
     return Slic3r::Line->new(@$l)->atan2_;
 }
 
+sub line_orientation {
+    my ($l) = @_;
+    return Slic3r::Line->new(@$l)->orientation;
+}
+
 sub line_direction {
     my ($l) = @_;
     return Slic3r::Line->new(@$l)->direction;

t/thin.t

@@ -1,4 +1,4 @@
-use Test::More tests => 13;
+use Test::More tests => 14;
 use strict;
 use warnings;
 
@@ -8,7 +8,7 @@ BEGIN {
 }
 
 use Slic3r;
-use List::Util qw(first);
+use List::Util qw(first sum);
 use Slic3r::Geometry qw(epsilon scale unscale scaled_epsilon Y);
 use Slic3r::Test;
 
@@ -121,4 +121,16 @@ if (0) {
     ok $len > 80*2 && $len < 100*2, 'medial axis has reasonable length';
 }
 
+{
+    my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new(
+        [-203064906,-51459966],[-219312231,-51459966],[-219335477,-51459962],[-219376095,-51459962],[-219412047,-51459966],
+        [-219572094,-51459966],[-219624814,-51459962],[-219642183,-51459962],[-219656665,-51459966],[-220815482,-51459966],
+        [-220815482,-37738966],[-221117540,-37738966],[-221117540,-51762024],[-203064906,-51762024],
+    ));
+    my $polylines = $expolygon->medial_axis(819998, 102499.75);
+    
+    my $perimeter = $expolygon->contour->split_at_first_point->length;
+    ok sum(map $_->length, @$polylines) > $perimeter/2/4*3, 'medial axis has a reasonable length';
+}
+
 __END__

xs/src/Geometry.cpp

@@ -270,15 +270,14 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
         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);
+        double angle = fabs(segment2.orientation() - segment1.orientation());
         
         // 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;
-        
+        if (angle < PI*2/3 ) {
+            return false;
+        }
         
         // each vertex is equidistant to both cell segments
         // but such distance might differ between the two vertices;

xs/src/Line.cpp

@@ -100,6 +100,14 @@ Line::atan2_() const
     return atan2(this->b.y - this->a.y, this->b.x - this->a.x);
 }
 
+double
+Line::orientation() const
+{
+    double angle = this->atan2_();
+    if (angle < 0) angle = 2*PI + angle;
+    return angle;
+}
+
 double
 Line::direction() const
 {

xs/src/Line.hpp

@@ -31,6 +31,7 @@ class Line
     bool parallel_to(double angle) const;
     bool parallel_to(const Line &line) const;
     double atan2_() const;
+    double orientation() const;
     double direction() const;
     Vector vector() const;
     

xs/xsp/Line.xsp

@@ -23,6 +23,7 @@
     void translate(double x, double y);
     double length();
     double atan2_();
+    double orientation();
     double direction();
     bool parallel_to(double angle);
     bool parallel_to_line(Line* line)