Add Point Linestring Distance

cpp/CMakeLists.txt

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

cpp/benchmarks/pairwise_linestring_distance.cu

@@ -25,7 +25,9 @@
 #include <rmm/device_vector.hpp>
 #include <rmm/exec_policy.hpp>
 
+#include <thrust/execution_policy.h>
 #include <thrust/iterator/counting_iterator.h>
+#include <thrust/scan.h>
 
 #include <memory>
 

cpp/benchmarks/points_in_range.cu

@@ -31,6 +31,7 @@
 #include <thrust/random/linear_congruential_engine.h>
 #include <thrust/random/normal_distribution.h>
 #include <thrust/random/uniform_int_distribution.h>
+#include <thrust/tabulate.h>
 
 #include <memory>
 

cpp/include/cuspatial/distance/point_linestring_distance.hpp

@@ -0,0 +1,65 @@
+/*
+ * 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.
+ */
+
+#include <cudf/column/column_view.hpp>
+
+namespace cuspatial {
+
+/**
+ * @brief Compute distance between pairs of points and linestrings (a.k.a. polylines)
+ *
+ * The distance between a point and a linestring is defined as the minimum distance
+ * between the point and any segment of the linestring. For each input point, this
+ * function returns the distance between the point and the corresponding linestring.
+ *
+ * The following example contains 2 pairs of points and linestrings.
+ * ```
+ * First pair:
+ * Point: (0, 0)
+ * Linestring: (0, 1) -> (1, 0) -> (2, 0)
+ *
+ * Second pair:
+ * Point: (1, 1)
+ * Linestring: (0, 0) -> (1, 1) -> (2, 0) -> (3, 0) -> (3, 1)
+ *
+ * The input of the abbove example is:
+ * points_xy: {0, 1, 0, 1}
+ * linestring_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}
+ * ```
+ *
+ * @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 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.
+ * @throws cuspatial::logic_error if the any of the point arrays have mismatched types.
+ */
+std::unique_ptr<cudf::column> pairwise_point_linestring_distance(
+  cudf::column_view const& points_xy,
+  cudf::device_span<cudf::size_type const> linestring_offsets,
+  cudf::column_view const& linestring_points_xy,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

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

@@ -24,7 +24,12 @@
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>
 
+#include <thrust/advance.h>
 #include <thrust/binary_search.h>
+#include <thrust/distance.h>
+#include <thrust/execution_policy.h>
+#include <thrust/fill.h>
+#include <thrust/memory.h>
 
 #include <cuda/atomic>
 

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

@@ -26,8 +26,12 @@
 #include <rmm/exec_policy.hpp>
 
 #include <thrust/binary_search.h>
+#include <thrust/distance.h>
+#include <thrust/execution_policy.h>
+#include <thrust/fill.h>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/transform_iterator.h>
+#include <thrust/memory.h>
 
 #include <iterator>
 #include <limits>

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

@@ -23,6 +23,8 @@
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>
 
+#include <thrust/transform.h>
+
 #include <type_traits>
 
 namespace cuspatial {

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

@@ -0,0 +1,150 @@
+/*
+ * 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/detail/utility/traits.hpp>
+#include <cuspatial/error.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/advance.h>
+#include <thrust/binary_search.h>
+#include <thrust/distance.h>
+#include <thrust/execution_policy.h>
+#include <thrust/fill.h>
+#include <thrust/memory.h>
+
+#include <iterator>
+#include <limits>
+#include <memory>
+#include <type_traits>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @internal
+ * @brief The kernel to compute point to 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.
+ *
+ * @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
+ * [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] 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
+ * 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
+ * points
+ * @param[in] linestring_points_end 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,
+                                                   Cart2dItB linestring_points_first,
+                                                   Cart2dItB linestring_points_last,
+                                                   OutputIterator distances)
+{
+  using T = iterator_vec_base_type<Cart2dItA>;
+
+  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);
+    // 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));
+    cartesian_2d<T> const a = linestring_points_first[idx];
+    cartesian_2d<T> const b = linestring_points_first[idx + 1];
+    cartesian_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)));
+  }
+}
+
+}  // 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)
+{
+  using T = detail::iterator_vec_base_type<Cart2dItA>;
+
+  static_assert(detail::is_same_floating_point<T, detail::iterator_vec_base_type<Cart2dItB>>(),
+                "Inputs must have same floating point value type.");
+
+  static_assert(detail::is_same<cartesian_2d<T>,
+                                detail::iterator_value_type<Cart2dItA>,
+                                detail::iterator_value_type<Cart2dItB>>(),
+                "Inputs must be cuspatial::cartesian_2d");
+
+  auto const num_pairs = thrust::distance(points_first, points_last);
+
+  if (num_pairs == 0) { return; }
+
+  auto const num_linestring_points =
+    thrust::distance(linestring_points_first, linestring_points_last);
+
+  thrust::fill_n(
+    rmm::exec_policy(stream), distances_first, num_pairs, std::numeric_limits<T>::max());
+
+  std::size_t constexpr threads_per_block = 256;
+  std::size_t const num_blocks =
+    (num_linestring_points + threads_per_block - 1) / threads_per_block;
+
+  detail::pairwise_point_linestring_distance<<<num_blocks, threads_per_block, 0, stream.value()>>>(
+    points_first,
+    linestring_offsets_first,
+    linestring_offsets_first + num_pairs + 1,
+    linestring_points_first,
+    linestring_points_first + num_linestring_points,
+    distances_first);
+
+  CUSPATIAL_CUDA_TRY(cudaGetLastError());
+}
+
+}  // namespace cuspatial

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

@@ -23,6 +23,10 @@
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>
 
+#include <thrust/copy.h>
+#include <thrust/count.h>
+#include <thrust/iterator/iterator_traits.h>
+
 #include <type_traits>
 
 namespace cuspatial {

cpp/include/cuspatial/experimental/point_linestring_distance.cuh

@@ -0,0 +1,67 @@
+/*
+ * 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 {
+
+/**
+ * @ingroup distance
+ * @copybrief cuspatial::pairwise_point_linestring_distance
+ *
+ * The number of distances computed is `std::distance(points_first, points_last)`.
+ *
+ * @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 OffsetIterator iterator type for 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
+ * 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 stream The CUDA stream to use for device memory operations and kernel launches.
+ *
+ * @pre all input iterators for coordinates must have `cuspatial::cartesian_2d` type.
+ * @pre all scalar types must be floating point types, and must be the same type for all input
+ * iterators and output iterators.
+ *
+ * [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);
+
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/point_linestring_distance.cuh>