Header only API for polygon-polygon distance (#1065)

This PR contains 2 major additions:
1. Range casting methods. Developer can now cast a `multipolygon_range` to a `multilinestring_range` or a `multipoint_range`. This change is included in `multipolygon_range.cuh` and `multipolygon_range_test.cu`.
2. Pairwise polygon-polygon distance. This change is separated in two parts:
    1. linestring-linestring compute kernel is refactored into algorithm/linetring_distance.cuh. 
    2. This kernel is then reused to compute polygon ring distances.

Closes #1052

Authors:
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - Mark Harris (https://github.com/harrism)

URL: #1065

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

@@ -0,0 +1,74 @@
+/*
+ * Copyright (c) 2023, 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 <cuspatial/detail/utility/device_atomics.cuh>
+#include <cuspatial/detail/utility/linestring.cuh>
+
+#include <rmm/device_uvector.hpp>
+
+#include <thrust/optional.h>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @internal
+ * @brief The kernel to compute linestring to linestring distance
+ *
+ * Each thread of the kernel computes the distance between a segment in a linestring in pair 1 to a
+ * linestring in pair 2. For a segment in pair 1, the linestring index is looked up from the offset
+ * array and mapped to the linestring in the pair 2. The segment is then computed with all segments
+ * in the corresponding linestring in pair 2. This forms a local minima of the shortest distance,
+ * which is then combined with other segment results via an atomic operation to form the global
+ * minimum distance between the linestrings.
+ *
+ * `intersects` is an optional pointer to a boolean range where the `i`th element indicates the
+ * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if
+ * set to nullopt, no distance computation will be bypassed.
+ */
+template <class MultiLinestringRange1, class MultiLinestringRange2, class OutputIt>
+__global__ void linestring_distance(MultiLinestringRange1 multilinestrings1,
+                                    MultiLinestringRange2 multilinestrings2,
+                                    thrust::optional<uint8_t*> intersects,
+                                    OutputIt distances_first)
+{
+  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);
+
+    if (intersects.has_value() && intersects.value()[geometry_idx]) {
+      distances_first[geometry_idx] = 0;
+      continue;
+    }
+
+    auto [a, b]            = multilinestrings1.segment(idx);
+    T min_distance_squared = std::numeric_limits<T>::max();
+
+    for (auto const& linestring2 : multilinestrings2[geometry_idx]) {
+      for (auto [c, d] : linestring2) {
+        min_distance_squared = min(min_distance_squared, squared_segment_distance(a, b, c, d));
+      }
+    }
+    atomicMin(&distances_first[geometry_idx], static_cast<T>(sqrt(min_distance_squared)));
+  }
+}
+
+}  // namespace detail
+}  // namespace cuspatial

cpp/include/cuspatial/detail/distance_utils.cuh

@@ -14,6 +14,7 @@
  * limitations under the License.
  */
 
+#include <cuspatial/detail/algorithm/is_point_in_polygon.cuh>
 #include <cuspatial/detail/utility/zero_data.cuh>
 #include <cuspatial/geometry/vec_2d.hpp>
 #include <cuspatial/iterator_factory.cuh>

cpp/include/cuspatial/detail/linestring_distance.cuh

@@ -16,8 +16,7 @@
 
 #pragma once
 
-#include <cuspatial/detail/utility/device_atomics.cuh>
-#include <cuspatial/detail/utility/linestring.cuh>
+#include <cuspatial/detail/algorithm/linestring_distance.cuh>
 #include <cuspatial/error.hpp>
 #include <cuspatial/geometry/vec_2d.hpp>
 #include <cuspatial/traits.hpp>
@@ -26,49 +25,12 @@
 #include <rmm/exec_policy.hpp>
 
 #include <thrust/fill.h>
+#include <thrust/optional.h>
 
 #include <limits>
 #include <type_traits>
 
 namespace cuspatial {
-namespace detail {
-
-/**
- * @internal
- * @brief The kernel to compute linestring to linestring distance
- *
- * Each thread of the kernel computes the distance between a segment in a linestring in pair 1 to a
- * linestring in pair 2. For a segment in pair 1, the linestring index is looked up from the offset
- * array and mapped to the linestring in the pair 2. The segment is then computed with all segments
- * in the corresponding linestring in pair 2. This forms a local minima of the shortest distance,
- * which is then combined with other segment results via an atomic operation to form the globally
- * minimum distance between the linestrings.
- */
-template <class MultiLinestringRange1, class MultiLinestringRange2, class OutputIt>
-__global__ void pairwise_linestring_distance_kernel(MultiLinestringRange1 multilinestrings1,
-                                                    MultiLinestringRange2 multilinestrings2,
-                                                    OutputIt distances_first)
-{
-  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);
-    T min_distance_squared  = std::numeric_limits<T>::max();
-
-    for (auto const& linestring2 : multilinestrings2[geometry_idx]) {
-      for (auto [c, d] : linestring2) {
-        min_distance_squared = min(min_distance_squared, squared_segment_distance(a, b, c, d));
-      }
-    }
-    atomicMin(&distances_first[geometry_idx], static_cast<T>(sqrt(min_distance_squared)));
-  }
-}
-
-}  // namespace detail
 
 template <class MultiLinestringRange1, class MultiLinestringRange2, class OutputIt>
 OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1,
