Iterative, not recursive, version of the Douglas-Peucker-Ramer algorithm

based on the work by @fuchstraumer
https://github.com/slic3r/Slic3r/pull/3825
https://gist.github.com/fuchstraumer/9421573fc281b946e5f561758961212a
which was based on
http://anis-moussa.blogspot.com/2014/03/ramer-douglas-peucker-algorithm-for.html

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 */
+// Distance to the closest point of line.
-double Line::distance_to(const Point &point) const
+double Line::distance_to_squared(const Point &point, const Point &a, const Point &b)
 {
-    const Line   &line = *this;
+    const Vec2d   v  = (b - a).cast<double>();
-    const Vec2d   v  = (line.b - line.a).cast<double>();
+    const Vec2d   va = (point  - a).cast<double>();
-    const Vec2d   va = (point  - line.a).cast<double>();
     const double  l2 = v.squaredNorm();  // avoid a sqrt
     if (l2 == 0.0) 
-        // line.a == line.b case
+        // a == b case
-        return va.norm();
+        return va.squaredNorm();
-    // Consider the line extending the segment, parameterized as line.a + t (line.b - line.a).
+    // 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
+    if (t < 0.0)      return va.squaredNorm();  // beyond the 'a' end of the segment
-    else if (t > 1.0) return (point - line.b).cast<double>().norm();  // beyond the 'b' end of the segment
+    else if (t > 1.0) return (point - b).cast<double>().squaredNorm();  // beyond the 'b' end of the segment
-    return (t * v - va).norm();
+    return (t * v - va).squaredNorm();
 }
 
 double Line::perp_distance_to(const Point &point) const

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;
 };

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.
+std::vector<Point> MultiPoint::_douglas_peucker(const std::vector<Point>& pts, const double tolerance)
-// The recursive algorithm shall run in place, not allocating temporary data in each recursion.
-Points
-MultiPoint::_douglas_peucker(const Points &points, const double tolerance)
 {
-    assert(points.size() >= 2);
+    std::vector<Point> result_pts;
-    Points results;
+    if (! pts.empty()) {
-    double dmax = 0;
+        const Point  *anchor      = &pts.front();
-    size_t index = 0;
+        size_t        anchor_idx  = 0;
-    Line full(points.front(), points.back());
+        const Point  *floater     = &pts.back();
-    for (Points::const_iterator it = points.begin() + 1; it != points.end(); ++it) {
+        size_t        floater_idx = pts.size() - 1;
-        // we use shortest distance, not perpendicular distance
+        result_pts.reserve(pts.size());
-        double d = full.distance_to(*it);
+        result_pts.emplace_back(*anchor);
-        if (d > dmax) {
+        if (anchor_idx != floater_idx) {
-            index = it - points.begin();
+            assert(pts.size() > 1);
-            dmax = d;
+            std::vector<size_t> 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) {
+    return result_pts;
-        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;
 }
 
 // Visivalingam simplification algorithm https://github.com/slic3r/Slic3r/pull/3825