C++ pairwise linestring distance

cpp/CMakeLists.txt

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

cpp/include/cuspatial/distances/linestring_distance.hpp

@@ -0,0 +1,109 @@
+/*
+ * 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>
+#include <cudf/utilities/span.hpp>
+
+#include <memory>
+
+namespace cuspatial {
+
+/**
+ * @brief Compute shortest distance between pairs of linestrings (a.k.a. polylines)
+ *
+ * The shortest distance between two linestrings is defined as the shortest distance
+ * between all pairs of segments of the two linestrings. If any of the segments intersect,
+ * the distance is 0.
+ *
+ * The following example contains 4 pairs of linestrings.
+ *
+ * First pair:
+ * (0, 1) -> (1, 0) -> (-1, 0)
+ * (1, 1) -> (2, 1) -> (2, 0) -> (3, 0)
+ *
+ *     |
+ *     *   #---#
+ *     | \     |
+ * ----O---*---#---#
+ *     | /
+ *     *
+ *     |
+ *
+ * The shortest distance between the two linestrings is the distance
+ * from point (1, 1) to segment (0, 1) -> (1, 0), which is sqrt(2)/2.
+ *
+ * Second pair:
+ *
+ * (0, 0) -> (0, 1)
+ * (1, 0) -> (1, 1) -> (1, 2)
+ *
+ * These linestrings are parallel. Their distance is 1 (point (0, 0) to point (1, 0)).
+ *
+ * Third pair:
+ *
+ * (0, 0) -> (2, 2) -> (-2, 0)
+ * (2, 0) -> (0, 2)
+ *
+ * These linestrings intersect, so their distance is 0.
+ *
+ * Forth pair:
+ *
+ * (2, 2) -> (-2, -2)
+ * (1, 1) -> (5, 5) -> (10, 0)
+ *
+ * These linestrings contain colinear and overlapping sections, so
+ * their distance is 0.
+ *
+ * The input of above example is:
+ * linestring1_offsets:  {0, 3, 5, 8}
+ * linestring1_points_x: {0, 1, -1, 0, 0, 0, 2, -2, 2, -2}
+ * linestring1_points_y: {1, 0, 0, 0, 1, 0, 2, 0, 2, -2}
+ * linestring2_offsets:  {0, 4, 7, 9}
+ * linestring2_points_x: {1, 2, 2, 3, 1, 1, 1, 2, 0, 1, 5, 10}
+ * linestring2_points_y: {1, 1, 0, 0, 0, 1, 2, 0, 2, 1, 5, 0}
+ *
+ * Result: {sqrt(2.0)/2, 1, 0, 0}
+ *
+ * @param linestring1_offsets Indices of the first point of the first linestring of each pair.
+ * @param linestring1_points_x x-components of points in the first linestring of each pair.
+ * @param linestring1_points_y y-component of points in the first linestring of each pair.
+ * @param linestring2_offsets Indices of the first point of the second linestring of each pair.
+ * @param linestring2_points_x x-component of points in the first linestring of each pair.
+ * @param linestring2_points_y y-component of points in the first linestring of each pair.
+ * @param mr Device memory resource used to allocate the returned column's device memory.
+ * @return A column of shortest distances between each pair of linestrings.
+ *
+ * @note If any of the linestring contains less than 2 points, the behavior is undefined.
+ *
+ * @throw cuspatial::logic_error if `linestring1_offsets.size() != linestring2_offsets.size()`
+ * @throw cuspatial::logic_error if there is a size mismatch between the x- and y-coordinates of the
+ * linestring points.
+ * @throw cuspatial::logic_error if any of the point arrays have mismatched types.
+ * @throw cuspatial::logic_error if any linestring has fewer than 2 points.
+ *
+ */
+std::unique_ptr<cudf::column> pairwise_linestring_distance(
+  cudf::device_span<cudf::size_type const> linestring1_offsets,
+  cudf::column_view const& linestring1_points_x,
+  cudf::column_view const& linestring1_points_y,
+  cudf::device_span<cudf::size_type const> linestring2_offsets,
+  cudf::column_view const& linestring2_points_x,
+  cudf::column_view const& linestring2_points_y,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource());
+
+}  // namespace cuspatial

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