@@ -98,8 +60,8 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1,
   std::size_t const num_blocks =
     (multilinestrings1.num_points() + threads_per_block - 1) / threads_per_block;
 
-  detail::pairwise_linestring_distance_kernel<<<num_blocks, threads_per_block, 0, stream.value()>>>(
-    multilinestrings1, multilinestrings2, distances_first);
+  detail::linestring_distance<<<num_blocks, threads_per_block, 0, stream.value()>>>(
+    multilinestrings1, multilinestrings2, thrust::nullopt, distances_first);
 
   CUSPATIAL_CUDA_TRY(cudaGetLastError());
   return distances_first + multilinestrings1.size();

cpp/include/cuspatial/detail/polygon_distance.cuh

@@ -0,0 +1,93 @@
+/*
+ * Copyright (c) 2023, 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 "distance_utils.cuh"
+#include "linestring_distance.cuh"
+
+#include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/error.hpp>
+#include <cuspatial/geometry/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/device_uvector.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/logical.h>
+#include <thrust/transform.h>
+
+#include <cstdint>
+#include <limits>
+#include <type_traits>
+
+namespace cuspatial {
+
+/**
+ * @brief Implementation of pairwise distance between two multipolygon ranges.
+ *
+ * All points in lhs and rhs are tested for intersection its corresponding pair,
+ * and if any intersection is found, the distance between the two polygons is 0.
+ * Otherwise, the distance is the minimum distance between any two segments in the
+ * multipolygon pair.
+ */
+template <class MultipolygonRangeA, class MultipolygonRangeB, class OutputIt>
+OutputIt pairwise_polygon_distance(MultipolygonRangeA lhs,
+                                   MultipolygonRangeB rhs,
+                                   OutputIt distances_first,
+                                   rmm::cuda_stream_view stream)
+{
+  using T = typename MultipolygonRangeA::element_t;
+
+  CUSPATIAL_EXPECTS(lhs.size() == rhs.size(), "Must have the same number of input rows.");
+
+  if (lhs.size() == 0) return distances_first;
+
+  auto lhs_as_multipoints = lhs.as_multipoint_range();
+  auto rhs_as_multipoints = rhs.as_multipoint_range();
+
+  auto intersects = [&]() {
+    auto lhs_in_rhs = point_polygon_intersects(lhs_as_multipoints, rhs, stream);
+    auto rhs_in_lhs = point_polygon_intersects(rhs_as_multipoints, lhs, stream);
+
+    rmm::device_uvector<uint8_t> intersects(lhs_in_rhs.size(), stream);
+    thrust::transform(rmm::exec_policy(stream),
+                      lhs_in_rhs.begin(),
+                      lhs_in_rhs.end(),
+                      rhs_in_lhs.begin(),
+                      intersects.begin(),
+                      thrust::logical_or<uint8_t>{});
+    return intersects;
+  }();
+
+  auto lhs_as_multilinestrings = lhs.as_multilinestring_range();
+  auto rhs_as_multilinestrings = rhs.as_multilinestring_range();
+
+  thrust::fill(rmm::exec_policy(stream),
+               distances_first,
+               distances_first + lhs.size(),
+               std::numeric_limits<T>::max());
+
+  auto [threads_per_block, num_blocks] = grid_1d(lhs.num_points());
+
+  detail::linestring_distance<<<num_blocks, threads_per_block, 0, stream.value()>>>(
+    lhs_as_multilinestrings, rhs_as_multilinestrings, intersects.begin(), distances_first);
+
+  CUSPATIAL_CUDA_TRY(cudaGetLastError());
+  return distances_first + lhs.size();
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/detail/range/multipolygon_range.cuh

@@ -23,6 +23,8 @@
 #include <cuspatial/geometry/segment.cuh>
 #include <cuspatial/geometry/vec_2d.hpp>
 #include <cuspatial/geometry_collection/multipolygon_ref.cuh>
+#include <cuspatial/range/multilinestring_range.cuh>
+#include <cuspatial/range/multipoint_range.cuh>
 #include <cuspatial/traits.hpp>
 
 #include <thrust/binary_search.h>
@@ -441,4 +443,39 @@ multipolygon_range<GeometryIterator, PartIterator, RingIterator, VecIterator>::
   return point_idx == _ring_begin[_part_begin[_geometry_begin[geometry_idx]]];
 }
 
