3 files changed, 73 insertions, 17 deletions

diff --git a/server/src/api.cpp b/server/src/api.cpp
index 137666b..e35e1b3 100644
--- a/server/src/api.cpp
+++ b/server/src/api.cpp
@@ -95,7 +95,7 @@ auto tag_invoke(json::value_from_tag, json::value& jv, sysinfo const& info)
     {"using_publication_of",
      std::format("{:%FT%TZ}", info.using_publication_of)},
     {"blse_index_size", info.blse_index_size},
-    {"p_index_size", info.p_index_size},
+    {"bpe_index_size", info.bpe_index_size},
   };
 }
 
@@ -124,7 +124,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub)
       }
       for (auto p : rc->relevant_points)
       {
-        p_index_.insert(p_index_value{p, rc});
+        bpe_index_.insert(p, rc);
       }
     }
 
@@ -141,7 +141,7 @@ handler::handler(log::logger const& l, datex2::situation_publication pub)
   auto const dur_build_s = chrono::duration<double>{after_build - before_build};
   l_.info("Indexes built in {}", dur_build_s);
   l_.info("BLSE index size: {}", blse_index_.size());
-  l_.info("Point index size: {}", p_index_.size());
+  l_.info("Point index size: {}", bpe_index_.size());
 }
 
 auto handler::process_gpx(gpx::file&& gpx_file)
@@ -194,17 +194,12 @@ auto handler::process_gpx(gpx::file&& gpx_file)
       //     relevant_road_closures.emplace(rc);
       ls_checked++;
     }
-    for (auto it = p_index_.qbegin(bgeo::index::intersects(part_box));
-         it != p_index_.qend(); it++)
+    for (auto const& rc : bpe_index_.intersection(part_zone_lss))
     {
-      // Cannot use structured bindings here for the same reason as above.
-      geo::utm::zonable_wgs84_point const& p = std::get<0>(*it);
-      std::shared_ptr<datex2::road_closure> const& rc = std::get<1>(*it);
       if (rc->validity
           && rc->validity->intersect(check_periods).periods().empty())
         continue;
-      if (bgeo::distance(p, part, vincenty_strategy) < min_distance_)
-        relevant_road_closures.emplace(rc);
+      relevant_road_closures.emplace(rc);
       p_checked++;
     }
   }
@@ -285,7 +280,7 @@ auto handler::sysinfo() -> struct sysinfo
   return {
     .using_publication_of = pub_.publication_time,
     .blse_index_size = blse_index_.size(),
-    .p_index_size = p_index_.size(),
+    .bpe_index_size = bpe_index_.size(),
   };
 }
 
diff --git a/server/src/api.cppm b/server/src/api.cppm
index 77e0f85..0bb0868 100644
--- a/server/src/api.cppm
+++ b/server/src/api.cppm
@@ -57,7 +57,7 @@ struct sysinfo
 {
   time::timestamp using_publication_of;
   std::size_t blse_index_size;
-  std::size_t p_index_size;
+  std::size_t bpe_index_size;
 };
 
 auto tag_invoke(
@@ -82,15 +82,13 @@ class handler
 
   using blse_index_value = std::shared_ptr<datex2::road_closure>;
   using blse_index = geo::multizonal::linestring_rtree<blse_index_value>;
-  using p_index_value = std::pair<
-      geo::utm::zonable_wgs84_point, std::shared_ptr<datex2::road_closure>
-  >;
-  using p_index = bgeo::index::rtree<p_index_value, bgeo::index::quadratic<16>>;
+  using bpe_index_value = std::shared_ptr<datex2::road_closure>;
+  using bpe_index = geo::multizonal::point_rtree<bpe_index_value>;
 
   log::logger l_;
   datex2::situation_publication pub_;
   blse_index blse_index_;
-  p_index p_index_;
+  bpe_index bpe_index_;
 
 public:
   explicit handler(log::logger const& l, datex2::situation_publication pub);
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