From 6b3815a7666d49fe108ce12b3eadfcb085367f50 Mon Sep 17 00:00:00 2001
From: Gabriel Wicke <gswicke@gmail.com>
Date: Wed, 17 Jun 2020 20:10:18 -0700
Subject: [PATCH] [path] Implement Ramer-Douglas-Peucker line simplification

Implement an iterative version of the Ramer-Douglas-Peucker line
simplification algorithm
(https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm),
which reduces line complexity to a limited linear deviation from the
original polyline. The ability to reason about linear deflection is the
key improvement over the previous linear implementation.

Worst case complexity is O(n^2), but expected complexity for typical
cases is O(n log n). A potentially faster alternative would be to call
out to libclipper, treating the line as a closed polygon. However, in
practice, performance of this implementation seems good enough. A
complex 3d surface operation optimizes in a few seconds, and reduces
output gcode size from about 220MB with the previous implementation to
10MB.
---
 src/Mod/Path/PathScripts/PathSurface.py | 21 +++-----------
 src/Mod/Path/PathScripts/PathUtils.py   | 38 +++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 17 deletions(-)

diff --git a/src/Mod/Path/PathScripts/PathSurface.py b/src/Mod/Path/PathScripts/PathSurface.py
index 5e05ee625288..44debb89b4a7 100644
--- a/src/Mod/Path/PathScripts/PathSurface.py
+++ b/src/Mod/Path/PathScripts/PathSurface.py
@@ -1078,7 +1078,8 @@ def _planarDropCutSingle(self, JOB, obj, pdc, safePDC, depthparams, SCANDATA):
 
     def _planarSinglepassProcess(self, obj, points):
         if obj.OptimizeLinearPaths:
-            points = self._optimizeLinearSegments(points)
+            points = PathUtils.simplify3dLine(points,
+                    tolerance=obj.LinearDeflection.Value)
         # Begin processing ocl points list into gcode
         commands = []
         for pnt in points:
@@ -2085,28 +2086,14 @@ def setOclCutter(self, obj, safe=False):
             PathLog.warning("Defaulting cutter to standard end mill.")
             return ocl.CylCutter(diam_1, (CEH + lenOfst))
 
-    def _optimizeLinearSegments(self, line):
-        """Eliminate collinear interior segments"""
-        if len(line) > 2:
-            prv, pnt = line[0:2]
-            pts = [prv]
-            for nxt in line[2:]:
-                if not pnt.isOnLineSegment(prv, nxt):
-                    pts.append(pnt)
-                    prv = pnt
-                pnt = nxt
-            pts.append(line[-1])
-            return pts
-        else:
-            return line
-
     def _getTransitionLine(self, pdc, p1, p2, obj):
         """Use an OCL PathDropCutter to generate a safe transition path between
         two points in the x/y plane."""
         p1xy, p2xy = ((p1.x, p1.y), (p2.x, p2.y))
         pdcLine = self._planarDropCutScan(pdc, p1xy, p2xy)
         if obj.OptimizeLinearPaths:
-            pdcLine = self._optimizeLinearSegments(pdcLine)
+            pdcLine = PathUtils.simplify3dLine(
+                pdcLine, tolerance=obj.LinearDeflection.Value)
         zs = [obj.z for obj in pdcLine]
         # PDC z values are based on the model, and do not take into account
         # any remaining stock / multi layer paths. Adjust raw PDC z values to
diff --git a/src/Mod/Path/PathScripts/PathUtils.py b/src/Mod/Path/PathScripts/PathUtils.py
index f5c8b0416f29..35ed2808dbd8 100644
--- a/src/Mod/Path/PathScripts/PathUtils.py
+++ b/src/Mod/Path/PathScripts/PathUtils.py
@@ -871,3 +871,41 @@ def __fixed_steps(self, start, stop, size):
             return depths
         else:
             return [stop] + depths
+
+
+def simplify3dLine(line, tolerance=1e-4):
+    """Simplify a line defined by a list of App.Vectors, while keeping the
+    maximum deviation from the original line within the defined tolerance.
+    Implementation of
+    https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm"""
+    stack = [(0, len(line) - 1)]
+    results = []
+
+    def processRange(start, end):
+        """Internal worker. Process a range of Vector indices within the
+        line."""
+        if end - start < 2:
+            results.extend(line[start:end])
+            return
+        # Find point with maximum distance
+        maxIndex, maxDistance = 0, 0.0
+        startPoint, endPoint = (line[start], line[end])
+        for i in range(start + 1, end):
+            v = line[i]
+            distance = v.distanceToLineSegment(startPoint, endPoint).Length
+            if distance > maxDistance:
+                maxDistance = distance
+                maxIndex = i
+        if maxDistance > tolerance:
+            # Push second branch first, to be executed last
+            stack.append((maxIndex, end))
+            stack.append((start, maxIndex))
+        else:
+            results.append(line[start])
+
+    while len(stack):
+        processRange(*stack.pop())
+    # Each segment only appended its start point to the final result, so fill in
+    # the last point.
+    results.append(line[-1])
+    return results