Pairwise Point to Point Distance; Rename Folder distances to `dista…

…nce` (#558)

Implements 2D point-point L2 distance. Also renamed folder `distances` into `distance`

Contributes to #231

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

Approvers:
  - H. Thomson Comer (https://github.com/thomcom)
  - Robert Maynard (https://github.com/robertmaynard)
  - Mark Harris (https://github.com/harrism)

URL: #558

cpp/CMakeLists.txt

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

cpp/benchmarks/hausdorff_benchmark.cpp

@@ -14,7 +14,7 @@
  * limitations under the License.
  */
 
-#include <cuspatial/hausdorff.hpp>
+#include <cuspatial/distance/hausdorff.hpp>
 
 #include <benchmarks/fixture/benchmark_fixture.hpp>
 #include <benchmarks/synchronization/synchronization.hpp>

cpp/include/cuspatial/distance/point_distance.hpp

@@ -0,0 +1,41 @@
+/*
+ * 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 <cudf/column/column_view.hpp>
+
+namespace cuspatial {
+
+/**
+ * @ingroup distance
+ * @brief Compute pairwise point-to-point Cartesian distance
+ *
+ * @param points1_x Column of x-coordinates of the first point in each pair
+ * @param points1_y Column of y-coordinates of the first point in each pair
+ * @param points2_x Column of x-coordinates of the second point in each pair
+ * @param points2_y Column of y-coordinates of the second point in each pair
+ * @param stream The CUDA stream to use for device memory operations and kernel launches
+ * @return Column of distances between each pair of input points
+ */
+std::unique_ptr<cudf::column> pairwise_point_distance(
+  cudf::column_view const& points1_x,
+  cudf::column_view const& points1_y,
+  cudf::column_view const& points2_x,
+  cudf::column_view const& points2_y,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

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

@@ -152,15 +152,13 @@ void pairwise_linestring_distance(OffsetIterator linestring1_offsets_first,
                               typename std::iterator_traits<OutputIt>::value_type>(),
     "Inputs and output must have floating point value type.");
 
-  static_assert(detail::is_same<T,
-                                typename std::iterator_traits<Cart2dItB>::value_type::value_type,
-                                typename std::iterator_traits<OutputIt>::value_type>(),
+  static_assert(detail::is_same<T, typename std::iterator_traits<OutputIt>::value_type>(),
                 "Inputs and output must have the same value type.");
 
   static_assert(detail::is_same<cartesian_2d<T>,
                                 typename std::iterator_traits<Cart2dItA>::value_type,
                                 typename std::iterator_traits<Cart2dItB>::value_type>(),
-                "Inputs must be cuspatial::cartesian_2d");
+                "All input types must be cuspatial::cartesian_2d with the same value type");
 
   auto const num_linestring_pairs =
     thrust::distance(linestring1_offsets_first, linestring1_offsets_last);

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

@@ -0,0 +1,64 @@
+/*
+ * 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/traits.hpp>
+#include <cuspatial/error.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <type_traits>
+
+namespace cuspatial {
+
+template <class Cart2dItA, class Cart2dItB, class OutputIt>
+OutputIt pairwise_point_distance(Cart2dItA points1_first,
+                                 Cart2dItA points1_last,
+                                 Cart2dItB points2_first,
+                                 OutputIt distances_first,
+                                 rmm::cuda_stream_view stream)
+{
+  using T = typename std::iterator_traits<Cart2dItA>::value_type::value_type;
+
+  static_assert(
+    detail::is_floating_point<T,
+                              typename std::iterator_traits<Cart2dItB>::value_type::value_type,
+                              typename std::iterator_traits<OutputIt>::value_type>(),
+    "Inputs and output must be floating point types.");
+
+  static_assert(detail::is_same<T, typename std::iterator_traits<OutputIt>::value_type>(),
+                "Inputs and output must be the same types.");
+
+  static_assert(detail::is_same<cartesian_2d<T>,
+                                typename std::iterator_traits<Cart2dItA>::value_type,
+                                typename std::iterator_traits<Cart2dItB>::value_type>(),
+                "All Input types must be cuspatial::cartesian_2d with the same value type");
+
+  return thrust::transform(rmm::exec_policy(stream),
+                           points1_first,
+                           points1_last,
+                           points2_first,
+                           distances_first,
+                           [] __device__(auto p1, auto p2) {
+                             auto v = p1 - p2;
+                             return std::sqrt(dot(v, v));
+                           });
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/linestring_distance.cuh

@@ -21,6 +21,7 @@
 namespace cuspatial {
 
 /**
+ * @ingroup distance
  * @copybrief cuspatial::pairwise_linestring_distance
  *
  * The shortest distance between two linestrings is defined as the shortest distance

cpp/include/cuspatial/experimental/point_distance.cuh

@@ -0,0 +1,57 @@
+/*
+ * 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_distance
+ *
+ * @tparam Cart2dItA iterator type for point array of the first point of each pair. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam Cart2dItB iterator type for point array of the second point of each pair. 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], be mutable and be device-accessible.
+ *
+ * @param points1_first beginning of range of the first point of each pair
+ * @param points1_last end of range of the first point of each pair
+ * @param points2_first beginning of range of the second point of each pair
+ * @param distances_first beginning iterator to output
+ * @param stream The CUDA stream to use for device memory operations and kernel launches
+ * @return Output iterator to one past the last element in the output range
+ *
+ * @pre all input iterators for coordinates must have a `value_type` of `cuspatial::cartesian_2d`.
+ * @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 OutputIt>
+OutputIt pairwise_point_distance(Cart2dItA points1_first,
+                                 Cart2dItA points1_last,
+                                 Cart2dItB points2_first,
+                                 OutputIt distances_first,
+                                 rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/point_distance.cuh>

cpp/include/doxygen_groups.h

@@ -48,6 +48,7 @@
  *      @{
  *          @brief Distance computation APIs
  *
+ *          @file point_distance.hpp
  *          @file linestring_distance.hpp
  *          @file linestring_distance.cuh
  *          @file hausdorff.hpp

cpp/src/spatial/point_distance.cu

@@ -0,0 +1,113 @@
+/*
+ * 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 <cuspatial/error.hpp>
+#include <cuspatial/experimental/point_distance.cuh>
+#include <cuspatial/experimental/type_utils.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <cudf/column/column_factories.hpp>
+#include <cudf/column/column_view.hpp>
+#include <cudf/copying.hpp>
+#include <cudf/utilities/type_dispatcher.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <memory>
+#include <type_traits>
+
+namespace cuspatial {
+namespace detail {
+struct pairwise_point_distance_functor {
+  template <typename T>
+  std::enable_if_t<std::is_floating_point<T>::value, std::unique_ptr<cudf::column>> operator()(
+    cudf::column_view const& points1_x,
+    cudf::column_view const& points1_y,
+    cudf::column_view const& points2_x,
+    cudf::column_view const& points2_y,
+    rmm::cuda_stream_view stream,
+    rmm::mr::device_memory_resource* mr)
+  {
+    auto distances = cudf::make_numeric_column(cudf::data_type{cudf::type_to_id<T>()},
+                                               points1_x.size(),
+                                               cudf::mask_state::UNALLOCATED,
+                                               stream,
+                                               mr);
+
+    auto points1_it = make_cartesian_2d_iterator(points1_x.begin<T>(), points1_y.begin<T>());
+    auto points2_it = make_cartesian_2d_iterator(points2_x.begin<T>(), points2_y.begin<T>());
+
+    pairwise_point_distance(points1_it,
+                            points1_it + points1_x.size(),
+                            points2_it,
+                            distances->mutable_view().begin<T>(),
+                            stream);
+
+    return distances;
+  }
+
+  template <typename T, typename... Args>
+  std::enable_if_t<not std::is_floating_point<T>::value, std::unique_ptr<cudf::column>> operator()(
+    Args&&...)
+  {
+    CUSPATIAL_FAIL("Point distances only supports floating point coordinates.");
+  }
+};
+
+std::unique_ptr<cudf::column> pairwise_point_distance(cudf::column_view const& points1_x,
+                                                      cudf::column_view const& points1_y,
+                                                      cudf::column_view const& points2_x,
+                                                      cudf::column_view const& points2_y,
+                                                      rmm::cuda_stream_view stream,
+                                                      rmm::mr::device_memory_resource* mr)
+{
+  CUSPATIAL_EXPECTS(points1_x.size() == points1_y.size() and
+                      points2_x.size() == points2_y.size() and points1_x.size() == points2_x.size(),
+                    "Mismatch number of coordinate or number of points.");
+
+  CUSPATIAL_EXPECTS(points1_x.type() == points1_y.type() and
+                      points2_x.type() == points2_y.type() and points1_x.type() == points2_x.type(),
+                    "The types of point coordinates arrays mismatch.");
+  CUSPATIAL_EXPECTS(not points1_x.has_nulls() and not points1_y.has_nulls() and
+                      not points2_x.has_nulls() and not points2_y.has_nulls(),
+                    "The coordinate columns cannot have nulls.");
+
+  if (points1_x.size() == 0) { return cudf::empty_like(points1_x); }
+
+  return cudf::type_dispatcher(points1_x.type(),
+                               pairwise_point_distance_functor{},
+                               points1_x,
+                               points1_y,
+                               points2_x,
+                               points2_y,
+                               stream,
+                               mr);
+}
+
+}  // namespace detail
+
+std::unique_ptr<cudf::column> pairwise_point_distance(cudf::column_view const& points1_x,
+                                                      cudf::column_view const& points1_y,
+                                                      cudf::column_view const& points2_x,
+                                                      cudf::column_view const& points2_y,
+                                                      rmm::mr::device_memory_resource* mr)
+{
+  return detail::pairwise_point_distance(
+    points1_x, points1_y, points2_x, points2_y, rmm::cuda_stream_default, mr);
+}
+
+}  // namespace cuspatial

cpp/tests/CMakeLists.txt

@@ -77,6 +77,9 @@ ConfigureTest(POLYLINE_BBOX_TEST
 ConfigureTest(POLYGON_BBOX_TEST
     spatial/polygon_bbox_test.cu)
 
+ConfigureTest(POINT_DISTANCE_TEST
+    spatial/point_distance_test.cpp)
+
 ConfigureTest(LINESTRING_DISTANCE_TEST
     spatial/linestring_distance_test.cpp spatial/linestring_distance_test_medium.cpp)
 
@@ -105,6 +108,9 @@ ConfigureTest(SPATIAL_WINDOW_POINT_TEST
 ConfigureTest(HAVERSINE_TEST_EXP
     experimental/spatial/haversine_test.cu)
 
+ConfigureTest(POINT_DISTANCE_TEST_EXP
+    experimental/spatial/point_distance_test.cu)
+
 ConfigureTest(HAUSDORFF_TEST_EXP
     experimental/spatial/hausdorff_test.cu)