Misc updates (#16)

district10 · web-flow · commit 6abc08892aa7 · 2023-03-11T12:25:45.000+08:00
* add iter version of douglas simplify, sync from https://github.com/cubao/pybind11-rdp * add range to segment index, t * fix * test for dir * fix polyline dir, update test * update version * update header * update header * update pybind11 * fix CI * fix CI * fix * fix * update headers * test douglas, update docs
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -32,7 +32,6 @@ repos:
   - id: mixed-line-ending
   - id: requirements-txt-fixer
   - id: trailing-whitespace
-  - id: end-of-file-fixer
 
 # Black, the code formatter, natively supports pre-commit
 - repo: https://github.com/psf/black
diff --git a/docs/about/release-notes.md b/docs/about/release-notes.md
@@ -10,6 +10,12 @@ To upgrade `polyline-ruler` to the latest version, use pip:
 pip install -U polyline-ruler
 ```
 
+## Version 0.0.4 (2023-03-11)
+
+*   Misc updates
+    *   Update dir logic
+    *   Integrate douglas_simplify iter
+
 ## Version 0.0.3 (2023-03-05)
 
 *   Use GitHub workflow to release to pypi
diff --git a/headers b/headers
@@ -1 +1 @@
-Subproject commit 512ca351c24b03b85c0034d418f1959ef5c4a219
+Subproject commit 2810720d3726d170cd15abcb782012e893f42470
diff --git a/mkdocs.yml b/mkdocs.yml
@@ -15,4 +15,4 @@ nav:
         - License: about/license.md
 plugins:
   - include-markdown
-copyright: Copyright &copy; 2023 <a href="https://github.com/cubao">cubao team</a>.
+copyright: Copyright &copy; 2023 <a href="https://cubao.readthedocs.io">cubao team</a>.
diff --git a/pybind11 b/pybind11
@@ -1 +1 @@
-Subproject commit 914c06fb252b6cc3727d0eedab6736e88a3fcb01
+Subproject commit 442261da585536521ff459b1457b2904895f23b4
diff --git a/setup.py b/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="polyline_ruler",
-    version="0.0.3",
+    version="0.0.4",
     author="tzx",
     author_email="dvorak4tzx@gmail.com",
     url="https://polyline-ruler.readthedocs.io",
diff --git a/src/polyline_ruler.hpp b/src/polyline_ruler.hpp
@@ -13,6 +13,7 @@
 
 #include <Eigen/Core>
 #include <optional>
+#include <queue>
 
 #include "cheap_ruler.hpp"
 #include "crs_transform.hpp"
@@ -291,6 +292,22 @@ struct PolylineRuler
         return ranges[seg_idx] * (1.0 - t) + ranges[seg_idx + 1] * t;
     }
 
+    int segment_index(double range) const
+    {
+        const double *ranges = this->ranges().data();
+        int I = std::upper_bound(ranges, ranges + N_, range) - ranges;
+        return std::min(std::max(0, I - 1), N_ - 2);
+    }
+
+    std::pair<int, double> segment_index_t(double range) const
+    {
+        const double *ranges = this->ranges().data();
+        int I = std::upper_bound(ranges, ranges + N_, range) - ranges;
+        int i = std::min(std::max(0, I - 1), N_ - 2);
+        double t = (range - ranges[i]) / (ranges[i + 1] - ranges[i]);
+        return {i, t};
+    }
+
     double length() const { return ranges()[N_ - 1]; }
 
     static RowVectors dirs(const Eigen::Ref<const RowVectors> &polyline,
@@ -352,42 +369,31 @@ struct PolylineRuler
         return *dirs_;
     }
 
+    Eigen::Vector3d dir(int pt_index) const
+    {
+        return dirs().row(std::min(pt_index, N_ - 2));
+    }
+
     Eigen::Vector3d dir(double range, bool smooth_joint = true) const
     {
+        if (!smooth_joint) {
+            return dir(segment_index(range));
+        }
+        auto [i, t] = segment_index_t(range);
         auto &dirs = this->dirs();
-        if (range <= 0.0) {
+        if (i == 0) {
             return dirs.row(0);
-        }
-        auto &ranges = this->ranges();
-        int i = 0;
-        while (i + 1 < N_ && ranges[i + 1] < range) {
-            ++i;
-        }
-        if (smooth_joint && i + 1 < N_ && ranges[i + 1] == range) {
-            Eigen::Vector3d dir = dirs.row(i + 1) + dirs.row(i);
+        } else if (t == 0) {
+            Eigen::Vector3d dir = dirs.row(i - 1) + dirs.row(i);
             return dir / dir.norm();
+        } else {
+            return dirs.row(i);
         }
-        return dirs.row(std::min(i, (int)dirs.rows() - 1));
     }
 
     Eigen::Vector3d extended_along(double range) const
     {
-        auto &ranges = this->ranges();
-        if (range <= 0.0) {
-            double t = range / ranges[1];
-            return interpolate(polyline_.row(0), polyline_.row(1), t,
-                               is_wgs84_);
-        } else if (range >= length()) {
-            double t =
-                (range - ranges[N_ - 2]) / (ranges[N_ - 1] - ranges[N_ - 2]);
-            return interpolate(polyline_.row(N_ - 2), polyline_.row(N_ - 1), t,
-                               is_wgs84_);
-        }
-        int i = 0;
-        while (i + 1 < N_ && ranges[i + 1] < range) {
-            ++i;
-        }
-        double t = (range - ranges[i]) / (ranges[i + 1] - ranges[i]);
+        auto [i, t] = segment_index_t(range);
         return interpolate(polyline_.row(i), polyline_.row(i + 1), t,
                            is_wgs84_);
     }
@@ -742,31 +748,103 @@ inline void douglas_simplify(const Eigen::Ref<const RowVectors> &coords,
     douglas_simplify(coords, to_keep, max_index, j, epsilon);
 }
 
-inline Eigen::VectorXi
-douglas_simplify_mask(const Eigen::Ref<const RowVectors> &coords,
-                      double epsilon, bool is_wgs84 = false)
+void douglas_simplify_iter(const Eigen::Ref<const RowVectors> &coords,
+                           Eigen::VectorXi &to_keep, const double epsilon)
+{
+    std::queue<std::pair<int, int>> q;
+    q.push({0, to_keep.size() - 1});
+    while (!q.empty()) {
+        int i = q.front().first;
+        int j = q.front().second;
+        q.pop();
+        to_keep[i] = to_keep[j] = 1;
+        if (j - i <= 1) {
+            continue;
+        }
+        LineSegment line(coords.row(i), coords.row(j));
+        double max_dist2 = 0.0;
+        int max_index = 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;
+            }
+        }
+        if (max_dist2 <= epsilon * epsilon) {
+            continue;
+        }
+        q.push({i, max_index});
+        q.push({max_index, j});
+    }
+}
+
+inline Eigen::VectorXi douglas_simplify_mask(const RowVectors &coords,
+                                             double epsilon,        //
+                                             bool is_wgs84 = false, //
+                                             bool recursive = true)
 {
     if (is_wgs84) {
-        return douglas_simplify_mask(lla2enu(coords), epsilon, !is_wgs84);
+        return douglas_simplify_mask(lla2enu(coords), epsilon, //
+                                     false, recursive);
     }
     Eigen::VectorXi mask(coords.rows());
     mask.setZero();
-    douglas_simplify(coords, mask, 0, mask.size() - 1, epsilon);
+    if (recursive) {
+        douglas_simplify(coords, mask, 0, mask.size() - 1, epsilon);
+    } else {
+        douglas_simplify_iter(coords, mask, epsilon);
+    }
     return mask;
 }
 
-inline Eigen::VectorXi
-douglas_simplify_indexes(const Eigen::Ref<const RowVectors> &coords,
-                         double epsilon, bool is_wgs84 = false)
+inline Eigen::VectorXi douglas_simplify_indexes(const RowVectors &coords,
+                                                double epsilon,        //
+                                                bool is_wgs84 = false, //
+                                                bool recursive = true)
 {
-    return mask2indexes(douglas_simplify_mask(coords, epsilon, is_wgs84));
+    return mask2indexes(
+        douglas_simplify_mask(coords, epsilon, is_wgs84, recursive));
 }
 
-inline RowVectors douglas_simplify(const Eigen::Ref<const RowVectors> &coords,
-                                   double epsilon, bool is_wgs84 = false)
+inline RowVectors douglas_simplify(const RowVectors &coords,
+                                   double epsilon,        //
+                                   bool is_wgs84 = false, //
+                                   bool recursive = true)
+{
+    return select_by_mask(
+        coords, //
+        douglas_simplify_mask(coords, epsilon, is_wgs84, recursive));
+}
+
+// Nx2
+inline Eigen::VectorXi
+douglas_simplify_mask(const Eigen::Ref<const RowVectorsNx2> &coords,
+                      double epsilon,        //
+                      bool is_wgs84 = false, //
+                      bool recursive = true)
+{
+    return douglas_simplify_mask(to_Nx3(coords), epsilon, is_wgs84, recursive);
+}
+inline Eigen::VectorXi
+douglas_simplify_indexes(const Eigen::Ref<const RowVectorsNx2> &coords,
+                         double epsilon,        //
+                         bool is_wgs84 = false, //
+                         bool recursive = true)
+{
+    return douglas_simplify_indexes(to_Nx3(coords), epsilon, is_wgs84,
+                                    recursive);
+}
+inline RowVectorsNx2
+douglas_simplify(const Eigen::Ref<const RowVectorsNx2> &coords,
+                 double epsilon,        //
+                 bool is_wgs84 = false, //
+                 bool recursive = true)
 {
-    return select_by_mask(coords,
-                          douglas_simplify_mask(coords, epsilon, is_wgs84));
+    RowVectorsNx2 ret =
+        douglas_simplify(to_Nx3(coords), epsilon, is_wgs84, recursive)
+            .leftCols(2);
+    return ret;
 }
 
 } // namespace cubao
diff --git a/src/pybind11_polyline_ruler.hpp b/src/pybind11_polyline_ruler.hpp
@@ -93,6 +93,9 @@ CUBAO_INLINE void bind_polyline_ruler(py::module &m)
              py::overload_cast<int, double>(&PolylineRuler::range, py::const_),
              py::kw_only(), "segment_index"_a, "t"_a)
         //
+        .def("segment_index", &PolylineRuler::segment_index, "range"_a)
+        .def("segment_index_t", &PolylineRuler::segment_index_t, "range"_a)
+        //
         .def("length", &PolylineRuler::length)
         //
         .def_static(
@@ -103,17 +106,19 @@ CUBAO_INLINE void bind_polyline_ruler(py::module &m)
         .def("dirs", py::overload_cast<>(&PolylineRuler::dirs, py::const_),
              rvp::reference_internal)
         //
+        .def("dir", py::overload_cast<int>(&PolylineRuler::dir, py::const_),
+             py::kw_only(), "point_index"_a)
         .def("dir",
              py::overload_cast<double, bool>(&PolylineRuler::dir, py::const_),
-             "range"_a, py::kw_only(), "smooth_joint"_a = true)
+             py::kw_only(), "range"_a, "smooth_joint"_a = true)
         .def("extended_along", &PolylineRuler::extended_along, "range"_a)
         .def("at", py::overload_cast<double>(&PolylineRuler::at, py::const_),
-             "range"_a)
+             py::kw_only(), "range"_a)
         .def("at", py::overload_cast<int>(&PolylineRuler::at, py::const_),
-             "segment_index"_a)
+             py::kw_only(), "segment_index"_a)
         .def("at",
              py::overload_cast<int, double>(&PolylineRuler::at, py::const_),
-             "segment_index"_a, py::kw_only(), "t"_a)
+             py::kw_only(), "segment_index"_a, "t"_a)
         .def("arrow", &PolylineRuler::arrow, "range"_a, //
              py::kw_only(), "smooth_joint"_a = true)
         .def("arrows",
@@ -201,5 +206,47 @@ CUBAO_INLINE void bind_polyline_ruler(py::module &m)
                     "A"_a, "B"_a, py::kw_only(), "t"_a, "is_wgs84"_a = false)
         //
         ;
+
+    m.def("douglas_simplify",
+          py::overload_cast<const RowVectors &, double, bool,
+                            bool>(&douglas_simplify), //
+          "coords"_a, "epsilon"_a,                    //
+          py::kw_only(),                              //
+          "is_wgs84"_a = false,                       //
+          "recursive"_a = true);
+    m.def(
+        "douglas_simplify",
+        py::overload_cast<const Eigen::Ref<const RowVectorsNx2> &, double, bool,
+                          bool>(&douglas_simplify), //
+        "coords"_a, "epsilon"_a,                    //
+        py::kw_only(),                              //
+        "is_wgs84"_a = false,                       //
+        "recursive"_a = true);
+    m.def("douglas_simplify_mask",
+          py::overload_cast<const RowVectors &, double, bool,
+                            bool>(&douglas_simplify_mask), //
+          "coords"_a, "epsilon"_a, py::kw_only(),          //
+          "is_wgs84"_a = false,                            //
+          "recursive"_a = true);
+    m.def(
+        "douglas_simplify_mask",
+        py::overload_cast<const Eigen::Ref<const RowVectorsNx2> &, double, bool,
+                          bool>(&douglas_simplify_mask), //
+        "coords"_a, "epsilon"_a, py::kw_only(),          //
+        "is_wgs84"_a = false,                            //
+        "recursive"_a = true);
+    m.def("douglas_simplify_indexes",
+          py::overload_cast<const RowVectors &, double, bool,
+                            bool>(&douglas_simplify_indexes), //
+          "coords"_a, "epsilon"_a, py::kw_only(),             //
+          "is_wgs84"_a = false,                               //
+          "recursive"_a = true);
+    m.def(
+        "douglas_simplify_indexes",
+        py::overload_cast<const Eigen::Ref<const RowVectorsNx2> &, double, bool,
+                          bool>(&douglas_simplify_indexes), //
+        "coords"_a, "epsilon"_a, py::kw_only(),             //
+        "is_wgs84"_a = false,                               //
+        "recursive"_a = true);
 }
 } // namespace cubao
diff --git a/tests/test_basic.py b/tests/test_basic.py
