diff --git a/src/libslic3r/Line.cpp b/src/libslic3r/Line.cpp index 35cfa2b76..02f1cb7c2 100644 --- a/src/libslic3r/Line.cpp +++ b/src/libslic3r/Line.cpp @@ -34,23 +34,22 @@ bool Line::intersection_infinite(const Line &other, Point* point) const return true; } -/* distance to the closest point of line */ -double Line::distance_to(const Point &point) const +// Distance to the closest point of line. +double Line::distance_to_squared(const Point &point, const Point &a, const Point &b) { - const Line &line = *this; - const Vec2d v = (line.b - line.a).cast(); - const Vec2d va = (point - line.a).cast(); + const Vec2d v = (b - a).cast(); + const Vec2d va = (point - a).cast(); const double l2 = v.squaredNorm(); // avoid a sqrt if (l2 == 0.0) - // line.a == line.b case - return va.norm(); - // Consider the line extending the segment, parameterized as line.a + t (line.b - line.a). + // a == b case + return va.squaredNorm(); + // Consider the line extending the segment, parameterized as a + t (b - 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 + // It falls where t = [(this-a) . (b-a)] / |b-a|^2 const double t = va.dot(v) / l2; - if (t < 0.0) return va.norm(); // beyond the 'a' end of the segment - else if (t > 1.0) return (point - line.b).cast().norm(); // beyond the 'b' end of the segment - return (t * v - va).norm(); + if (t < 0.0) return va.squaredNorm(); // beyond the 'a' end of the segment + else if (t > 1.0) return (point - b).cast().squaredNorm(); // beyond the 'b' end of the segment + return (t * v - va).squaredNorm(); } double Line::perp_distance_to(const Point &point) const diff --git a/src/libslic3r/Line.hpp b/src/libslic3r/Line.hpp index 36e02247c..559ca946a 100644 --- a/src/libslic3r/Line.hpp +++ b/src/libslic3r/Line.hpp @@ -31,7 +31,8 @@ public: Point midpoint() const { return (this->a + this->b) / 2; } bool intersection_infinite(const Line &other, Point* point) const; bool operator==(const Line &rhs) const { return this->a == rhs.a && this->b == rhs.b; } - double distance_to(const Point &point) const; + double distance_to_squared(const Point &point) const { return distance_to_squared(point, this->a, this->b); } + double distance_to(const Point &point) const { return distance_to(point, this->a, this->b); } double perp_distance_to(const Point &point) const; bool parallel_to(double angle) const; bool parallel_to(const Line &line) const { return this->parallel_to(line.direction()); } @@ -43,6 +44,9 @@ public: bool intersection(const Line& line, Point* intersection) const; double ccw(const Point& point) const { return point.ccw(*this); } + static double distance_to_squared(const Point &point, const Point &a, const Point &b); + static double distance_to(const Point &point, const Point &a, const Point &b) { return sqrt(distance_to_squared(point, a, b)); } + Point a; Point b; }; diff --git a/src/libslic3r/MultiPoint.cpp b/src/libslic3r/MultiPoint.cpp index 93c63c4cf..795ee38c2 100644 --- a/src/libslic3r/MultiPoint.cpp +++ b/src/libslic3r/MultiPoint.cpp @@ -162,45 +162,51 @@ bool MultiPoint::first_intersection(const Line& line, Point* intersection) const return found; } -//FIXME This is very inefficient in term of memory use. -// The recursive algorithm shall run in place, not allocating temporary data in each recursion. -Points -MultiPoint::_douglas_peucker(const Points &points, const double tolerance) +std::vector MultiPoint::_douglas_peucker(const std::vector& pts, const double tolerance) { - assert(points.size() >= 2); - Points results; - double dmax = 0; - 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 = full.distance_to(*it); - if (d > dmax) { - index = it - points.begin(); - dmax = d; + std::vector result_pts; + if (! pts.empty()) { + const Point *anchor = &pts.front(); + size_t anchor_idx = 0; + const Point *floater = &pts.back(); + size_t floater_idx = pts.size() - 1; + result_pts.reserve(pts.size()); + result_pts.emplace_back(*anchor); + if (anchor_idx != floater_idx) { + assert(pts.size() > 1); + std::vector dpStack; + dpStack.reserve(pts.size()); + dpStack.emplace_back(floater_idx); + for (;;) { + double max_distSq = 0.0; + size_t furthest_idx = anchor_idx; + // find point furthest from line seg created by (anchor, floater) and note it + for (size_t i = anchor_idx + 1; i < floater_idx; ++ i) { + double dist = Line::distance_to_squared(pts[i], *anchor, *floater); + if (dist > max_distSq) { + max_distSq = dist; + furthest_idx = i; + } + } + // remove point if less than tolerance + if (max_distSq <= tolerance) { + result_pts.emplace_back(*floater); + anchor_idx = floater_idx; + anchor = floater; + assert(dpStack.back() == floater_idx); + dpStack.pop_back(); + if (dpStack.empty()) + break; + floater_idx = dpStack.back(); + } else { + floater_idx = furthest_idx; + dpStack.emplace_back(floater_idx); + } + floater = &pts[floater_idx]; + } } } - if (dmax >= tolerance) { - Points dp0; - dp0.reserve(index + 1); - dp0.insert(dp0.end(), points.begin(), points.begin() + index + 1); - // Recursive call. - Points dp1 = MultiPoint::_douglas_peucker(dp0, tolerance); - results.reserve(results.size() + dp1.size() - 1); - results.insert(results.end(), dp1.begin(), dp1.end() - 1); - - dp0.clear(); - dp0.reserve(points.size() - index); - dp0.insert(dp0.end(), points.begin() + index, points.end()); - // Recursive call. - dp1 = MultiPoint::_douglas_peucker(dp0, tolerance); - results.reserve(results.size() + dp1.size()); - results.insert(results.end(), dp1.begin(), dp1.end()); - } else { - results.push_back(points.front()); - results.push_back(points.back()); - } - return results; + return result_pts; } // Visivalingam simplification algorithm https://github.com/slic3r/Slic3r/pull/3825