fix polylines: iterative rdp, doc XY-only simplify, py39 compat, validate threshold

Author: jf---

src/compas_cgal/polylines.py

@@ -1,7 +1,5 @@
 """Polyline utilities using CGAL."""
 
-from __future__ import annotations
-
 from typing import TYPE_CHECKING
 from typing import List
 from typing import Union
@@ -20,9 +18,12 @@
 __all__ = ["simplify_polylines", "simplify_polyline", "closest_points_on_polyline"]
 
 
-def simplify_polylines(polylines: list[PointsList], threshold: float) -> list[NDArray]:
+def simplify_polylines(polylines: List[PointsList], threshold: float) -> List[NDArray]:
     """Simplify multiple polylines using Douglas-Peucker algorithm.
 
+    Simplification is performed in the XY plane only. For 3D polylines,
+    Z coordinates are preserved but not considered in distance calculations.
+
     Parameters
     ----------
     polylines : list of array-like
@@ -45,13 +46,18 @@ def simplify_polylines(polylines: list[PointsList], threshold: float) -> list[ND
     3
 
     """
+    if threshold < 0:
+        raise ValueError("threshold must be non-negative")
     arrays = [np.asarray(p, dtype=np.float64) for p in polylines]
     return _simplify(arrays, threshold)
 
 
 def simplify_polyline(polyline: PointsList, threshold: float) -> NDArray:
     """Simplify a single polyline using Douglas-Peucker algorithm.
 
+    Simplification is performed in the XY plane only. For 3D polylines,
+    Z coordinates are preserved but not considered in distance calculations.
+
     Parameters
     ----------
     polyline : array-like

src/polylines.cpp

@@ -33,37 +33,46 @@ static double perpendicular_distance(
     return std::sqrt((px - proj_x) * (px - proj_x) + (py - proj_y) * (py - proj_y));
 }
 
-// Recursive Douglas-Peucker
+// Iterative Douglas-Peucker (avoids stack overflow on large polylines)
 static void rdp_simplify(
     const compas::RowMatrixXd& polyline,
     int start, int end,
     double threshold,
     std::vector<bool>& keep)
 {
-    if (end <= start + 1) return;
-
-    double max_dist = 0.0;
-    int max_idx = start;
-
-    double x1 = polyline(start, 0);
-    double y1 = polyline(start, 1);
-    double x2 = polyline(end, 0);
-    double y2 = polyline(end, 1);
-
-    for (int i = start + 1; i < end; ++i) {
-        double dist = perpendicular_distance(
-            polyline(i, 0), polyline(i, 1),
-            x1, y1, x2, y2);
-        if (dist > max_dist) {
-            max_dist = dist;
-            max_idx = i;
+    std::vector<std::pair<int, int>> stack;
+    stack.reserve(128);  // Reasonable default to avoid reallocs
+    stack.push_back({start, end});
+
+    while (!stack.empty()) {
+        auto [seg_start, seg_end] = stack.back();
+        stack.pop_back();
+
+        if (seg_end <= seg_start + 1) continue;
+
+        double max_dist = 0.0;
+        int max_idx = seg_start;
+
+        double x1 = polyline(seg_start, 0);
+        double y1 = polyline(seg_start, 1);
+        double x2 = polyline(seg_end, 0);
+        double y2 = polyline(seg_end, 1);
+
+        for (int i = seg_start + 1; i < seg_end; ++i) {
+            double dist = perpendicular_distance(
+                polyline(i, 0), polyline(i, 1),
+                x1, y1, x2, y2);
+            if (dist > max_dist) {
+                max_dist = dist;
+                max_idx = i;
+            }
         }
-    }
 
-    if (max_dist > threshold) {
-        keep[max_idx] = true;
-        rdp_simplify(polyline, start, max_idx, threshold, keep);
-        rdp_simplify(polyline, max_idx, end, threshold, keep);
+        if (max_dist > threshold) {
+            keep[max_idx] = true;
+            stack.push_back({seg_start, max_idx});
+            stack.push_back({max_idx, seg_end});
+        }
     }
 }
 

src/polylines.h

@@ -8,10 +8,13 @@
 #include <CGAL/AABB_segment_primitive.h>
 
 /**
- * @brief Simplify multiple 2D polylines using Douglas-Peucker algorithm.
+ * @brief Simplify polylines using Douglas-Peucker algorithm.
+ *
+ * Simplification is performed in the XY plane only. For 3D polylines,
+ * Z coordinates are preserved but not considered in distance calculations.
  *
  * @param polylines Vector of polylines, each as Nx2 or Nx3 matrix (float64)
- * @param threshold Squared distance threshold for simplification
+ * @param threshold Distance threshold for simplification
  * @return std::vector<compas::RowMatrixXd> Simplified polylines
  */
 std::vector<compas::RowMatrixXd>