[REVIEW] Directed Polygon Point Distance

CHANGELOG.md

@@ -14,6 +14,7 @@
 - PR #246 Hausdorff performance improvement
 - PR #253 Update conda upload versions for new supported CUDA/Python
 - PR #250 cartesian product iterator + more Hausdorff performance improvements.
+- PR #251 directed polygon distance
 
 ## Bug Fixes
 - PR #244 Restrict gdal version

cpp/CMakeLists.txt

@@ -311,6 +311,7 @@ add_library(cuspatial SHARED
             src/join/quadtree_poly_filtering.cu
             src/spatial/polygon_bounding_box.cu
             src/spatial/polyline_bounding_box.cu
+            src/spatial/polygon_distance.cu
             src/spatial/point_in_polygon.cu
             src/spatial_window/spatial_window.cu
             src/spatial/haversine.cu

cpp/benchmarks/CMakeLists.txt

@@ -75,4 +75,7 @@ link_directories("${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}" # CMAKE_CUDA_IMPLICIT
 set(HAUSDORFF_BENCH_SRC
     "${CMAKE_CURRENT_SOURCE_DIR}/hausdorff_benchmark.cpp")
 
-ConfigureBench(HAUSDORFF_BENCH "${HAUSDORFF_BENCH_SRC}")
+set(POLYGON_DISTANCE_BENCH_SRC
+    "${CMAKE_CURRENT_SOURCE_DIR}/polygon_distance_benchmark.cpp")
+
+ConfigureBench(POLYGON_DISTANCE_BENCH "${POLYGON_DISTANCE_BENCH_SRC}")

cpp/benchmarks/polygon_distance_benchmark.cpp

@@ -0,0 +1,66 @@
+/*
+ * Copyright (c) 2020, 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/polygon_distance.hpp>
+
+#include <benchmarks/fixture/benchmark_fixture.hpp>
+#include <benchmarks/synchronization/synchronization.hpp>
+
+#include <tests/utilities/column_wrapper.hpp>
+
+#include <thrust/iterator/constant_iterator.h>
+
+static void BM_polygon(benchmark::State& state)
+{
+  int32_t num_points           = state.range(1) - 1;
+  int32_t num_spaces_asked     = state.range(0) - 1;
+  int32_t num_spaces           = std::min(num_points, num_spaces_asked);
+  int32_t num_points_per_space = num_points / num_spaces;
+
+  auto counting_iter = thrust::counting_iterator<int32_t>();
+  auto zero_iter     = thrust::make_transform_iterator(counting_iter, [](auto idx) { return 0; });
+
+  auto space_offset_iter = thrust::make_transform_iterator(
+    counting_iter, [num_points_per_space](int32_t idx) { return idx * num_points_per_space; });
+
+  auto xs = cudf::test::fixed_width_column_wrapper<double>(zero_iter, zero_iter + num_points);
+  auto ys = cudf::test::fixed_width_column_wrapper<double>(zero_iter, zero_iter + num_points);
+
+  auto space_offsets = cudf::test::fixed_width_column_wrapper<int32_t>(
+    space_offset_iter, space_offset_iter + num_spaces);
+
+  for (auto _ : state) {
+    cuda_event_timer raii(state, true);
+    cuspatial::directed_polygon_distance(xs, ys, space_offsets);
+  }
+
+  state.SetItemsProcessed(state.iterations() * num_points * num_points);
+}
+
+class PolygonDistanceBenchmark : public cuspatial::benchmark {
+};
+
+#define DUMMY_BM_BENCHMARK_DEFINE(name)                                          \
+  BENCHMARK_DEFINE_F(PolygonDistanceBenchmark, name)(::benchmark::State & state) \
+  {                                                                              \
+    BM_polygon(state);                                                           \
+  }                                                                              \
+  BENCHMARK_REGISTER_F(PolygonDistanceBenchmark, name)                           \
+    ->Ranges({{1 << 10, 1 << 14}, {1 << 10, 1 << 15}})                           \
+    ->UseManualTime()                                                            \
+    ->Unit(benchmark::kMillisecond);
+
+DUMMY_BM_BENCHMARK_DEFINE(polygon);

cpp/include/cuspatial/hausdorff.hpp

@@ -71,18 +71,19 @@ namespace cuspatial {
  *
  * @param[in] xs: x component of points
  * @param[in] ys: y component of points
- * @param[in] points_per_space: number of points in each space
+ * @param[in] offsets: number of points in each space
+ * @param[in] mr: Device memory resource used to allocate the returned memory
  *
  * @returns Hausdorff distances for each pair of spaces
  *
- * @throw cudf::cuda_error if `points_per_space` contains negative values
+ * @throw cudf::cuda_error if `offsets` contains negative values
  *
  * @note Hausdorff distances are asymmetrical
  */
 std::unique_ptr<cudf::column> directed_hausdorff_distance(
   cudf::column_view const& xs,
   cudf::column_view const& ys,
-  cudf::column_view const& points_per_space,
+  cudf::column_view const& offsets,
   rmm::mr::device_memory_resource* mr = rmm::mr::get_default_resource());
 
 }  // namespace cuspatial

cpp/include/cuspatial/polygon_distance.hpp

@@ -0,0 +1,45 @@
+/*
+ * Copyright (c) 2020, 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 point_b_y 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>
+
+namespace cuspatial {
+
+/**
+ * @brief calculates minimum segment-to-point distance between all shapes
+ *
+ * Element `i + j*n` is the minimum distance from any segment in shape i to any point in shape j.
+ * The minimum of value of elements `[i + j*n]` and `j + i*n` is equal to the euclidian distance
+ * between points i and j.
+ *
+ *
+ * @param[in] xs: x component of points
+ * @param[in] ys: y component of points
+ * @param[in] offsets: number of points in each space
+ * @param[in] mr: Device memory resource used to allocate the returned memory
+ * @return std::unique_ptr<cudf::column>
+ */
+std::unique_ptr<cudf::column> directed_polygon_distance(
+  cudf::column_view const& xs,
+  cudf::column_view const& ys,
+  cudf::column_view const& offsets,
+  rmm::mr::device_memory_resource* mr = rmm::mr::get_default_resource());
+
+}  // namespace cuspatial