+template <typename GeometryIterator,
+          typename PartIterator,
+          typename RingIterator,
+          typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto
+multipolygon_range<GeometryIterator, PartIterator, RingIterator, VecIterator>::as_multipoint_range()
+{
+  auto multipoint_geometry_it = thrust::make_permutation_iterator(
+    _ring_begin, thrust::make_permutation_iterator(_part_begin, _geometry_begin));
+
+  return multipoint_range{multipoint_geometry_it,
+                          multipoint_geometry_it + thrust::distance(_geometry_begin, _geometry_end),
+                          _point_begin,
+                          _point_end};
+}
+
+template <typename GeometryIterator,
+          typename PartIterator,
+          typename RingIterator,
+          typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto
+multipolygon_range<GeometryIterator, PartIterator, RingIterator, VecIterator>::
+  as_multilinestring_range()
+{
+  auto multilinestring_geometry_it =
+    thrust::make_permutation_iterator(_part_begin, _geometry_begin);
+  return multilinestring_range{
+    multilinestring_geometry_it,
+    multilinestring_geometry_it + thrust::distance(_geometry_begin, _geometry_end),
+    _ring_begin,
+    _ring_end,
+    _point_begin,
+    _point_end};
+}
+
 }  // namespace cuspatial

cpp/include/cuspatial/polygon_distance.cuh

@@ -0,0 +1,47 @@
+/*
+ * Copyright (c) 2023, 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
+ * @brief Computes pairwise multipolygon to multipolygon distance
+ *
+ * @tparam MultiPolygonRangeA An instance of template type `multipolygon_range`
+ * @tparam MultiPolygonRangeB An instance of template type `multipolygon_range`
+ * @tparam OutputIt iterator type for output array. Must meet the requirements of [LRAI](LinkLRAI).
+ * Must be an iterator to type convertible from floating points.
+ *
+ * @param lhs The first multipolygon range to compute distance from
+ * @param rhs The second multipolygon range to compute distance to
+ * @param stream The CUDA stream on which to perform computations
+ * @return Output Iterator past the last distance computed
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <class MultipolygonRangeA, class MultipolygonRangeB, class OutputIt>
+OutputIt pairwise_polygon_distance(MultipolygonRangeA lhs,
+                                   MultipolygonRangeB rhs,
+                                   OutputIt distances_first,
+                                   rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+}  // namespace cuspatial
+
+#include <cuspatial/detail/polygon_distance.cuh>

cpp/include/cuspatial/range/multipolygon_range.cuh

@@ -62,7 +62,9 @@ class multipolygon_range {
   using ring_it_t     = RingIterator;
   using point_it_t    = VecIterator;
   using point_t       = iterator_value_type<VecIterator>;
-  using element_t     = iterator_vec_base_type<VecIterator>;
+
+  using index_t   = iterator_value_type<GeometryIterator>;
+  using element_t = iterator_vec_base_type<VecIterator>;
 
   int64_t static constexpr INVALID_INDEX = -1;
 
@@ -180,6 +182,15 @@ class multipolygon_range {
   /// Returns an iterator to the end of the segment
   CUSPATIAL_HOST_DEVICE auto segment_end();
 
+  /// Range Casting
+
+  /// Cast the range of multipolygons as a range of multipoints, ignoring all edge connections and
+  /// ring relationships.
+  CUSPATIAL_HOST_DEVICE auto as_multipoint_range();
+
+  /// Cast the range of multipolygons as a range of multilinestrings, ignoring ring relationships.
+  CUSPATIAL_HOST_DEVICE auto as_multilinestring_range();
+
  protected:
   GeometryIterator _geometry_begin;
   GeometryIterator _geometry_end;

cpp/include/cuspatial_test/base_fixture.hpp

@@ -66,7 +66,9 @@ class BaseFixture : public RMMResourceMixin, public ::testing::Test {};
  * class MyTest : public cuspatial::test::BaseFixtureWithParam {};
  *
  * TEST_P(MyTest, TestParamterGet) {
- *  auto [a, b, c] = GetParam();
+ *  auto a = std::get<0>(GetParam());
+ *  auto b = std::get<1>(GetParam());
+ *  auto c = std::get<2>(GetParam());
  *  ...
  * }
  *

cpp/include/cuspatial_test/test_util.cuh

@@ -100,7 +100,14 @@ void print_device_range(Iter begin,
 }
 
 /**
- * @brief
+ * @brief Print a device vector.
+ *
+ * @note Copies the device vector to host before printing.
+ *
+ * @tparam Vector The device vector type
+ * @param vec The device vector to print
+ * @param pre String to print before the device vector
+ * @param post String to print after the device vector
  */
 template <typename Vector>
 void print_device_vector(Vector const& vec, std::string_view pre = "", std::string_view post = "\n")