Header-only refactoring of points_in_spatial_window

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,13 +33,34 @@ 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()
 {
   return std::conjunction_v<std::is_floating_point<Ts>...>;
 }
 
+/**
+ * @internal
+ * @brief returns true if T and all types Ts... are the same floating point type.
+ */
+template <typename T, typename... Ts>
+constexpr bool is_same_floating_point()
+{
+  return std::conjunction_v<std::is_same<T, Ts>...> and
+         std::conjunction_v<std::is_floating_point<Ts>...>;
+}
+
 }  // namespace detail
 }  // namespace cuspatial

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/spatial_window.cuh

@@ -0,0 +1,96 @@
+/*
+ * 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 window
+// This is passed to thrust::copy_if
+template <typename T>
+struct spatial_window_filter {
+  spatial_window_filter(vec_2d<T> window_min, vec_2d<T> window_max)
+    : min{std::min(window_min.x, window_max.x),   // support mirrored rectangles
+          std::min(window_min.y, window_max.y)},  // where specified min > max
+      max{std::max(window_min.x, window_max.x), std::max(window_min.y, window_max.y)}
+  {
+  }
+
+  __device__ inline bool operator()(vec_2d<T> point)
+  {
+    return point.x > min.x && point.x < max.x && point.y > min.y && point.y < max.y;
+  }
+
+ protected:
+  vec_2d<T> min;
+  vec_2d<T> max;
+};
+
+}  // namespace detail
+
+template <class InputIt, class T>
+typename thrust::iterator_traits<InputIt>::difference_type count_points_in_spatial_window(
+  vec_2d<T> window_min,
+  vec_2d<T> window_max,
+  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::spatial_window_filter{window_min, window_max});
+}
+
+template <class InputIt, class OutputIt, class T>
+OutputIt points_in_spatial_window(vec_2d<T> window_min,
+                                  vec_2d<T> window_max,
+                                  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 window coordinates must have the same value type.");
+
+  return thrust::copy_if(rmm::exec_policy(stream),
+                         points_first,
+                         points_last,
+                         output_points_first,
+                         detail::spatial_window_filter{window_min, window_max});
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/spatial_window.cuh

@@ -0,0 +1,104 @@
+/*
+ * 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 rectangular query window.
+ *
+ * @ingroup spatial_relationship
+ *
+ * A point (x, y) is in the window if `x > window_min_x && x < window_max_x && y > window_min_y &&
+ * y < window_max_y`.
+ *
+ * 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`.
+ *
+ * @param[in] window_min lower-left (x, y) coordinate of the query window
+ * @param[in] window_max upper-right (x, y) coordinate of the query window
+ * @param[in] points_first beginning of range of (x, y) coordinates of points to be queried
+ * @param[in] points_last end of range 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 window.
+ */
+template <class InputIt, class T>
+typename thrust::iterator_traits<InputIt>::difference_type count_points_in_spatial_window(
+  vec_2d<T> window_min,
+  vec_2d<T> window_max,
+  InputIt points_first,
+  InputIt points_last,
+  rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+/**
+ * @brief Find all points (x,y) that fall within a rectangular query window.
+ *
+ * @ingroup spatial_relationship
+ *
+ * A point (x, y) is in the window if `x > window_min_x && x < window_max_x && y > window_min_y &&
+ * y < window_max_y`.
+ *
+ * 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`.
+ *
+ * @param[in] window_min lower-left (x, y) coordinate of the query window
+ * @param[in] window_max upper-right (x, y) coordinate of the query window
+ * @param[in] points_first beginning of range of (x, y) coordinates of points to be queried
+ * @param[in] points_last end of range 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 window
+ * @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.
+ */
+template <class InputIt, class OutputIt, class T>
+OutputIt points_in_spatial_window(vec_2d<T> window_min,
+                                  vec_2d<T> window_max,
+                                  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/spatial_window.cuh>

cpp/src/spatial_window/spatial_window.cu

@@ -15,11 +15,14 @@
  */
 
 #include <cuspatial/error.hpp>
+#include <cuspatial/experimental/spatial_window.cuh>
+#include <cuspatial/experimental/type_utils.hpp>
 
-#include <cudf/column/column_device_view.cuh>
+#include <cudf/column/column.hpp>
 #include <cudf/column/column_view.hpp>
 #include <cudf/copying.hpp>
-#include <cudf/detail/copy_if.cuh>
+#include <cudf/detail/copy.hpp>
+#include <cudf/table/table.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
 
@@ -28,41 +31,6 @@
 
 namespace {
 
-// Functor to filter out points that are not inside the query window
-// This is passed to cudf::detail::copy_if
-template <typename T>
-struct spatial_window_filter {
-  spatial_window_filter(T window_min_x,
-                        T window_max_x,
-                        T window_min_y,
-                        T window_max_y,
-                        cudf::column_device_view const& x,
-                        cudf::column_device_view const& y)
-    : min_x{std::min(window_min_x, window_max_x)},  // support mirrored rectangles
-      max_x{std::max(window_min_x, window_max_x)},  // where specified min > max
-      min_y{std::min(window_min_y, window_max_y)},
-      max_y{std::max(window_min_y, window_max_y)},
-      points_x{x},
-      points_y{y}
-  {
-  }
-
-  __device__ inline bool operator()(cudf::size_type i)
-  {
-    auto x = points_x.element<T>(i);
-    auto y = points_y.element<T>(i);
-    return x > min_x && x < max_x && y > min_y && y < max_y;
-  }
-
- protected:
-  T min_x;
-  T max_x;
-  T min_y;
-  T max_y;
-  cudf::column_device_view points_x;
-  cudf::column_device_view points_y;
-};
-
 // Type-dispatch functor that creates the spatial window filter of the correct type.
 // Only floating point types are supported.
 struct spatial_window_dispatch {
@@ -76,17 +44,32 @@ struct spatial_window_dispatch {
                                           rmm::cuda_stream_view stream,
                                           rmm::mr::device_memory_resource* mr)
   {
-    auto device_x = cudf::column_device_view::create(x, stream);
-    auto device_y = cudf::column_device_view::create(y, stream);
-    return cudf::detail::copy_if(cudf::table_view{{x, y}},
-                                 spatial_window_filter<T>{static_cast<T>(window_min_x),
-                                                          static_cast<T>(window_max_x),
-                                                          static_cast<T>(window_min_y),
-                                                          static_cast<T>(window_max_y),
-                                                          *device_x,
-                                                          *device_y},
-                                 stream,
-                                 mr);
+    auto points_begin = cuspatial::make_cartesian_2d_iterator(x.begin<T>(), y.begin<T>());
+
+    auto window_min =
+      cuspatial::vec_2d<T>{static_cast<T>(window_min_x), static_cast<T>(window_min_y)};
+    auto window_max =
+      cuspatial::vec_2d<T>{static_cast<T>(window_max_x), static_cast<T>(window_max_y)};
+
+    auto output_size = cuspatial::count_points_in_spatial_window(
+      window_min, window_max, points_begin, points_begin + x.size(), stream);
+
+    std::vector<std::unique_ptr<cudf::column>> cols{};
+    cols.reserve(2);
+    auto mask_policy = cudf::mask_allocation_policy::NEVER;
+    cols.push_back(cudf::detail::allocate_like(x, output_size, mask_policy, stream, mr));
+    cols.push_back(cudf::detail::allocate_like(y, output_size, mask_policy, stream, mr));
+
+    auto& output_x = cols[0];
+    auto& output_y = cols[1];
+
+    auto output_zip = cuspatial::make_zipped_cartesian_2d_output_iterator(
+      output_x->mutable_view().begin<T>(), output_y->mutable_view().begin<T>());
+
+    cuspatial::points_in_spatial_window(
+      window_min, window_max, points_begin, points_begin + x.size(), output_zip, stream);
+
+    return std::make_unique<cudf::table>(std::move(cols));
   }
 
   template <typename T,

cpp/tests/CMakeLists.txt

@@ -108,9 +108,11 @@ ConfigureTest(HAVERSINE_TEST_EXP
 ConfigureTest(HAUSDORFF_TEST_EXP
     experimental/spatial/hausdorff_test.cu)
 
-
 ConfigureTest(LINESTRING_DISTANCE_TEST_EXP
     experimental/spatial/linestring_distance_test.cu)
 
 ConfigureTest(COORDINATE_TRANSFORM_TEST_EXP
     experimental/spatial/coordinate_transform_test.cu)
+
+ConfigureTest(SPATIAL_WINDOW_POINT_TEST_EXP
+    experimental/spatial/spatial_window_test.cu)