Refactor Hausdorff distance to header-only API

cpp/CMakeLists.txt

@@ -103,6 +103,7 @@ rapids_cpm_init()
 # find or add cuDF
 include(cmake/thirdparty/CUSPATIAL_GetCUDF.cmake)
 
+
 ###################################################################################################
 # - library targets -------------------------------------------------------------------------------
 

cpp/include/cuspatial/experimental/coordinate_transform.cuh

@@ -43,6 +43,10 @@ namespace cuspatial {
  * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible and mutable.
  * @tparam T the floating-point coordinate value type of input longitude/latitude coordinates.
  *
+ * @pre `lonlat_first` may equal `xy_first`, but the range `[lonlat_first, lonlat_last)`
+ * shall not otherwise overlap the range `[xy_first, xy_first + std::distance(lonlat_first,
+ * lonlat_last))`.
+ *
  * @return Output iterator to the element past the last x/y coordinate computed.
  *
  * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator

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

@@ -0,0 +1,178 @@
+/*
+ * 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/error.hpp>
+#include <cuspatial/utility/device_atomics.cuh>
+#include <cuspatial/utility/traits.hpp>
+#include <cuspatial/utility/vec_2d.hpp>
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/binary_search.h>
+
+#include <cuda/atomic>
+
+#include <type_traits>
+
+namespace cuspatial {
+
+namespace detail {
+
+template <typename T>
+constexpr auto magnitude_squared(T a, T b)
+{
+  return a * a + b * b;
+}
+
+/**
+ * @brief computes Hausdorff distance by equally dividing up work on a per-thread basis.
+ *
+ * Each thread is responsible for computing the distance from a single point in the input against
+ * all other points in the input. Because points in the input can originate from different spaces,
+ * each thread must know which spaces it is comparing. For the LHS argument, the point is always
+ * the same for any given thread and is determined once for that thread using a binary search of
+ * the provided space_offsets. Therefore if space 0 contains 10 points, the first 10 threads will
+ * know that the LHS space is 0. The 11th thread will know the LHS space is 1, and so on depending
+ * on the sizes/offsets of each space. Each thread then loops over each space, and uses an inner
+ * loop to loop over each point within that space, thereby knowing the RHS space and RHS point.
+ * the thread computes the minimum distance from it's LHS point to _any_ point in the RHS space, as
+ * this is the first step to computing Hausdorff distance. The second step of computing Hausdorff
+ * distance is to determine the maximum of these minimums, which is done by each thread writing
+ * it's minimum to the output using atomicMax. This is done once per thread per RHS space. Once
+ * all threads have run to completion, all "maximums of the minumum distances" (aka, directed
+ * Hausdorff distances) reside in the output.
+ *
+ * @tparam T type of coordinate, either float or double.
+ * @param num_points number of total points in the input (sum of points from all spaces)
+ * @param points x/y points to compute the distances between
+ * @param num_spaces number of spaces in the input
+ * @param space_offsets starting position of first point in each space
+ * @param results directed Hausdorff distances computed by kernel
+ */
+template <typename T,
+          typename Index,
+          typename PointsIter,
+          typename OffsetsIter,
+          typename OutputIter>
+__global__ void kernel_hausdorff(Index num_points,
+                                 PointsIter points,
+                                 Index num_spaces,
+                                 OffsetsIter space_offsets,
+                                 OutputIter results)
+{
+  using Point = typename std::iterator_traits<PointsIter>::value_type;
+
+  // determine the LHS point this thread is responsible for.
+  auto const thread_idx = blockIdx.x * blockDim.x + threadIdx.x;
+  Index const lhs_p_idx = thread_idx;
+
+  if (lhs_p_idx >= num_points) { return; }
+
+  // determine the LHS space this point belongs to.
+  Index const lhs_space_idx =
+    thrust::distance(
+      space_offsets,
+      thrust::upper_bound(thrust::seq, space_offsets, space_offsets + num_spaces, lhs_p_idx)) -
+    1;
+
+  // get the coordinates of this LHS point.
+  Point const lhs_p = points[lhs_p_idx];
+
+  // loop over each RHS space, as determined by space_offsets
+  for (uint32_t rhs_space_idx = 0; rhs_space_idx < num_spaces; rhs_space_idx++) {
+    // determine the begin/end offsets of points contained within this RHS space.
+    Index const rhs_p_idx_begin = space_offsets[rhs_space_idx];
+    Index const rhs_p_idx_end =
+      (rhs_space_idx + 1 == num_spaces) ? num_points : space_offsets[rhs_space_idx + 1];
+
+    // each space must contain at least one point, this initial value is just an identity value to
+    // simplify calculations. If a space contains <= 0 points, then this initial value will be
+    // written to the output, which can serve as a signal that the input is ill-formed.
+    auto min_distance_squared = std::numeric_limits<T>::max();
+
+    // loop over each point in the current RHS space
+    for (uint32_t rhs_p_idx = rhs_p_idx_begin; rhs_p_idx < rhs_p_idx_end; rhs_p_idx++) {
+      // get the x and y coordinate of this RHS point
+      Point const rhs_p = thrust::raw_reference_cast(points[rhs_p_idx]);
+
+      // get distance between the LHS and RHS point
+      auto const distance_squared = magnitude_squared(rhs_p.x - lhs_p.x, rhs_p.y - lhs_p.y);
+
+      // remember only smallest distance from this LHS point to any RHS point.
+      min_distance_squared = min(min_distance_squared, distance_squared);
+    }
+
+    // determine the output offset for this pair of spaces (LHS, RHS)
+    Index output_idx = lhs_space_idx * num_spaces + rhs_space_idx;
+
+    // use atomicMax to find the maximum of the minimum distance calculated for each space pair.
+    atomicMax(&thrust::raw_reference_cast(*(results + output_idx)),
+              static_cast<T>(std::sqrt(min_distance_squared)));
+  }
+}
+
+}  // namespace detail
+
+template <class PointIt, class OffsetIt, class OutputIt>
+OutputIt directed_hausdorff_distance(PointIt points_first,
+                                     PointIt points_last,
+                                     OffsetIt space_offsets_first,
+                                     OffsetIt space_offsets_last,
+                                     OutputIt distance_first,
+                                     rmm::cuda_stream_view stream)
+{
+  using Point   = typename std::iterator_traits<PointIt>::value_type;
+  using Index   = typename std::iterator_traits<OffsetIt>::value_type;
+  using T       = typename Point::value_type;
+  using OutputT = typename std::iterator_traits<OutputIt>::value_type;
+
+  static_assert(std::is_convertible_v<Point, cuspatial::vec_2d<T>>,
+                "Input points must be convertible to cuspatial::vec_2d");
+  static_assert(detail::is_floating_point<T, OutputT>(),
+                "Hausdorff supports only floating-point coordinates.");
+  static_assert(std::is_integral_v<Index>, "Indices must be integral");
+
+  auto const num_points = std::distance(points_first, points_last);
+  auto const num_spaces = std::distance(space_offsets_first, space_offsets_last);
+
+  CUSPATIAL_EXPECTS(num_points >= num_spaces, "At least one point is required for each space");
+  CUSPATIAL_EXPECTS(num_spaces < (1 << 15), "Total number of spaces must be less than 2^16");
+
+  auto const num_results = num_spaces * num_spaces;
+
+  if (num_results > 0) {
+    // Due to hausdorff kernel using `atomicMax` for output, the output must be initialized to <= 0
+    // here the output is being initialized to -1, which should always be overwritten. If -1 is
+    // found in the output, there is a bug where the output is not being written to in the hausdorff
+    // kernel.
+    thrust::fill_n(rmm::exec_policy(stream), distance_first, num_results, -1);
+
+    auto const threads_per_block = 64;
+    auto const num_tiles         = (num_points + threads_per_block - 1) / threads_per_block;
+
+    detail::kernel_hausdorff<T, decltype(num_points)>
+      <<<num_tiles, threads_per_block, 0, stream.value()>>>(
+        num_points, points_first, num_spaces, space_offsets_first, distance_first);
+
+    CUSPATIAL_CUDA_TRY(cudaGetLastError());
+  }
+
+  return distance_first + num_results;
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/hausdorff.cuh

@@ -0,0 +1,106 @@
+/*
+ * 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/types.hpp>
+#include <memory>
+
+#include <rmm/cuda_stream_view.hpp>
+
+namespace cuspatial {
+
+/**
+ * @ingroup distance
+ * @brief Computes Hausdorff distances for all pairs in a collection of spaces
+ *
+ * https://en.wikipedia.org/wiki/Hausdorff_distance
+ *
+ * Example in 1D (this function operates in 2D):
+ * ```
+ * spaces
+ * [0 2 5] [9] [3 7]
+ *
+ * spaces represented as points per space and concatenation of all points
+ * [0 2 5 9 3 7] [3 1 2]
+ *
+ * note: the following matrices are visually separated to highlight the relationship of a pair of
+ * points with the pair of spaces from which it is produced
+ *
+ * cartesian product of all
+ * points by pair of spaces     distance between points
+ * +----------+----+-------+    +---------+---+------+
+ * : 00 02 05 : 09 : 03 07 :    : 0  2  5 : 9 : 3  7 :
+ * : 20 22 25 : 29 : 23 27 :    : 2  0  3 : 7 : 1  5 :
+ * : 50 52 55 : 59 : 53 57 :    : 5  3  0 : 4 : 2  2 :
+ * +----------+----+-------+    +---------+---+------+
+ * : 90 92 95 : 99 : 93 97 :    : 9  7  4 : 0 : 6  2 :
+ * +----------+----+-------+    +---------+---+------+
+ * : 30 32 35 : 39 : 33 37 :    : 3  1  2 : 6 : 0  4 :
+ * : 70 72 75 : 79 : 73 77 :    : 7  5  2 : 2 : 4  0 :
+ * +----------+----+-------+    +---------+---+------+
+
+ * minimum distance from
+ * every point in one           Hausdorff distance is
+ * space to any point in        the maximum of the
+ * the other space              minimum distances
+ * +----------+----+-------+    +---------+---+------+
+ * :  0       :  9 :  3    :    : 0       : 9 : 3    :
+ * :     0    :  7 :  1    :    :         :   :      :
+ * :        0 :  4 :  2    :    :         :   :      :
+ * +----------+----+-------+    +---------+---+------+
+ * :        4 :  0 :     2 :    :       4 : 0 :    2 :
+ * +----------+----+-------+    +---------+---+------+
+ * :     1    :  6 :  0    :    :         : 6 : 0    :
+ * :        2 :  2 :     0 :    :       2 :   :      :
+ * +----------+----+-------+    +---------+---+------+
+ *
+ * returned as concatenation of columns
+ * [0 2 4 3 0 2 9 6 0]
+ * ```
+ *
+ * @param[in] points_first: xs: beginning of range of (x,y) points
+ * @param[in] points_lasts: xs: end of range of (x,y) points
+ * @param[in] space_offsets_first: beginning of range of indices to each space.
+ * @param[in] space_offsets_first: end of range of indices to each space. Last index is the last
+ * @param[in] distance_first: beginning of range of output Hausdorff distance for each pair of
+ * spaces
+ *
+ * @tparam PointIt Iterator to input points. Points must be of a type that is convertible to
+ * `cuspatial::vec_2d<T>`. Must meet the requirements of [LegacyRandomAccessIterator][LinkLRAI] and
+ * be device-accessible.
+ * @tparam OffsetIt Iterator to space offsets. Value type must be integral. Must meet the
+ * requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OutputIt Output iterator. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible and mutable.
+ *
+ * @pre All iterators must have the same underlying floating-point value type.
+ *
+ * @return Output iterator to the element past the last distance computed.
+ *
+ * @note Hausdorff distances are asymmetrical
+ */
+template <class PointIt, class OffsetIt, class OutputIt>
+OutputIt directed_hausdorff_distance(PointIt points_first,
+                                     PointIt points_last,
+                                     OffsetIt space_offsets_first,
+                                     OffsetIt space_offsets_last,
+                                     OutputIt distance_first,
+                                     rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/hausdorff.cuh>

cpp/include/cuspatial/experimental/haversine.cuh

@@ -56,12 +56,6 @@ namespace cuspatial {
  * `cuspatial::lonlat_2d<T>`.
  * @tparam T The underlying coordinate type. Must be a floating-point type.
  *
- * @pre `a_lonlat_first` may equal `distance_first`, but the range `[a_lonlat_first, a_lonlat_last)`
- * shall not overlap the range `[distance_first, distance_first + (a_lonlat_last - a_lonlat_last))
- * otherwise.
- * @pre `b_lonlat_first` may equal `distance_first`, but the range `[b_lonlat_first, b_lonlat_last)`
- * shall not overlap the range `[distance_first, distance_first + (b_lonlat_last - b_lonlat_last))
- * otherwise.
  * @pre All iterators must have the same `Location` type, with  the same underlying floating-point
  * coordinate type (e.g. `cuspatial::lonlat_2d<float>`).
  *

cpp/include/cuspatial/experimental/type_utils.hpp

@@ -14,6 +14,8 @@
  * limitations under the License.
  */
 
+#include <cuspatial/utility/vec_2d.hpp>
+
 #include <thrust/iterator/transform_iterator.h>
 #include <thrust/iterator/transform_output_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
@@ -105,6 +107,13 @@ auto make_vec_2d_iterator(FirstIter first, SecondIter second)
  * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
  * "LegacyRandomAccessIterator"
  */
+template <typename FirstIter, typename SecondIter>
+auto make_vec_2d_iterator(FirstIter first, SecondIter second)
+{
+  using T = typename std::iterator_traits<FirstIter>::value_type;
+  return make_vec_2d_iterator<vec_2d<T>>(first, second);
+}
+
 template <typename FirstIter, typename SecondIter>
 auto make_lonlat_iterator(FirstIter first, SecondIter second)
 {

cpp/include/cuspatial/utility/device_atomics.cuh

@@ -14,6 +14,8 @@
  * limitations under the License.
  */
 
+#pragma once
+
 namespace cuspatial {
 namespace detail {
 
@@ -30,7 +32,7 @@ namespace detail {
  * @param val The value to compare
  * @return The old value stored in `addr`.
  */
-__device__ double atomicMin(double* addr, double val)
+__device__ inline double atomicMin(double* addr, double val)
 {
   unsigned long long int* address_as_ll = reinterpret_cast<unsigned long long int*>(addr);
   unsigned long long int old            = __double_as_longlong(*addr);
@@ -60,7 +62,7 @@ __device__ double atomicMin(double* addr, double val)
  * @param val The value to compare
  * @return The old value stored in `addr`.
  */
-__device__ float atomicMin(float* addr, float val)
+__device__ inline float atomicMin(float* addr, float val)
 {
   unsigned int* address_as_ui = reinterpret_cast<unsigned int*>(addr);
   unsigned int old            = __float_as_uint(*addr);
@@ -90,7 +92,7 @@ __device__ float atomicMin(float* addr, float val)
  * @param val The value to compare
  * @return The old value stored in `addr`.
  */
-__device__ double atomicMax(double* addr, double val)
+__device__ inline double atomicMax(double* addr, double val)
 {
   unsigned long long int* address_as_ll = reinterpret_cast<unsigned long long int*>(addr);
   unsigned long long int old            = __double_as_longlong(*addr);
@@ -120,7 +122,7 @@ __device__ double atomicMax(double* addr, double val)
  * @param val The value to compare
  * @return The old value stored in `addr`.
  */
-__device__ float atomicMax(float* addr, float val)
+__device__ inline float atomicMax(float* addr, float val)
 {
   unsigned int* address_as_ui = reinterpret_cast<unsigned int*>(addr);
   unsigned int old            = __float_as_uint(*addr);

cpp/include/doxygen_groups.h

@@ -49,6 +49,7 @@
  *
  *          @file linestring_distance.hpp
  *          @file hausdorff.hpp
+ *          @file experimental/hausdorff.cuh
  *          @file haversine.hpp
  *          @file haversine.cuh
  *      @}