Add Internal Function intersection_count_upper_bound

cpp/include/cuspatial/cuda_utils.hpp

@@ -15,9 +15,23 @@
  */
 
 #pragma once
+#include <utility>
 
 #ifdef __CUDACC__
 #define CUSPATIAL_HOST_DEVICE __host__ __device__
 #else
 #define CUSPATIAL_HOST_DEVICE
 #endif
+
+/**
+ * @brief Return kernel launch parameters for 1D grid with total `n` threads.
+ *
+ * @tparam threads_per_block Number of threads per block
+ * @param n Number of threads
+ * @return Threads per block and number of blocks
+ */
+template <std::size_t threads_per_block = 256>
+std::pair<std::size_t, std::size_t> constexpr grid_1d(std::size_t const n)
+{
+  return {threads_per_block, (n + threads_per_block - 1) / threads_per_block};
+}

cpp/include/cuspatial/detail/utility/device_atomics.cuh

@@ -16,7 +16,11 @@
 
 #pragma once
 
+#include <thrust/detail/raw_reference_cast.h>
 #include <thrust/device_ptr.h>
+#include <thrust/device_reference.h>
+
+#include <cuda/atomic>
 
 #include <type_traits>
 
@@ -260,5 +264,24 @@ __device__ inline float atomicMax(thrust::device_ptr<float> ptr, float val)
   return atomicMax(thrust::raw_pointer_cast(ptr), val);
 }
 
+/**
+ * @brief Factory method to create atomic_ref from a thrust::device_reference
+ */
+template <typename T, cuda::thread_scope Scope = cuda::thread_scope_device>
+auto __device__ make_atomic_ref(thrust::device_reference<T> ref)
+{
+  T& raw_ref = thrust::raw_reference_cast(ref);
+  return cuda::atomic_ref<T, Scope>{raw_ref};
+}
+
+/**
+ * @brief Factory method to create atomic_ref from raw reference
+ */
+template <typename T, cuda::thread_scope Scope = cuda::thread_scope_device>
+auto __device__ make_atomic_ref(T& ref)
+{
+  return cuda::atomic_ref<T, Scope>{ref};
+}
+
 }  // namespace detail
 }  // namespace cuspatial

cpp/include/cuspatial/experimental/detail/linestring_intersection_count.cuh

@@ -0,0 +1,93 @@
+/*
+ * Copyright (c) 2022, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/detail/utility/device_atomics.cuh>
+#include <cuspatial/detail/utility/linestring.cuh>
+
+#include <rmm/cuda_stream_view.hpp>
+
+#include <thrust/tuple.h>
+
+namespace cuspatial {
+namespace detail {
+
+template <typename MultiLinestringRange1,
+          typename MultiLinestringRange2,
+          typename OutputIt1,
+          typename OutputIt2>
+__global__ void count_intersection_and_overlaps_simple(MultiLinestringRange1 multilinestrings1,
+                                                       MultiLinestringRange2 multilinestrings2,
+                                                       OutputIt1 point_count_it,
+                                                       OutputIt2 segment_count_it)
+{
+  using T = typename MultiLinestringRange1::element_t;
+  for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multilinestrings1.num_points();
+       idx += gridDim.x * blockDim.x) {
+    auto const part_idx = multilinestrings1.part_idx_from_point_idx(idx);
+    if (!multilinestrings1.is_valid_segment_id(idx, part_idx)) continue;
+    auto const geometry_idx = multilinestrings1.geometry_idx_from_part_idx(part_idx);
+    auto [a, b]             = multilinestrings1.segment(idx);
+
+    for (auto const& linestring2 : multilinestrings2[geometry_idx]) {
+      for (auto [c, d] : linestring2) {
+        auto [point_opt, segment_opt] = segment_intersection(segment<T>{a, b}, segment<T>{c, d});
+        if (point_opt.has_value()) {
+          auto r = make_atomic_ref(*(point_count_it + geometry_idx));
+          r.fetch_add(1);
+        } else if (segment_opt.has_value()) {
+          auto r = make_atomic_ref(*(segment_count_it + geometry_idx));
+          r.fetch_add(1);
+        }
+      }
+    }
+  }
+}
+
+/**
+ * @internal
+ * @brief Count the upper bound of intersecting linestrings between a pair of multilinestring range
+ *
+ * @tparam MultiLinestringRange1 Type of first multilinestring range
+ * @tparam MultiLinestringRange2 Type of second multilinestring range
+ * @tparam OutputIt1 Type of intersecting point count iterator
+ * @tparam OutputIt2 Type of overlapping segment count iterator
+ * @param multilinestrings1 The first multilinestring range
+ * @param multilinestrings2 The second multilinestring range
+ * @param points_count_it Integral iterator to the number of intersecting points
+ * @param segments_count_it Integral iterator to the number of overlapping segments
+ * @param stream The CUDA stream for device memory operations
+ */
+template <typename MultiLinestringRange1,
+          typename MultiLinestringRange2,
+          typename OutputIt1,
+          typename OutputIt2>
+void pairwise_linestring_intersection_upper_bound_count(MultiLinestringRange1 multilinestrings1,
+                                                        MultiLinestringRange2 multilinestrings2,
+                                                        OutputIt1 points_count_it,
+                                                        OutputIt2 segments_count_it,
+                                                        rmm::cuda_stream_view stream)
+{
+  auto [threads_per_block, num_blocks] = grid_1d(multilinestrings1.size());
+  detail::
+    count_intersection_and_overlaps_simple<<<num_blocks, threads_per_block, 0, stream.value()>>>(
+      multilinestrings1, multilinestrings2, points_count_it, segments_count_it);
+}
+
+}  // namespace detail
+}  // namespace cuspatial

