Add Multi-Geometry support to point_linestring_distance and build python bindings

cpp/include/cuspatial/distance/point_linestring_distance.hpp

@@ -35,30 +35,65 @@ namespace cuspatial {
  * Point: (1, 1)
  * Linestring: (0, 0) -> (1, 1) -> (2, 0) -> (3, 0) -> (3, 1)
  *
- * The input of the abbove example is:
+ * The input of the above example is:
+ * multipoint_geometry_offsets: nullopt
  * points_xy: {0, 1, 0, 1}
- * linestring_offsets: {0, 3, 8}
+ * multilinestring_geometry_offsets: nullopt
+ * linestring_part_offsets: {0, 3, 8}
  * linestring_xy: {0, 1, 1, 0, 2, 0, 0, 0, 1, 1, 2, 0, 3, 0, 3, 1}
  *
  * Result: {sqrt(2)/2, 0}
  * ```
  *
+ * The following example contains 3 pairs of MultiPoint and MultiLinestring.
+ * ```
+ * First pair:
+ * MultiPoint: (0, 1)
+ * MultiLinestring: (0, -1) -> (-2, -3), (-4, -5) -> (-5, -6)
+ *
+ * Second pair:
+ * MultiPoint: (2, 3), (4, 5)
+ * MultiLinestring: (7, 8) -> (8, 9)
+ *
+ * Third pair:
+ * MultiPoint: (6, 7), (8, 9)
+ * MultiLinestring: (9, 10) -> (10, 11)
+
+ * The input of the above example is:
+ * multipoint_geometry_offsets: {0, 1, 3, 5}
+ * points_xy: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
+ * multilinestring_geometry_offsets: {0, 2, 3, 5}
+ * linestring_part_offsets: {0, 2, 4, 6, 8}
+ * linestring_points_xy: {0, -1, -2, -3, -4, -5, -5, -6, 7, 8, 8, 9, 9, 10, 10 ,11}
+ *
+ * Result: {2.0, 4.24264, 1.41421}
+ * ```
+ *
+ * @param multipoint_geometry_offsets Beginning and ending indices to each geometry in the
+ * multi-point
  * @param points_xy Interleaved x, y-coordinates of points
- * @param linestring_offsets Beginning and ending indices for each linestring in the `linestring_x`
- * and `linestring_y` arrays. Must satisfy `linestring_offsets.size() + 1 == points_xy.size()`.
+ * @param multilinestring_geometry_offsets Beginning and ending indices to each geometry in the
+ * multi-linestring
+ * @param linestring_part_offsets Beginning and ending indices for each linestring in the point
+ * array. Because the coordinates are interleaved, the actual starting position for the coordinate
+ * of linestring `i` is `2*linestring_part_offsets[i]`.
  * @param linestring_points_xy Interleaved x, y-coordinates of linestring points.
  * @param mr Device memory resource used to allocate the returned column.
  * @return A column containing the distance between each pair of corresponding points and
  * linestrings.
  *
- * @throws cuspatial::logic_error if the number of points and linestrings do not match.
- * @throws cuspatial::logic_error if there is a size mismatch between the x- and y-coordinates of
- * the points or linestring points.
+ * @note Any optional geometry indices, if is `nullopt`, implies the underlying geometry contains
+ * only one component. Otherwise, it contains multiple components.
+ *
+ * @throws cuspatial::logic_error if the number of (multi)points and (multi)linestrings do not
+ * match.
  * @throws cuspatial::logic_error if the any of the point arrays have mismatched types.
  */
 std::unique_ptr<cudf::column> pairwise_point_linestring_distance(
+  std::optional<cudf::device_span<cudf::size_type const>> multipoint_geometry_offsets,
   cudf::column_view const& points_xy,
-  cudf::device_span<cudf::size_type const> linestring_offsets,
+  std::optional<cudf::device_span<cudf::size_type const>> multilinestring_geometry_offsets,
+  cudf::device_span<cudf::size_type const> linestring_part_offsets,
   cudf::column_view const& linestring_points_xy,
   rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
 

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

@@ -42,38 +42,60 @@ namespace detail {
 
 /**
  * @internal
- * @brief The kernel to compute point to linestring distance
+ * @brief The kernel to compute multi-point to multi-linestring distance
  *
  * Each thread computes the distance between a line segment in the linestring and the
- * corresponding point in the pair. The shortest distance is computed in the output
- * array via an atomic operation.
+ * corresponding multi-point part in the pair. The shortest distance is computed in the
+ * output array via an atomic operation.
  *
  * @tparam Cart2dItA Iterator to 2d cartesian coordinates. Must meet requirements of
  * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
  * @tparam Cart2dItB Iterator to 2d cartesian coordinates. Must meet requirements of
  * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
- * @tparam OffsetIterator Iterator to linestring offsets. Must meet requirements of
+ * @tparam OffsetIteratorA Iterator to offsets. Must meet requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorB Iterator to offsets. Must meet requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorC Iterator to offsets. Must meet requirements of
  * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
  * @tparam OutputIterator Iterator to output distances. Must meet requirements of
  * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible and mutable.
  *
+ * @param[in] point_geometry_offset_first Iterator to the beginning of the range of the multipoint
+ * parts
+ * @param[in] point_geometry_offset_last Iterator to the end of the range of the multipoint parts
  * @param[in] points_first Iterator to the beginning of the range of the points
  * @param[in] points_last  Iterator to the end of the range of the points
- * @param[in] linestring_offsets_begin Iterator to the beginning of the range of the linestring
+ * @param[in] linestring_geometry_offset_first Iterator to the beginning of the range of the
+ * linestring parts
+ * @param[in] linestring_geometry_offset_last Iterator to the end of the range of the linestring
+ * parts
+ * @param[in] linestring_part_offsets_first Iterator to the beginning of the range of the linestring
  * offsets
- * @param[in] linestring_offsets_end Iterator to the end of the range of the linestring offsets
- * @param[in] linestring_points_begin Iterator to the beginning of the range of the linestring
+ * @param[in] linestring_part_offsets_last Iterator to the beginning of the range of the linestring
+ * offsets
+ * @param[in] linestring_points_first Iterator to the beginning of the range of the linestring
  * points
- * @param[in] linestring_points_end Iterator to the end of the range of the linestring points
+ * @param[in] linestring_points_last Iterator to the end of the range of the linestring points
  * @param[out] distances Iterator to the output range of shortest distances between pairs.
  *
  * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
  * "LegacyRandomAccessIterator"
  */
-template <typename Cart2dItA, typename Cart2dItB, typename OffsetIterator, typename OutputIterator>
-void __global__ pairwise_point_linestring_distance(Cart2dItA points_first,
-                                                   OffsetIterator linestring_offsets_first,
-                                                   OffsetIterator linestring_offsets_last,
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIterator>
+void __global__ pairwise_point_linestring_distance(OffsetIteratorA point_geometry_offset_first,
+                                                   OffsetIteratorA point_geometry_offset_last,
+                                                   Cart2dItA points_first,
+                                                   Cart2dItA points_last,
+                                                   OffsetIteratorB linestring_geometry_offset_first,
+                                                   OffsetIteratorB linestring_geometry_offset_last,
+                                                   OffsetIteratorC linestring_part_offsets_first,
+                                                   OffsetIteratorC linestring_part_offsets_last,
                                                    Cart2dItB linestring_points_first,
                                                    Cart2dItB linestring_points_last,
                                                    OutputIterator distances)
@@ -83,34 +105,69 @@ void __global__ pairwise_point_linestring_distance(Cart2dItA points_first,
   for (auto idx = threadIdx.x + blockIdx.x * blockDim.x;
        idx < std::distance(linestring_points_first, thrust::prev(linestring_points_last));
        idx += gridDim.x * blockDim.x) {
-    auto offsets_iter =
-      thrust::upper_bound(thrust::seq, linestring_offsets_first, linestring_offsets_last, idx);
+    // Search from the part offsets array to determine the part idx of current linestring point
+    auto linestring_part_offsets_iter = thrust::upper_bound(
+      thrust::seq, linestring_part_offsets_first, linestring_part_offsets_last, idx);
+
     // Pointer to the last point in the linestring, skip iteration.
     // Note that the last point for the last linestring is guarded by the grid-stride loop.
-    if (offsets_iter != linestring_offsets_last and *offsets_iter - 1 == idx) { continue; }
-
-    auto point_idx    = thrust::distance(linestring_offsets_first, thrust::prev(offsets_iter));
-    vec_2d<T> const a = linestring_points_first[idx];
-    vec_2d<T> const b = linestring_points_first[idx + 1];
-    vec_2d<T> const c = points_first[point_idx];
-
-    auto const distance_squared = point_to_segment_distance_squared(c, a, b);
-
-    atomicMin(&thrust::raw_reference_cast(*(distances + point_idx)),
-              static_cast<T>(std::sqrt(distance_squared)));
+    if (linestring_part_offsets_iter != linestring_part_offsets_last &&
+        *linestring_part_offsets_iter - 1 == idx) {
+      continue;
+    }
+
+    auto part_offsets_idx =
+      thrust::distance(linestring_part_offsets_first, thrust::prev(linestring_part_offsets_iter));
+
+    // Search from the linestring geometry offsets array to determine the geometry idx of current
+    // linestring point
+    auto geometry_offsets_iter = thrust::upper_bound(thrust::seq,
+                                                     linestring_geometry_offset_first,
+                                                     linestring_geometry_offset_last,
+                                                     part_offsets_idx);
+    // geometry_idx is also the index to corresponding multipoint in the pair
+    auto geometry_idx =
+      thrust::distance(linestring_geometry_offset_first, thrust::prev(geometry_offsets_iter));
+
+    // Reduce the minimum distance between different parts of the multi-point.
+    vec_2d<T> const a      = linestring_points_first[idx];
+    vec_2d<T> const b      = linestring_points_first[idx + 1];
+    T min_distance_squared = std::numeric_limits<T>::max();
+
+    for (auto point_idx = point_geometry_offset_first[geometry_idx];
+         point_idx < point_geometry_offset_first[geometry_idx + 1];
+         point_idx++) {
+      vec_2d<T> const c = points_first[point_idx];
+
+      // TODO: reduce redundant computation only related to `a`, `b` in this helper.
+      auto const distance_squared = point_to_segment_distance_squared(c, a, b);
+      min_distance_squared        = std::min(distance_squared, min_distance_squared);
+    }
+
+    atomicMin(&thrust::raw_reference_cast(*(distances + geometry_idx)),
+              static_cast<T>(std::sqrt(min_distance_squared)));
   }
 }
 
 }  // namespace detail
 
-template <class Cart2dItA, class Cart2dItB, class OffsetIterator, class OutputIt>
-void pairwise_point_linestring_distance(Cart2dItA points_first,
-                                        Cart2dItA points_last,
-                                        OffsetIterator linestring_offsets_first,
-                                        Cart2dItB linestring_points_first,
-                                        Cart2dItB linestring_points_last,
-                                        OutputIt distances_first,
-                                        rmm::cuda_stream_view stream)
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIt>
+OutputIt pairwise_point_linestring_distance(OffsetIteratorA point_geometry_offset_first,
+                                            OffsetIteratorA point_geometry_offset_last,
+                                            Cart2dItA points_first,
+                                            Cart2dItA points_last,
+                                            OffsetIteratorB linestring_geometry_offset_first,
+                                            OffsetIteratorC linestring_part_offsets_first,
+                                            OffsetIteratorC linestring_part_offsets_last,
+                                            Cart2dItB linestring_points_first,
+                                            Cart2dItB linestring_points_last,
+                                            OutputIt distances_first,
+                                            rmm::cuda_stream_view stream)
 {
   using T = detail::iterator_vec_base_type<Cart2dItA>;
 
@@ -120,11 +177,12 @@ void pairwise_point_linestring_distance(Cart2dItA points_first,
   static_assert(detail::is_same<vec_2d<T>,
                                 detail::iterator_value_type<Cart2dItA>,
                                 detail::iterator_value_type<Cart2dItB>>(),
-                "Inputs must be cuspatial::vec_2d<T>");
+                "Inputs must be cuspatial::vec_2d");
 
-  auto const num_pairs = thrust::distance(points_first, points_last);
+  auto const num_pairs =
+    thrust::distance(point_geometry_offset_first, point_geometry_offset_last) - 1;
 
-  if (num_pairs == 0) { return; }
+  if (num_pairs == 0) { return distances_first; }
 
   auto const num_linestring_points =
     thrust::distance(linestring_points_first, linestring_points_last);
@@ -137,14 +195,21 @@ void pairwise_point_linestring_distance(Cart2dItA points_first,
     (num_linestring_points + threads_per_block - 1) / threads_per_block;
 
   detail::pairwise_point_linestring_distance<<<num_blocks, threads_per_block, 0, stream.value()>>>(
+    point_geometry_offset_first,
+    point_geometry_offset_last,
     points_first,
-    linestring_offsets_first,
-    linestring_offsets_first + num_pairs + 1,
+    points_last,
+    linestring_geometry_offset_first,
+    linestring_geometry_offset_first + num_pairs + 1,
+    linestring_part_offsets_first,
+    linestring_part_offsets_last,
     linestring_points_first,
-    linestring_points_first + num_linestring_points,
+    linestring_points_last,
     distances_first);
 
   CUSPATIAL_CUDA_TRY(cudaGetLastError());
+
+  return distances_first + num_pairs;
 }
 
 }  // namespace cuspatial

cpp/include/cuspatial/experimental/point_linestring_distance.cuh

@@ -30,20 +30,30 @@ namespace cuspatial {
  * 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 OffsetIterator iterator type for offset array. 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 points_first beginning of the range of points making up the first element of each
- * pair
- * @param points_last end of the range of the points making up the first element of each pair
- * @param linestring_offsets_first beginning of the range of the offsets to the linestring element
- * of each pair
- * @param linestring_points_first beginning of the range of points of the linestring element of each
+ * @param point_geometry_offset_first beginning of the range of multipoint geometries of each
  * pair
- * @param linestring_points_last end of the range of points of the linestring element of each pair
- * @param distances_first beginning the output range of distances
+ * @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 distances_first beginning of the output range of distances
  * @param stream The CUDA stream to use for device memory operations and kernel launches.
  *
  * @pre all input iterators for coordinates must have `cuspatial::vec_2d` type.
@@ -53,14 +63,24 @@ namespace cuspatial {
  * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
  * "LegacyRandomAccessIterator"
  */
-template <class Cart2dItA, class Cart2dItB, class OffsetIterator, class OutputIt>
-void pairwise_point_linestring_distance(Cart2dItA points_first,
-                                        Cart2dItA points_last,
-                                        OffsetIterator linestring_offsets_first,
-                                        Cart2dItB linestring_points_first,
-                                        Cart2dItB linestring_points_last,
-                                        OutputIt distances_first,
-                                        rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OffsetIteratorC,
+          class OutputIt>
+OutputIt pairwise_point_linestring_distance(
+  OffsetIteratorA point_geometry_offset_first,
+  OffsetIteratorA point_geometry_offset_last,
+  Cart2dItA points_first,
+  Cart2dItA points_last,
+  OffsetIteratorB linestring_geometry_offset_first,
+  OffsetIteratorC linestring_part_offsets_first,
+  OffsetIteratorC linestring_part_offsets_last,
+  Cart2dItB linestring_points_first,
+  Cart2dItB linestring_points_last,
+  OutputIt distances_first,
+  rmm::cuda_stream_view stream = rmm::cuda_stream_default);
 
 }  // namespace cuspatial