Bugfix: Douglas-Peucker used perpendicular distance instead of shortest distance, thus clipping more than it should. #2474

xs/src/libslic3r/Geometry.cpp

@@ -328,8 +328,8 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
         // our skeleton
         Point v0( edge.vertex0()->x(), edge.vertex0()->y() );
         Point v1( edge.vertex1()->x(), edge.vertex1()->y() );
-        double dist0 = v0.distance_to(segment1);
-        double dist1 = v1.distance_to(segment1);
+        double dist0 = v0.perp_distance_to(segment1);
+        double dist1 = v1.perp_distance_to(segment1);
         
         /*
         double diff = fabs(dist1 - dist0);

xs/src/libslic3r/MultiPoint.cpp

@@ -90,6 +90,7 @@ MultiPoint::_douglas_peucker(const Points &points, const double tolerance)
     size_t index = 0;
     Line full(points.front(), points.back());
     for (Points::const_iterator it = points.begin() + 1; it != points.end(); ++it) {
+        // we use shortest distance, not perpendicular distance
         double d = it->distance_to(full);
         if (d > dmax) {
             index = it - points.begin();

xs/src/libslic3r/Point.cpp

@@ -1,6 +1,7 @@
 #include "Point.hpp"
 #include "Line.hpp"
 #include "MultiPoint.hpp"
+#include <algorithm>
 #include <cmath>
 #include <sstream>
 
@@ -130,8 +131,31 @@ Point::distance_to(const Point &point) const
     return sqrt(dx*dx + dy*dy);
 }
 
+/* distance to the closest point of line */
 double
 Point::distance_to(const Line &line) const
+{
+    const double dx = line.b.x - line.a.x;
+    const double dy = line.b.y - line.a.y;
+    
+    const double l2 = dx*dx + dy*dy;  // avoid a sqrt
+    if (l2 == 0.0) return this->distance_to(line.a);   // line.a == line.b case
+    
+    // Consider the line extending the segment, parameterized as line.a + t (line.b - line.a).
+    // We find projection of this point onto the line. 
+    // It falls where t = [(this-line.a) . (line.b-line.a)] / |line.b-line.a|^2
+    const double t = ((this->x - line.a.x) * dx + (this->y - line.a.y) * dy) / l2;
+    if (t < 0.0)      return this->distance_to(line.a);  // beyond the 'a' end of the segment
+    else if (t > 1.0) return this->distance_to(line.b);  // beyond the 'b' end of the segment
+    Point projection(
+        line.a.x + t * dx,
+        line.a.y + t * dy
+    );
+    return this->distance_to(projection);
+}
+
+double
+Point::perp_distance_to(const Line &line) const
 {
     if (line.a.coincides_with(line.b)) return this->distance_to(line.a);
     

xs/src/libslic3r/Point.hpp

@@ -48,6 +48,7 @@ class Point
     bool nearest_point(const Points &points, Point* point) const;
     double distance_to(const Point &point) const;
     double distance_to(const Line &line) const;
+    double perp_distance_to(const Line &line) const;
     double ccw(const Point &p1, const Point &p2) const;
     double ccw(const Line &line) const;
     double ccw_angle(const Point &p1, const Point &p2) const;

xs/t/03_point.t

@@ -4,7 +4,7 @@ use strict;
 use warnings;
 
 use Slic3r::XS;
-use Test::More tests => 15;
+use Test::More tests => 22;
 
 my $point = Slic3r::Point->new(10, 15);
 is_deeply [ @$point ], [10, 15], 'point roundtrip';
@@ -30,13 +30,28 @@ ok !$point->coincides_with($point2), 'coincides_with';
     ok $nearest->coincides_with($point2), 'nearest_point';
 }
 
+{
+    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(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(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(50,50)->distance_to_line($line), 50, 'distance_to_line()';
+}
+
+{
+    my $line = Slic3r::Line->new([50,50], [125,-25]);
+    is +Slic3r::Point->new(100,0)->distance_to_line($line),  0, 'distance_to_line()';
+}
+
 {
     my $line = Slic3r::Line->new(
         [18335846,18335845],
         [18335846,1664160],
     );
     $point = Slic3r::Point->new(1664161,18335848);
-    is $point->distance_to_line($line), 16671685, 'distance_to_line() does not overflow';
+    is $point->perp_distance_to_line($line), 16671685, 'perp_distance_to_line() does not overflow';
 }
 
 {

xs/t/09_polyline.t

@@ -4,7 +4,7 @@ use strict;
 use warnings;
 
 use Slic3r::XS;
-use Test::More tests => 14;
+use Test::More tests => 16;
 
 my $points = [
     [100, 100],
@@ -34,6 +34,14 @@ is_deeply $polyline->pp, [ @$points, @$points ], 'append_polyline';
     ok abs($polyline->length - ($len-($len/3))) < 1, 'clip_end';
 }
 
+{
+    my $polyline = Slic3r::Polyline->new(
+        [0,0], [20,0], [50,0], [80,0], [100,0],
+    );
+    $polyline->simplify(2);
+    is_deeply $polyline->pp, [ [0,0], [100,0] ], 'Douglas-Peucker';
+}
+
 {
     my $polyline = Slic3r::Polyline->new(
         [0,0], [50,50], [100,0], [125,-25], [150,50],
@@ -42,6 +50,14 @@ is_deeply $polyline->pp, [ @$points, @$points ], 'append_polyline';
     is_deeply $polyline->pp, [ [0, 0], [50, 50], [125, -25], [150, 50] ], 'Douglas-Peucker';
 }
 
+{
+    my $polyline = Slic3r::Polyline->new(
+        [0,0], [100,0], [50,10],
+    );
+    $polyline->simplify(25);
+    is_deeply $polyline->pp, [ [0,0], [100,0], [50,10] ], 'Douglas-Peucker uses shortest distance instead of perpendicular distance';
+}
+
 {
     my $polyline = Slic3r::Polyline->new(@$points);
     is $polyline->length, 100*2, 'length';

xs/xsp/Point.xsp

@@ -29,6 +29,8 @@
         %code{% RETVAL = THIS->distance_to(*point); %};
     double distance_to_line(Line* line)
         %code{% RETVAL = THIS->distance_to(*line); %};
+    double perp_distance_to_line(Line* line)
+        %code{% RETVAL = THIS->perp_distance_to(*line); %};
     double ccw(Point* p1, Point* p2)
         %code{% RETVAL = THIS->ccw(*p1, *p2); %};
     double ccw_angle(Point* p1, Point* p2)