diff --git a/xs/src/libslic3r/Geometry.cpp b/xs/src/libslic3r/Geometry.cpp index 996d3fe10..e30e73702 100644 --- a/xs/src/libslic3r/Geometry.cpp +++ b/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); diff --git a/xs/src/libslic3r/MultiPoint.cpp b/xs/src/libslic3r/MultiPoint.cpp index 24f210a99..e4944edca 100644 --- a/xs/src/libslic3r/MultiPoint.cpp +++ b/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(); diff --git a/xs/src/libslic3r/Point.cpp b/xs/src/libslic3r/Point.cpp index 77f4224d5..d40a8efc7 100644 --- a/xs/src/libslic3r/Point.cpp +++ b/xs/src/libslic3r/Point.cpp @@ -1,6 +1,7 @@ #include "Point.hpp" #include "Line.hpp" #include "MultiPoint.hpp" +#include #include #include @@ -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); diff --git a/xs/src/libslic3r/Point.hpp b/xs/src/libslic3r/Point.hpp index 368fb1af8..6fa216405 100644 --- a/xs/src/libslic3r/Point.hpp +++ b/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; diff --git a/xs/t/03_point.t b/xs/t/03_point.t index 0ec9df552..d983fafc9 100644 --- a/xs/t/03_point.t +++ b/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'; } { diff --git a/xs/t/09_polyline.t b/xs/t/09_polyline.t index 5ce8d87cb..99077aad3 100644 --- a/xs/t/09_polyline.t +++ b/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'; diff --git a/xs/xsp/Point.xsp b/xs/xsp/Point.xsp index e45c41213..871911239 100644 --- a/xs/xsp/Point.xsp +++ b/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)