handle degenerate case (#14)

* not ready

* just merge into master

---------

Co-authored-by: TANG ZHIXIONG <zhixiong.tang@momenta.ai>

Author: district10

Makefile

@@ -35,7 +35,7 @@ test_in_linux:
 
 PYTHON ?= python3
 python_install:
-	$(PYTHON) setup.py install
+	$(PYTHON) setup.py install --force
 python_build:
 	$(PYTHON) setup.py bdist_wheel
 python_sdist:

docs/about/release-notes.md

@@ -10,6 +10,15 @@ To upgrade `pybind11-rdp` to the latest version, use pip:
 pip install -U pybind11-rdp
 ```
 
+## Version 0.1.4 (2023-07-28)
+
+*   Handle degenerate case, related to <https://github.com/mapbox/geojson-vt/issues/104>
+*   How to test? Use `ulimit -s 100 && python3 test.py`?
+
+## Version 0.1.3 (2023-07-28)
+
+*   Update docs, update packaging
+
 ## Version 0.1.2 (2023-03-02)
 
 *   Identical API to rdp, notice difference

setup.py

@@ -122,7 +122,7 @@ def build_extension(self, ext):
 # logic and declaration, and simpler if you include description/version in a file.
 setup(
     name="pybind11_rdp",
-    version="0.1.3",
+    version="0.1.4",
     author="tzx",
     author_email="dvorak4tzx@gmail.com",
     url="https://github.com/cubao/pybind11-rdp",

src/main.cpp

@@ -58,11 +58,22 @@ void douglas_simplify(const Eigen::Ref<const RowVectors> &coords,
     LineSegment line(coords.row(i), coords.row(j));
     double max_dist2 = 0.0;
     int max_index = i;
+    int mid = i + (j - i) / 2;
+    int min_pos_to_mid = j - i;
     for (int k = i + 1; k < j; ++k) {
         double dist2 = line.distance2(coords.row(k));
         if (dist2 > max_dist2) {
             max_dist2 = dist2;
             max_index = k;
+        } else if (dist2 == max_dist2) {
+            // a workaround to ensure we choose a pivot close to the middle of
+            // the list, reducing recursion depth, for certain degenerate inputs
+            // https://github.com/mapbox/geojson-vt/issues/104
+            int pos_to_mid = std::fabs(k - mid);
+            if (pos_to_mid < min_pos_to_mid) {
+                min_pos_to_mid = pos_to_mid;
+                max_index = k;
+            }
         }
     }
     if (max_dist2 <= epsilon * epsilon) {
@@ -88,11 +99,22 @@ void douglas_simplify_iter(const Eigen::Ref<const RowVectors> &coords,
         LineSegment line(coords.row(i), coords.row(j));
         double max_dist2 = 0.0;
         int max_index = i;
+        int mid = i + (j - i) / 2;
+        int min_pos_to_mid = j - i;
         for (int k = i + 1; k < j; ++k) {
             double dist2 = line.distance2(coords.row(k));
             if (dist2 > max_dist2) {
                 max_dist2 = dist2;
                 max_index = k;
+            } else if (dist2 == max_dist2) {
+                // a workaround to ensure we choose a pivot close to the middle
+                // of the list, reducing recursion depth, for certain degenerate
+                // inputs https://github.com/mapbox/geojson-vt/issues/104
+                int pos_to_mid = std::fabs(k - mid);
+                if (pos_to_mid < min_pos_to_mid) {
+                    min_pos_to_mid = pos_to_mid;
+                    max_index = k;
+                }
             }
         }
         if (max_dist2 <= epsilon * epsilon) {

tests/test_basic.py

@@ -1,6 +1,8 @@
 import os
 import sys
+import time
 
+import numpy as np
 import pytest
 
 from pybind11_rdp import LineSegment, rdp
@@ -23,6 +25,26 @@ def test_rdp():
     assert rdp([[0, 0], [5, 1 - 1e-3], [10, 0]], epsilon=1).shape == (2, 2)
 
 
+def test_degenerate_case():
+    # https://github.com/mapbox/geojson-vt/issues/104
+    coords = []
+    for _ in range(14000):
+        coords.extend(
+            [
+                [0.0, 0.0],
+                [1.0, 0.0],
+                [1.0, 1.0],
+                [0.0, 1.0],
+            ]
+        )
+    coords = np.array(coords)
+    tick = time.time()
+    ret = rdp(coords, 2e-15, algo="recursive")
+    tock = time.time()
+    print(tock - tick, "secs")  # 4 sec
+    assert len(ret) == len(coords)
+
+
 def pytest_main(dir: str, *, test_file: str = None):
     os.chdir(dir)
     sys.exit(