Add C++ API for point_linestring_nearest_points

cpp/CMakeLists.txt

@@ -124,6 +124,7 @@ add_library(cuspatial
     src/spatial/linestring_distance.cu
     src/spatial/point_distance.cu
     src/spatial/point_linestring_distance.cu
+    src/spatial/point_linestring_nearest_points.cu
     src/spatial/lonlat_to_cartesian.cu
     src/trajectory/derive_trajectories.cu
     src/trajectory/trajectory_bounding_boxes.cu

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

@@ -16,6 +16,8 @@
 
 #pragma once
 
+#include <thrust/tuple.h>
+
 #include <cuspatial/vec_2d.hpp>
 
 namespace cuspatial {
@@ -40,24 +42,44 @@ endpoint_index_of_linestring(SizeType const& linestring_idx,
 
 /**
  * @internal
- * @brief Computes shortest distance between @p c and segment ab
+ * @brief Computes shortest distance and nearest point between @p c and segment
+ * ab
  */
 template <typename T>
-__forceinline__ T __device__ point_to_segment_distance_squared(vec_2d<T> const& c,
-                                                               vec_2d<T> const& a,
-                                                               vec_2d<T> const& b)
+__forceinline__ thrust::tuple<T, vec_2d<T>> __device__
+point_to_segment_distance_squared_nearest_point(vec_2d<T> const& c,
+                                                vec_2d<T> const& a,
+                                                vec_2d<T> const& b)
 {
   auto ab        = b - a;
   auto ac        = c - a;
   auto l_squared = dot(ab, ab);
-  if (l_squared == 0) { return dot(ac, ac); }
+  if (l_squared == 0) { return thrust::make_tuple(dot(ac, ac), a); }
   auto r  = dot(ac, ab);
   auto bc = c - b;
   // If the projection of `c` is outside of segment `ab`, compute point-point distance.
-  if (r <= 0 or r >= l_squared) { return std::min(dot(ac, ac), dot(bc, bc)); }
+  if (r <= 0 or r >= l_squared) {
+    auto dac = dot(ac, ac);
+    auto dbc = dot(bc, bc);
+    return dac < dbc ? thrust::make_tuple(dac, a) : thrust::make_tuple(dbc, b);
+  }
   auto p  = a + (r / l_squared) * ab;
   auto pc = c - p;
-  return dot(pc, pc);
+  return thrust::make_tuple(dot(pc, pc), p);
+}
+
+/**
+ * @internal
+ * @brief Computes shortest distance between @p c and segment ab
+ */
+template <typename T>
+__forceinline__ T __device__ point_to_segment_distance_squared(vec_2d<T> const& c,
+                                                               vec_2d<T> const& a,
+                                                               vec_2d<T> const& b)
+{
+  [[maybe_unused]] auto [distance_squared, _] =
+    point_to_segment_distance_squared_nearest_point(c, a, b);
+  return distance_squared;
 }
 
 /**

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

@@ -0,0 +1,173 @@
+/*
+ * 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/detail/utility/device_atomics.cuh>
+#include <cuspatial/detail/utility/linestring.cuh>
+#include <cuspatial/error.hpp>
+#include <cuspatial/traits.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/binary_search.h>
+
+#include <iterator>
+#include <limits>
+#include <memory>
+#include <type_traits>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @internal
+ * @brief Kernel to compute the nearest point between a multipoint and multilinestring
+ *
+ * See header only API for input parameter definitions.
+ *
+ * Each thread computes the nearest point between a pair of multipoint and multilinestring.
+ * The minimum distance between the geometries are stored in `min_distance_squared` and updated
+ * when smaller is encountered. `linestring.cuh::point_to_segment_distance_squared_nearest_point`
+ * is used to compute the nearest point on the segment and its distance to the test point.
+ */
+template <class Vec2dItA,
+          class Vec2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIt>
+void __global__
+pairwise_point_linestring_nearest_points_kernel(OffsetIteratorA points_geometry_offsets_first,
+                                                OffsetIteratorA points_geometry_offsets_last,
+                                                Vec2dItA points_first,
+                                                Vec2dItA points_last,
+                                                OffsetIteratorB linestring_geometry_offsets_first,
+                                                OffsetIteratorB linestring_geometry_offsets_last,
+                                                OffsetIteratorC linestring_part_offsets_first,
+                                                OffsetIteratorC linestring_part_offsets_last,
+                                                Vec2dItB linestring_points_first,
+                                                Vec2dItB linestring_points_last,
+                                                OutputIt output_first)
+{
+  using T         = iterator_vec_base_type<Vec2dItA>;
+  using IndexType = iterator_value_type<OffsetIteratorA>;
+
+  auto num_pairs =
+    thrust::distance(points_geometry_offsets_first, points_geometry_offsets_last) - 1;
+  auto num_linestring_points = thrust::distance(linestring_points_first, linestring_points_last);
+
+  for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_pairs;
+       idx += gridDim.x * blockDim.x) {
+    IndexType nearest_point_idx;
+    IndexType nearest_part_idx;
+    IndexType nearest_segment_idx;
+    vec_2d<T> nearest_point;
+
+    T min_distance_squared = std::numeric_limits<T>::max();
+    IndexType point_start  = points_geometry_offsets_first[idx];
+    IndexType point_end    = points_geometry_offsets_first[idx + 1];
+    for (auto point_idx = point_start; point_idx < point_end; point_idx++) {
+      IndexType linestring_parts_start = linestring_geometry_offsets_first[idx];
+      IndexType linestring_parts_end   = linestring_geometry_offsets_first[idx + 1];
+
+      for (auto part_idx = linestring_parts_start; part_idx < linestring_parts_end; part_idx++) {
+        IndexType segment_start = linestring_part_offsets_first[part_idx];
+        // The last point of the linestring does not form a segment
+        IndexType segment_end = linestring_part_offsets_first[part_idx + 1] - 1;
+
+        for (auto segment_idx = segment_start; segment_idx < segment_end; segment_idx++) {
+          vec_2d<T> c = points_first[point_idx];
+          vec_2d<T> a = linestring_points_first[segment_idx];
+          vec_2d<T> b = linestring_points_first[segment_idx + 1];
+
+          auto distance_nearest_point_pair =
+            point_to_segment_distance_squared_nearest_point(c, a, b);
+          auto distance_squared = thrust::get<0>(distance_nearest_point_pair);
+          if (distance_squared < min_distance_squared) {
+            min_distance_squared = distance_squared;
+            nearest_point_idx    = point_idx - point_start;
+            nearest_part_idx     = part_idx - linestring_parts_start;
+            nearest_segment_idx  = segment_idx - segment_start;
+            nearest_point        = thrust::get<1>(distance_nearest_point_pair);
+          }
+        }
+      }
+    }
+    output_first[idx] =
+      thrust::make_tuple(nearest_point_idx, nearest_part_idx, nearest_segment_idx, nearest_point);
+  }
+}
+
+}  // namespace detail
+
+template <class Vec2dItA,
+          class Vec2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIt>
+OutputIt pairwise_point_linestring_nearest_points(OffsetIteratorA points_geometry_offsets_first,
+                                                  OffsetIteratorA points_geometry_offsets_last,
+                                                  Vec2dItA points_first,
+                                                  Vec2dItA points_last,
+                                                  OffsetIteratorB linestring_geometry_offsets_first,
+                                                  OffsetIteratorC linestring_part_offsets_first,
+                                                  OffsetIteratorC linestring_part_offsets_last,
+                                                  Vec2dItB linestring_points_first,
+                                                  Vec2dItB linestring_points_last,
+                                                  OutputIt output_first,
+                                                  rmm::cuda_stream_view stream)
+{
+  using T = iterator_vec_base_type<Vec2dItA>;
+
+  static_assert(is_same_floating_point<T, iterator_vec_base_type<Vec2dItB>>(),
+                "Coordinates must be the same floating point type.");
+
+  static_assert(is_same<vec_2d<T>, iterator_value_type<Vec2dItA>, iterator_value_type<Vec2dItB>>(),
+                "Inputs must be cuspatial::vec_2d<T>");
+
+  auto num_pairs = std::distance(points_geometry_offsets_first, points_geometry_offsets_last) - 1;
+
+  if (num_pairs == 0) return output_first;
+
+  std::size_t constexpr threads_per_block = 256;
+  std::size_t const num_blocks            = (num_pairs + threads_per_block - 1) / threads_per_block;
+
+  detail::pairwise_point_linestring_nearest_points_kernel<<<num_blocks,
+                                                            threads_per_block,
+                                                            0,
+                                                            stream.value()>>>(
+    points_geometry_offsets_first,
+    points_geometry_offsets_last,
+    points_first,
+    points_last,
+    linestring_geometry_offsets_first,
+    linestring_geometry_offsets_first + num_pairs + 1,
+    linestring_part_offsets_first,
+    linestring_part_offsets_last,
+    linestring_points_first,
+    linestring_points_last,
+    output_first);
+
+  CUSPATIAL_CUDA_TRY(cudaGetLastError());
+
+  return output_first + num_pairs;
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/point_linestring_nearest_points.cuh

@@ -0,0 +1,103 @@
+/*
+ * 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 <rmm/cuda_stream_view.hpp>
+
+namespace cuspatial {
+
+/**
+ * @brief Compute the nearest points and geometry ID between pairs of multipoint and
+ * multilinestring
+ *
+ * The nearest point from a test point to a linestring is a point on the linestring that has
+ * the shortest distance to the test point compared to any other points on the linestring.
+ *
+ * The nearest point from a test multipoint to a multilinestring is the nearest point that
+ * has the shortest distance in all pairs of points and linestrings.
+ *
+ * In addition, this API writes these geometry and part ID where the nearest point locates to output
+ * iterators:
+ * - The point ID indicates which point in the multipoint is the nearest point.
+ * - The linestring ID is the offset within the multilinestring that contains the nearest point.
+ * - The segment ID is the offset within the linestring of the segment that contains the nearest
+ * point. It is the same as the ID of the starting point of the segment.
+ *
+ * @tparam Cart2dItA iterator type for point array of the point element of each pair. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam Cart2dItB iterator type for point array of the linestring element of each pair. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorA iterator type for `point_geometry_offset` array. Must meet the
+ * requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorB iterator type for `linestring_geometry_offset` array. Must meet the
+ * requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorC iterator type for `linestring_part_offset` array. Must meet the
+ * requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OutputIt iterator type for output array. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ *
+ * @param point_geometry_offset_first beginning of the range of multipoint geometries of each
+ * pair
+ * @param point_geometry_offset_last end of the range of multipoint geometries of each pair
+ * @param points_first beginning of the range of point values
+ * @param points_last end of the range of the point values
+ * @param linestring_geometry_offset_first beginning of the range of offsets to the multilinestring
+ * geometry of each pair, the end range is implied by linestring_geometry_offset_first +
+ * std::distance(`point_geometry_offset_first`, `point_geometry_offset_last`)
+ * @param linestring_offsets_first beginning of the range of offsets to the starting point
+ * of each linestring
+ * @param linestring_offsets_last end of the range of offsets to the starting point
+ * of each linestring
+ * @param linestring_points_first beginning of the range of linestring points
+ * @param linestring_points_last end of the range of linestring points
+ * @param output_first A zipped-iterator of 4 outputs. The first element should be compatible
+ * with iterator_value_type<OffsetIteratorA>, stores the geometry ID of the nearest point in
+ * multipoint. The second element should be compatible with iterator_value_type<OffsetIteratorB>,
+ * stores the geometry ID of the nearest linestring. The third element should be compatible with
+ * iterator_value_type<OffsetIteratorC>, stores the part ID to the nearest segment. The forth
+ * element should be compatible with vec_2d, stores the coordinate of the nearest point on the
+ * (multi)linestring.
+ * @param stream The CUDA stream to use for device memory operations and kernel launches.
+ * @return Output iterator to the element past the last tuple computed.
+ *
+ * @pre all input iterators for coordinates must have `cuspatial::vec_2d` type,
+ * and must have the same base floating point type.
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <class Vec2dItA,
+          class Vec2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIt>
+OutputIt pairwise_point_linestring_nearest_points(
+  OffsetIteratorA points_geometry_offsets_first,
+  OffsetIteratorA points_geometry_offsets_last,
+  Vec2dItA points_first,
+  Vec2dItA points_last,
+  OffsetIteratorB linestring_geometry_offsets_first,
+  OffsetIteratorC linestring_part_offsets_first,
+  OffsetIteratorC linestring_part_offsets_last,
+  Vec2dItB linestring_points_first,
+  Vec2dItB linestring_points_last,
+  OutputIt output_first,
+  rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/point_linestring_nearest_points.cuh>

cpp/include/cuspatial/experimental/type_utils.hpp

@@ -18,6 +18,7 @@
 #include <cuspatial/traits.hpp>
 #include <cuspatial/vec_2d.hpp>
 
+#include <thrust/iterator/permutation_iterator.h>
 #include <thrust/iterator/transform_iterator.h>
 #include <thrust/iterator/transform_output_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
@@ -132,6 +133,35 @@ auto interleaved_iterator_to_vec_2d_iterator(Iter d_points_begin)
   return detail::make_counting_transform_iterator(0, interleaved_to_vec_2d<T>{d_points_begin});
 }
 
+struct strided_functor {
+  std::size_t _stride;
+  strided_functor(std::size_t stride) : _stride(stride) {}
+  auto __device__ operator()(std::size_t i) { return i * _stride; }
+};
+
+/**
+ * @brief Create an output iterator to `vec_2d` data from an iterator to an interleaved array.
+ *
+ * @tparam Iter type of iterator to interleaved data. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI], be mutable and be device-accessible.
+ * @param d_points_begin Iterator to beginning of interleaved data.
+ * @return Iterator to `vec_2d`
+ */
+template <typename Iter>
+auto vec_2d_iterator_to_output_interleaved_iterator(Iter d_points_begin)
+{
+  using T                     = typename std::iterator_traits<Iter>::value_type;
+  auto fixed_stride_2_functor = strided_functor(2);
+  auto even_positions         = thrust::make_permutation_iterator(
+    d_points_begin, detail::make_counting_transform_iterator(0, fixed_stride_2_functor));
+  auto odd_positions = thrust::make_permutation_iterator(
+    thrust::next(d_points_begin),
+    detail::make_counting_transform_iterator(0, fixed_stride_2_functor));
+  auto zipped_outputs =
+    thrust::make_zip_iterator(thrust::make_tuple(even_positions, odd_positions));
+  return thrust::make_transform_output_iterator(zipped_outputs, detail::vec_2d_to_tuple<T>());
+}
+
 /**
  * @} // end of doxygen group
  */