1 files changed, 63 insertions, 0 deletions

diff --git a/server/src/geo/multizonal.cppm b/server/src/geo/multizonal.cppm
index ed78a24..7930367 100644
--- a/server/src/geo/multizonal.cppm
+++ b/server/src/geo/multizonal.cppm
@@ -241,4 +241,67 @@ public:
   }
 };
 
+template<class T>
+class point_rtree
+{
+public:
+  using index_value = std::tuple<utm::zone_local::prim::box, utm::zone_local::prim::point, T>;
+
+private:
+  multi_zone<bgeo::index::rtree<index_value, bgeo::index::quadratic<16>>> local_rtrees_;
+
+public:
+  template<std::convertible_to<T> U>
+  auto insert(utm::zonable_wgs84_point const& p, U&& arg) -> void
+  {
+    auto to_utm = wgs84::utm_transforms::instance();
+
+    auto p_zone = utm::zone::for_wgs84_point(p);
+    auto mp_alt_zone = std::optional<utm::zone>{};
+    if (auto [neighbor_zone, neighbor_zone_dist] = wgs84::neighbor_utm_zone(p);
+        neighbor_zone_dist < 30.0 /* m */)
+      mp_alt_zone = neighbor_zone;
+
+    auto insert = [&](utm::zone z) -> void
+    {
+      auto p_utm = to_utm[z].apply(p);
+      auto bpe = bgeo::return_buffer<utm::zone_local::prim::box>(bgeo::return_envelope<utm::zone_local::prim::box>(p_utm), 30.0 /* m */);
+      local_rtrees_[z].insert(index_value{bpe, p_utm, arg});
+    };
+
+    insert(p_zone);
+    if (mp_alt_zone)
+      insert(*mp_alt_zone);
+  }
+
+  auto size() const -> std::size_t
+  {
+    auto total_size = 0uz;
+    for (auto z = utm::zone::min(); z != utm::zone::max(); z = z.next())
+      total_size += local_rtrees_[z].size();
+    return total_size;
+  }
+
+  // Note: the same T may be generated more than once!
+  auto intersection(zoned_linestring_seg_seq const& lss) const -> std::generator<T const&>
+  {
+    for (auto const& ls : lss.segments)
+    {
+      // auto ls_box = bgeo::return_envelope<utm::zone_local::prim::box>(ls);
+      for (auto it = local_rtrees_[ls.zone].qbegin(bgeo::index::intersects(ls));
+           it != local_rtrees_[ls.zone].qend(); it++)
+      {
+        using multi_linestring = bgeo::model::multi_linestring<utm::zone_local::prim::linestring>;
+
+        auto relevant_ls_parts = multi_linestring{};
+        bgeo::intersection(std::get<0>(*it), static_cast<utm::zone_local::prim::linestring const&>(ls), relevant_ls_parts);
+        if (bgeo::distance(relevant_ls_parts, std::get<1>(*it)) < 30.0 /* m */)
+        {
+          co_yield std::get<2>(*it);
+        }
+      }
+    }
+  }
+};
+
 } // namespace routemon::geo::multizonal