Header-only refactoring of points_in_spatial_window (#579)

Fixes #564 

A couple of notes:
1. For clarity and generality, this name for the header-only version of the function is now `copy_points_in_range()` rather than `points_in_spatial_window()`. This is the reflect that it is a range filter, and we may want to enable constraining the dimensions of the range in the future.
2. I have renamed the parameters for the range delimiters and clarified their documentation. This API is not limited to Cartesian coordinates and the documents make that clear by specifying that the input points and the range limits are specified using points in the same coordinate system. It also uses "quadrilateral" rather than "rectangle", since in lon/lat (spherical) coordinates the range is not a rectangle, but it is a quadrilateral.
3. Note that the change to the iterator-based API means that the user is responsible for allocating the output (as they are with STL and Thrust). Therefore, they need a way to know how much to allocate. In addition to the header-only `copy_points_in_range()`, this PR also adds a new function, `cuspatial::count_points_in_range()`. 

Previously, the cuDF-based API used `cudf::copy_if`, which does the equivalent of both of the above functions and returns a properly sized pair of columns. The new column-based implementation calls `count_points_in_range`, then allocates the result columns, and then calls `cuspatial::copy_points_in_range()`.

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

Approvers:
  - Michael Wang (https://github.com/isVoid)
  - H. Thomson Comer (https://github.com/thomcom)

URL: #579

cpp/benchmarks/CMakeLists.txt

@@ -85,8 +85,8 @@ ConfigureBench(HAUSDORFF_BENCH
 ConfigureNVBench(DISTANCES_BENCH
     pairwise_linestring_distance.cu)
 
+ConfigureNVBench(POINTS_IN_RANGE_BENCH
+    points_in_range.cu)
+
 ConfigureNVBench(POINT_IN_POLYGON_BENCH
     point_in_polygon.cu)
-
-ConfigureNVBench(SPATIAL_WINDOW_BENCH
-    spatial_window.cu)

cpp/benchmarks/spatial_window.cu → cpp/benchmarks/points_in_range.cu

@@ -19,8 +19,7 @@
 #include <benchmarks/utility/random.cuh>
 
 #include <cuspatial/detail/iterator.hpp>
-#include <cuspatial/experimental/type_utils.hpp>
-#include <cuspatial/spatial_window.hpp>
+#include <cuspatial/experimental/points_in_range.cuh>
 #include <cuspatial/vec_2d.hpp>
 
 #include <rmm/device_uvector.hpp>
@@ -38,72 +37,71 @@
 using namespace cuspatial;
 
 /**
- * @brief Helper to generate random points within a rectangular window
+ * @brief Helper to generate random points within a range
  *
  * @p begin and @p end must be iterators to device-accessible memory
  *
  * @tparam PointsIter The type of the iterator to the output points container
  * @tparam T The floating point type for the coordinates
- * @param begin The start of the range of points to generate
- * @param end The end of the range of points to generate
+ * @param begin The start of the sequence of points to generate
+ * @param end The end of the sequence of points to generate
  *
- * @param window_min the lower left window corner
- * @param window_max the upper right window corner
+ * @param range the lower left range corner
+ * @param range the upper right range corner
  *
  */
 template <class PointsIter, typename T>
-void generate_points(PointsIter begin, PointsIter end, vec_2d<T> window_min, vec_2d<T> window_max)
+void generate_points(PointsIter begin, PointsIter end, vec_2d<T> range_min, vec_2d<T> range_max)
 {
   auto engine_x = deterministic_engine(std::distance(begin, end));
   auto engine_y = deterministic_engine(2 * std::distance(begin, end));
 
-  auto x_dist = make_uniform_dist(window_min.x, window_max.x);
-  auto y_dist = make_uniform_dist(window_min.y, window_max.y);
+  auto x_dist = make_uniform_dist(range_min.x, range_max.x);
+  auto y_dist = make_uniform_dist(range_min.y, range_max.y);
 
-  auto x_gen = value_generator{window_min.x, window_max.x, engine_x, x_dist};
-  auto y_gen = value_generator{window_min.y, window_max.y, engine_y, y_dist};
+  auto x_gen = value_generator{range_min.x, range_max.x, engine_x, x_dist};
+  auto y_gen = value_generator{range_min.y, range_max.y, engine_y, y_dist};
 
   thrust::tabulate(rmm::exec_policy(), begin, end, [x_gen, y_gen] __device__(size_t n) mutable {
     return vec_2d<T>{x_gen(n), y_gen(n)};
   });
 }
 
 template <typename T>
-void points_in_spatial_window_benchmark(nvbench::state& state, nvbench::type_list<T>)
+void points_in_range_benchmark(nvbench::state& state, nvbench::type_list<T>)
 {
   // TODO: to be replaced by nvbench fixture once it's ready
   cuspatial::rmm_pool_raii rmm_pool;
 
   auto const num_points{state.get_int64("NumPoints")};
 
-  auto window_min = vec_2d<T>{-100, -100};
-  auto window_max = vec_2d<T>{100, 100};
+  auto range_min = vec_2d<T>{-100, -100};
+  auto range_max = vec_2d<T>{100, 100};
 
-  auto range_min = vec_2d<T>{-200, -200};
-  auto range_max = vec_2d<T>{200, 200};
+  auto generate_min = vec_2d<T>{-200, -200};
+  auto generate_max = vec_2d<T>{200, 200};
 
-  auto d_x = rmm::device_uvector<T>(num_points, rmm::cuda_stream_default);
-  auto d_y = rmm::device_uvector<T>(num_points, rmm::cuda_stream_default);
+  auto points = rmm::device_uvector<cuspatial::vec_2d<T>>(num_points, rmm::cuda_stream_default);
 
-  auto d_points =
-    cuspatial::make_zipped_vec_2d_output_iterator<cuspatial::vec_2d<T>>(d_x.begin(), d_y.begin());
-
-  generate_points(d_points, d_points + num_points, range_min, range_max);
-
-  auto xs = cudf::column(cudf::data_type{cudf::type_to_id<T>()}, num_points, d_x.release());
-  auto ys = cudf::column(cudf::data_type{cudf::type_to_id<T>()}, num_points, d_y.release());
+  generate_points(points.begin(), points.end(), generate_min, generate_max);
 
   CUSPATIAL_CUDA_TRY(cudaDeviceSynchronize());
 
   state.add_element_count(num_points);
 
   state.exec(nvbench::exec_tag::sync, [&](nvbench::launch& launch) {
-    auto points_in =
-      points_in_spatial_window(window_min.x, window_max.x, window_min.y, window_max.y, xs, ys);
+    auto stream = rmm::cuda_stream_view(launch.get_stream());
+    auto num_points_in =
+      cuspatial::count_points_in_range(range_min, range_max, points.begin(), points.end(), stream);
+
+    auto result_points = rmm::device_uvector<cuspatial::vec_2d<T>>(num_points_in, stream);
+
+    cuspatial::copy_points_in_range(
+      range_min, range_max, points.begin(), points.end(), result_points.begin(), stream);
   });
 }
 
 using floating_point_types = nvbench::type_list<float, double>;
-NVBENCH_BENCH_TYPES(points_in_spatial_window_benchmark, NVBENCH_TYPE_AXES(floating_point_types))
+NVBENCH_BENCH_TYPES(points_in_range_benchmark, NVBENCH_TYPE_AXES(floating_point_types))
   .set_type_axes_names({"CoordsType"})
   .add_int64_axis("NumPoints", {100'000, 1'000'000, 10'000'000, 100'000'000});

cpp/include/cuspatial/detail/utility/traits.hpp

@@ -23,7 +23,7 @@ namespace detail {
 
 /**
  * @internal
- * @brief returns true if all types are the same.
+ * @brief returns true if all types Ts... are the same as T.
  */
 template <typename T, typename... Ts>
 constexpr bool is_same()
@@ -33,7 +33,17 @@ constexpr bool is_same()
 
 /**
  * @internal
- * @brief returns true if all types are floating point types.
+ * @brief returns true if all types Ts... are convertible to U.
+ */
+template <typename U, typename... Ts>
+constexpr bool is_convertible_to()
+{
+  return std::conjunction_v<std::is_convertible<Ts, U>...>;
+}
+
+/**
+ * @internal
+ * @brief returns true if all types Ts... are floating point types.
  */
 template <typename... Ts>
 constexpr bool is_floating_point()

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

@@ -17,6 +17,7 @@
 #pragma once
 
 #include <cuspatial/constants.hpp>
+#include <cuspatial/vec_2d.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>

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

@@ -0,0 +1,92 @@
+/*
+ * 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 {
+namespace detail {
+
+// Functor to filter out points that are not inside the query range. Passed to thrust::copy_if
+template <typename T>
+struct range_filter {
+  range_filter(vec_2d<T> vertex_1, vec_2d<T> vertex_2)
+    : v1{std::min(vertex_1.x, vertex_2.x), std::min(vertex_1.y, vertex_2.y)},
+      v2{std::max(vertex_1.x, vertex_2.x), std::max(vertex_1.y, vertex_2.y)}
+  {
+  }
+
+  __device__ inline bool operator()(vec_2d<T> point)
+  {
+    return point.x > v1.x && point.x < v2.x && point.y > v1.y && point.y < v2.y;
+  }
+
+ protected:
+  vec_2d<T> v1;
+  vec_2d<T> v2;
+};
+
+}  // namespace detail
+
+template <class InputIt, class T>
+typename thrust::iterator_traits<InputIt>::difference_type count_points_in_range(
+  vec_2d<T> vertex_1,
+  vec_2d<T> vertex_2,
+  InputIt points_first,
+  InputIt points_last,
+  rmm::cuda_stream_view stream)
+{
+  using Point = typename std::iterator_traits<InputIt>::value_type;
+
+  static_assert(detail::is_convertible_to<cuspatial::vec_2d<T>, Point>(),
+                "Input points must be convertible to cuspatial::vec_2d");
+
+  return thrust::count_if(
+    rmm::exec_policy(stream), points_first, points_last, detail::range_filter{vertex_1, vertex_2});
+}
+
+template <class InputIt, class OutputIt, class T>
+OutputIt copy_points_in_range(vec_2d<T> vertex_1,
+                              vec_2d<T> vertex_2,
+                              InputIt points_first,
+                              InputIt points_last,
+                              OutputIt output_points_first,
+                              rmm::cuda_stream_view stream)
+{
+  using Point = typename std::iterator_traits<InputIt>::value_type;
+
+  static_assert(detail::is_convertible_to<cuspatial::vec_2d<T>, Point>(),
+                "Input points must be convertible to cuspatial::vec_2d");
+
+  static_assert(detail::is_same_floating_point<T, typename Point::value_type>(),
+                "Inputs and Range coordinates must have the same value type.");
+
+  return thrust::copy_if(rmm::exec_policy(stream),
+                         points_first,
+                         points_last,
+                         output_points_first,
+                         detail::range_filter{vertex_1, vertex_2});
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/points_in_range.cuh

@@ -0,0 +1,118 @@
+/*
+ * 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/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+
+#include <thrust/iterator/iterator_traits.h>
+
+namespace cuspatial {
+
+/**
+ * @brief Count of points (x,y) that fall within a query range.
+ *
+ * @ingroup spatial_relationship
+ *
+ * The query range is defined by a pair of opposite vertices within the coordinate system of the
+ * input points, `v1` and `v2`. A point (x, y) is in the range if `x` lies between `v1.x` and `v2.x`
+ * and `y` lies between `v1.y` and `v2.y`. A point is only counted if it is strictly within the
+ * interior of the query range. Points exactly on an edge or vertex of the range are not counted.
+ *
+ * The query range vertices and the input points are assumed to be defined in the same coordinate
+ * system.
+ *
+ * @param[in] vertex_1 Vertex of the query range quadrilateral
+ * @param[in] vertex_2 Vertex of the query range quadrilateral opposite `vertex_1`
+ * @param[in] points_first beginning of sequence of (x, y) coordinates of points to be queried
+ * @param[in] points_last end of sequence of (x, y) coordinates of points to be queried
+ * @param[in] stream The CUDA stream on which to perform computations
+ *
+ * @tparam InputIt Iterator to input points. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam T The underlying coordinate type. Must be a floating-point type.
+ *
+ * @pre All iterators must have the same `value_type`, with the same underlying floating-point
+ * coordinate type (e.g. `cuspatial::vec_2d<float>`).
+ *
+ * @returns The number of input points that fall within the specified query range.
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <class InputIt, class T>
+typename thrust::iterator_traits<InputIt>::difference_type count_points_in_range(
+  vec_2d<T> vertex_1,
+  vec_2d<T> vertex_2,
+  InputIt points_first,
+  InputIt points_last,
+  rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+/**
+ * @brief Copies points (x,y) that fall within a query range.
+ *
+ * @ingroup spatial_relationship
+ *
+ * The query range is defined by a pair of opposite vertices of a quadrilateral within the
+ * coordinate system of the input points, `v1` and `v2`. A point (x, y) is in the range if `x` lies
+ * between `v1.x` and `v2.x` and `y` lies between `v1.y` and `v2.y`. A point is only counted if it
+ * is strictly within the interior of the query range. Points exactly on an edge or vertex of the
+ * range are not copied.
+ *
+ * The query range vertices and the input points are assumed to be defined in the same coordinate
+ * system.
+ *
+ * `output_points_first` must be an iterator to storage of sufficient size for the points that will
+ * be copied. cuspatial::count_points_in_range may be used to determine the size required.
+ *
+ * @param[in] vertex_1 Vertex of the query range quadrilateral
+ * @param[in] vertex_2 Vertex of the query range quadrilateral opposite `vertex_1`
+ * @param[in] points_first beginning of sequence of (x, y) coordinates of points to be queried
+ * @param[in] points_last end of sequence of (x, y) coordinates of points to be queried
+ * @param[out] output_points_first beginning of output range of (x, y) coordinates within the
+ * query range
+ * @param[in] stream The CUDA stream on which to perform computations and allocate memory.
+ *
+ * @tparam InputIt Iterator to input points. 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.
+ * @tparam T The underlying coordinate type. Must be a floating-point type.
+ *
+ * @pre The range `[points_first, points_last)` may equal the range `[output_points_first,
+ * output_points_first + std::distance(points_first, points_last)), but the ranges may not
+ * partially overlap.
+ * @pre All iterators must have the same `value_type`, with the same underlying floating-point
+ * coordinate type (e.g. `cuspatial::vec_2d<float>`).
+ *
+ * @returns Output iterator to the element past the last output point.
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <class InputIt, class OutputIt, class T>
+OutputIt copy_points_in_range(vec_2d<T> vertex_1,
+                              vec_2d<T> vertex_2,
+                              InputIt points_first,
+                              InputIt points_last,
+                              OutputIt output_points_first,
+                              rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/points_in_range.cuh>

cpp/include/cuspatial/spatial_window.hpp

@@ -37,6 +37,9 @@ namespace cuspatial {
  * Swaps `window_min_x` and `window_max_x` if `window_min_x > window_max_x`.
  * Swaps `window_min_y` and `window_max_y` if `window_min_y > window_max_y`.
  *
+ * The window coordinates and the (x, y) points to be tested are assumed to be defined in the same
+ * coordinate system.
+ *
  * @param[in] window_min_x lower x-coordinate of the query window
  * @param[in] window_max_x upper x-coordinate of the query window
  * @param[in] window_min_y lower y-coordinate of the query window