@@ -19,6 +19,7 @@
 #include <cuspatial/constants.hpp>
 #include <cuspatial/error.hpp>
 #include <cuspatial/types.hpp>
+#include <cuspatial/utility/vec_2d.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>

cpp/include/cuspatial/types.hpp

@@ -18,31 +18,13 @@
 
 #include <cstdint>
 
-namespace cuspatial {
+#ifdef __CUDACC__
+#define CUSPATIAL_HOST_DEVICE __host__ __device__
+#else
+#define CUSPATIAL_HOST_DEVICE
+#endif
 
-/**
- * @brief A 2D vector
- *
- * Used in cuspatial for both Longitude/Latitude (LonLat) coordinate pairs and Cartesian (X/Y)
- * coordinate pairs. For LonLat pairs, the `x` member represents Longitude, and `y` represents
- * Latitude.
- *
- * @tparam T the base type for the coordinates
- */
-template <typename T>
-struct alignas(2 * sizeof(T)) vec_2d {
-  using value_type = T;
-  value_type x;
-  value_type y;
-};
-
-template <typename T>
-struct alignas(2 * sizeof(T)) lonlat_2d : vec_2d<T> {
-};
-
-template <typename T>
-struct alignas(2 * sizeof(T)) cartesian_2d : vec_2d<T> {
-};
+namespace cuspatial {
 
 /**
  * @brief A timestamp

cpp/include/cuspatial/utility/vec_2d.hpp

@@ -0,0 +1,84 @@
+#pragma once
+#include <cuspatial/types.hpp>
+
+namespace cuspatial {
+
+/**
+ * @brief A 2D vector
+ *
+ * Used in cuspatial for both Longitude/Latitude (LonLat) coordinate pairs and Cartesian (X/Y)
+ * coordinate pairs. For LonLat pairs, the `x` member represents Longitude, and `y` represents
+ * Latitude.
+ *
+ * @tparam T the base type for the coordinates
+ */
+template <typename T>
+struct alignas(2 * sizeof(T)) vec_2d {
+  using value_type = T;
+  value_type x;
+  value_type y;
+};
+
+template <typename T>
+struct alignas(2 * sizeof(T)) lonlat_2d : vec_2d<T> {
+};
+
+template <typename T>
+struct alignas(2 * sizeof(T)) cartesian_2d : vec_2d<T> {
+};
+
+/**
+ * @brief Element-wise add of two 2d vectors.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE operator+(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+  return vec_2d<T>{a.x + b.x, a.y + b.y};
+}
+
+/**
+ * @brief Element-wise subtract of two 2d vectors.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE operator-(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+  return vec_2d<T>{a.x - b.x, a.y - b.y};
+}
+
+/**
+ * @brief Scale a 2d vector by ratio @p r.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE operator*(vec_2d<T> vec, T const& r)
+{
+  return vec_2d<T>{vec.x * r, vec.y * r};
+}
+
+/**
+ * @brief Scale a 2d vector by ratio @p r.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE operator*(T const& r, vec_2d<T> vec)
+{
+  return vec * r;
+}
+
+/**
+ * @brief Compute dot product of two 2d vectors.
+ */
+template <typename T>
+T CUSPATIAL_HOST_DEVICE dot(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+  return a.x * b.x + a.y * b.y;
+}
+
+/**
+ * @brief Compute 2d determinant of a 2x2 matrix with column vectors @p a and @p b.
+ */
+template <typename T>
+T CUSPATIAL_HOST_DEVICE det(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+  return a.x * b.y - a.y * b.x;
+}
+
+}  // namespace cuspatial