cpp/include/cuspatial_test/vector_factories.cuh

@@ -14,12 +14,64 @@
  * limitations under the License.
  */
 
+#include <cuspatial/experimental/ranges/multilinestring_range.cuh>
+#include <cuspatial/vec_2d.hpp>
+
 #include <rmm/device_vector.hpp>
 
 #include <initializer_list>
+#include <vector>
+
+namespace cuspatial {
+namespace test {
 
 template <typename T>
 auto make_device_vector(std::initializer_list<T> inl)
 {
   return rmm::device_vector<T>(inl.begin(), inl.end());
 }
+
+template <typename GeometryArray, typename PartArray, typename CoordinateArray>
+class multilinestring_array {
+ public:
+  multilinestring_array(GeometryArray geometry_offsets_array,
+                        PartArray part_offsets_array,
+                        CoordinateArray coordinate_offset_array)
+    : _geometry_offset_array(geometry_offsets_array),
+      _part_offset_array(part_offsets_array),
+      _coordinate_offset_array(coordinate_offset_array)
+  {
+  }
+
+  auto size() { return _geometry_offset_array.size() - 1; }
+
+  auto range()
+  {
+    return multilinestring_range(_geometry_offset_array.begin(),
+                                 _geometry_offset_array.end(),
+                                 _part_offset_array.begin(),
+                                 _part_offset_array.end(),
+                                 _coordinate_offset_array.begin(),
+                                 _coordinate_offset_array.end());
+  }
+
+ protected:
+  GeometryArray _geometry_offset_array;
+  PartArray _part_offset_array;
+  CoordinateArray _coordinate_offset_array;
+};
+
+template <typename T>
+auto make_multilinestring_array(std::initializer_list<std::size_t> geometry_inl,
+                                std::initializer_list<std::size_t> part_inl,
+                                std::initializer_list<vec_2d<T>> coord_inl)
+{
+  return multilinestring_array(
+    rmm::device_vector<std::size_t>(
+      std::vector<std::size_t>(geometry_inl.begin(), geometry_inl.end())),
+    rmm::device_vector<std::size_t>(std::vector<unsigned int>(part_inl.begin(), part_inl.end())),
+    rmm::device_vector<vec_2d<T>>(std::vector<vec_2d<T>>(coord_inl.begin(), coord_inl.end())));
+}
+
+}  // namespace test
+}  // namespace cuspatial

cpp/tests/CMakeLists.txt

@@ -126,6 +126,9 @@ ConfigureTest(LINESTRING_DISTANCE_TEST_EXP
     experimental/spatial/linestring_distance_test.cu
     experimental/spatial/linestring_distance_test_medium.cu)
 
+ConfigureTest(LINESTRING_INTERSECTION_TEST_EXP
+    experimental/spatial/linestring_intersection_count_test.cu)
+
 ConfigureTest(POINT_LINESTRING_NEAREST_POINT_TEST_EXP
     experimental/spatial/point_linestring_nearest_points_test.cu)