cpp/src/spatial/polygon_distance.cu

@@ -0,0 +1,169 @@
+/*
+ * Copyright (c) 2020, 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 point_b_y 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/cartesian_product_group_index_iterator.cuh>
+#include <cuspatial/error.hpp>
+
+#include <cudf/column/column_device_view.cuh>
+#include <cudf/column/column_factories.hpp>
+#include <cudf/types.hpp>
+#include <cudf/utilities/type_dispatcher.hpp>
+
+#include <thrust/iterator/discard_iterator.h>
+
+#include <limits>
+#include <memory>
+
+using size_type = cudf::size_type;
+
+namespace cuspatial {
+namespace detail {
+namespace {
+
+/**
+ * @brief Calculates segment-to-point distance
+ *
+ * Given a `cartesian_product_group_index` and two columns representing `x` and `y` coordinates,
+ * calculates the segment-to-point distance of a line segment in group `a` and a point in group `b`.
+ *
+ * If group `a` contains two or more points, calculate segment-to-point distance.
+ * If group `a` contains only a single point, calculate point-point distance.
+ */
+template <typename T>
+struct segment_point_distance_calculator {
+  cudf::column_device_view xs;
+  cudf::column_device_view ys;
+
+  T __device__ operator()(cartesian_product_group_index idx)
+  {
+    auto const a_idx_0 = idx.group_a.offset + (idx.element_a_idx);
+    auto const a_idx_1 = idx.group_a.offset + (idx.element_a_idx + 1) % idx.group_a.size;
+    auto const b_idx_0 = idx.group_b.offset + (idx.element_b_idx);
+
+    auto const origin_x = xs.element<T>(a_idx_0);
+    auto const origin_y = ys.element<T>(a_idx_0);
+    auto const edge_x   = xs.element<T>(a_idx_1) - origin_x;
+    auto const edge_y   = ys.element<T>(a_idx_1) - origin_y;
+    auto const point_x  = xs.element<T>(b_idx_0) - origin_x;
+    auto const point_y  = ys.element<T>(b_idx_0) - origin_y;
+
+    auto const magnitude = edge_x * edge_x + edge_y * edge_y;
+    auto const travel    = point_x * edge_x + point_y * edge_y;
+
+    // if point is projected within line segment bounds, use segment-to-point distance.
+    // if point is projected outside line segment bounds, use point-point distance.
+    if (0 < travel && travel < magnitude) {
+      return abs(point_y * edge_x - point_x * edge_y) * rhypot(edge_x, edge_y);
+    } else {
+      return hypot(point_x, point_y);
+    }
+  }
+};
+
+struct directed_polygon_distance_functor {
+  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("Non-floating point operation is not supported");
+  }
+
+  template <typename T>
+  std::enable_if_t<std::is_floating_point<T>::value, std::unique_ptr<cudf::column>> operator()(
+    cudf::column_view const& xs,
+    cudf::column_view const& ys,
+    cudf::column_view const& space_offsets,
+    rmm::mr::device_memory_resource* mr,
+    cudaStream_t stream)
+  {
+    size_type num_points  = xs.size();
+    size_type num_spaces  = space_offsets.size();
+    size_type num_results = num_spaces * num_spaces;
+
+    if (num_results == 0) {
+      return cudf::make_empty_column(cudf::data_type{cudf::type_to_id<T>()});
+    }
+
+    // ===== Make Separation and Key Iterators =====================================================
+
+    auto cartesian_iter = make_cartesian_product_group_index_iterator(
+      num_points, num_spaces, space_offsets.begin<cudf::size_type>());
+
+    auto cartesian_key_iter = thrust::make_transform_iterator(
+      cartesian_iter, [] __device__(cartesian_product_group_index idx) {
+        return thrust::make_pair(idx.group_a.idx, idx.group_b.idx);
+      });
+
+    auto d_xs = cudf::column_device_view::create(xs);
+    auto d_ys = cudf::column_device_view::create(ys);
+
+    auto separation_iter = thrust::make_transform_iterator(
+      cartesian_iter, segment_point_distance_calculator<T>{*d_xs, *d_ys});
+
+    // ===== Materialize ===========================================================================
+
+    auto result = cudf::make_fixed_width_column(cudf::data_type{cudf::type_to_id<T>()},
+                                                num_results,
+                                                cudf::mask_state::UNALLOCATED,
+                                                stream,
+                                                mr);
+
+    auto num_cartesian = num_points * num_points;
+
+    thrust::reduce_by_key(rmm::exec_policy(stream)->on(stream),
+                          cartesian_key_iter,
+                          cartesian_key_iter + num_cartesian,
+                          separation_iter,
+                          thrust::make_discard_iterator(),
+                          result->mutable_view().begin<T>(),
+                          thrust::equal_to<thrust::pair<int32_t, int32_t>>(),
+                          thrust::minimum<T>());
+
+    return result;
+  }
+};
+
+}  // namespace
+
+std::unique_ptr<cudf::column> directed_polygon_distance(cudf::column_view const& xs,
+                                                        cudf::column_view const& ys,
+                                                        cudf::column_view const& offsets,
+                                                        cudaStream_t stream,
+                                                        rmm::mr::device_memory_resource* mr)
+{
+  CUSPATIAL_EXPECTS(xs.type() == ys.type(), "Inputs `xs` and `ys` must have same type.");
+  CUSPATIAL_EXPECTS(xs.size() == ys.size(), "Inputs `xs` and `ys` must have same length.");
+
+  CUSPATIAL_EXPECTS(not xs.has_nulls() and not ys.has_nulls() and not offsets.has_nulls(),
+                    "Inputs must not have nulls.");
+
+  CUSPATIAL_EXPECTS(xs.size() >= offsets.size(), "At least one point is required for each space");
+
+  return cudf::type_dispatcher(
+    xs.type(), detail::directed_polygon_distance_functor(), xs, ys, offsets, mr, stream);
+}
+
+}  // namespace detail
+
+std::unique_ptr<cudf::column> directed_polygon_distance(cudf::column_view const& xs,
+                                                        cudf::column_view const& ys,
+                                                        cudf::column_view const& offsets,
+                                                        rmm::mr::device_memory_resource* mr)
+{
+  return detail::directed_polygon_distance(xs, ys, offsets, 0, mr);
+}
+
+}  // namespace cuspatial

cpp/tests/CMakeLists.txt

@@ -107,6 +107,10 @@ set(HAUSDORFF_TEST_SRC
    "${CMAKE_CURRENT_SOURCE_DIR}/spatial/hausdorff_test.cpp")
 ConfigureTest(HAUSDORFF_TEST "${HAUSDORFF_TEST_SRC}")
 
+set(POLYGON_DISTANCE_TEST_SRC
+   "${CMAKE_CURRENT_SOURCE_DIR}/spatial/polygon_distance_test.cpp")
+ConfigureTest(POLYGON_DISTANCE_TEST "${POLYGON_DISTANCE_TEST_SRC}")
+
 set(POINT_IN_POLYGON_TEST_SRC
     "${CMAKE_CURRENT_SOURCE_DIR}/spatial/point_in_polygon_test.cpp")
 ConfigureTest(POINT_IN_POLYGON_TEST "${POINT_IN_POLYGON_TEST_SRC}")