Header-only trajectory_bounding_boxes

cpp/include/cuspatial/detail/utility/device_atomics.cuh

@@ -29,7 +29,7 @@ namespace detail {
 template <typename T>
 const T __device__ min_(const T a, const T b)
 {
-  return min(a, b);
+  return ::min(a, b);
 }
 
 /**
@@ -40,7 +40,7 @@ const T __device__ min_(const T a, const T b)
 template <typename T>
 const T __device__ max_(const T a, const T b)
 {
-  return max(a, b);
+  return ::max(a, b);
 }
 
 /**

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

@@ -121,7 +121,7 @@ __global__ void kernel_hausdorff(
       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);
+      min_distance_squared = ::min(min_distance_squared, distance_squared);
     }
 
     // determine the output offset for this pair of spaces (LHS, RHS)

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

@@ -0,0 +1,74 @@
+/*
+ * 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/traits.hpp>
+#include <cuspatial/vec_2d.hpp>
+
+#include <rmm/cuda_stream_view.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <thrust/iterator/discard_iterator.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/reduce.h>
+#include <thrust/zip_function.h>
+
+namespace cuspatial {
+
+namespace detail {
+
+template <typename T>
+struct box_minmax {
+  using point_tuple = thrust::tuple<cuspatial::vec_2d<T>, cuspatial::vec_2d<T>>;
+  __host__ __device__ point_tuple operator()(point_tuple const& a, point_tuple const& b)
+  {
+    // structured binding doesn't seem to work with thrust::tuple
+    vec_2d<T> p1, p2, p3, p4;
+    thrust::tie(p1, p2) = a;
+    thrust::tie(p3, p4) = b;
+    return {box_min(box_min(p1, p2), p3), box_max(box_max(p1, p2), p4)};
+  }
+};
+
+}  // namespace detail
+
+template <typename IdInputIt, typename PointInputIt, typename BoundingBoxOutputIt>
+BoundingBoxOutputIt trajectory_bounding_boxes(IdInputIt ids_first,
+                                              IdInputIt ids_last,
+                                              PointInputIt points_first,
+                                              BoundingBoxOutputIt bounding_boxes_first,
+                                              rmm::cuda_stream_view stream)
+{
+  using T      = iterator_vec_base_type<PointInputIt>;
+  using IdType = iterator_value_type<IdInputIt>;
+
+  auto points_zipped_first = thrust::make_zip_iterator(points_first, points_first);
+
+  [[maybe_unused]] auto [_, bounding_boxes_last] =
+    thrust::reduce_by_key(rmm::exec_policy(stream),
+                          ids_first,
+                          ids_last,
+                          points_zipped_first,
+                          thrust::make_discard_iterator(),
+                          bounding_boxes_first,
+                          thrust::equal_to<IdType>(),
+                          detail::box_minmax<T>{});
+
+  return bounding_boxes_last;
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/iterator_factory.cuh

@@ -198,7 +198,7 @@ auto make_vec_2d_output_iterator(FirstIter first, SecondIter second)
 {
   using T         = typename std::iterator_traits<FirstIter>::value_type;
   auto zipped_out = thrust::make_zip_iterator(thrust::make_tuple(first, second));
-  return thrust::make_transform_output_iterator(zipped_out, detail::vec_2d_to_tuple<T>());
+  return thrust::transform_output_iterator(zipped_out, detail::vec_2d_to_tuple<T>());
 }
 
 /**

cpp/include/cuspatial/experimental/trajectory_bounding_boxes.cuh

@@ -0,0 +1,70 @@
+/*
+ * 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 <rmm/cuda_stream_view.hpp>
+
+namespace cuspatial {
+
+/**
+ * @addtogroup trajectory_api
+ * @{
+ */
+
+/**
+ * @brief Compute the spatial bounding boxes of trajectories.
+ *
+ * Computes a tight bounding box around all points within each trajectory (points with the same ID).
+ *
+ * @note Assumes Object IDs and points are presorted by ID. This can be done using
+ * cuspatial::derive_trajectories.
+ *
+ * @tparam IdInputIt Iterator over object IDs. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-readable.
+ * @tparam PointInputIt Iterator over points. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-readable.
+ * @tparam BoundingBoxOutputIt Iterator over output bounding boxes. Each element is a tuple of two
+ * points representing corners of the axis-aligned bounding box. The type of the points is the same
+ * as the `value_type` of PointInputIt. Must meet the requirements of
+ * [LegacyRandomAccessIterator][LinkLRAI] and be device-writeable.
+ *
+ * @param ids_first beginning of the range of input object ids
+ * @param ids_last end of the range of input object ids
+ * @param points_first beginning of the range of input point (x,y) coordinates
+ * @param bounding_boxes_first beginning of the range of output bounding boxes, one per trajectory
+ * @param stream the CUDA stream on which to perform computations.
+ *
+ * @return An iterator to the end of the range of output bounding boxes.
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <typename IdInputIt, typename PointInputIt, typename BoundingBoxOutputIt>
+BoundingBoxOutputIt trajectory_bounding_boxes(
+  IdInputIt ids_first,
+  IdInputIt ids_last,
+  PointInputIt points_first,
+  BoundingBoxOutputIt bounding_boxes_first,
+  rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+/**
+ * @} // end of doxygen group
+ */
+
+}  // namespace cuspatial
+
+#include "detail/trajectory_bounding_boxes.cuh"

cpp/include/cuspatial/vec_2d.hpp

@@ -18,6 +18,7 @@
 
 #include <cuspatial/cuda_utils.hpp>
 
+#include <algorithm>
 #include <ostream>
 
 namespace cuspatial {
@@ -115,6 +116,32 @@ T CUSPATIAL_HOST_DEVICE det(vec_2d<T> const& a, vec_2d<T> const& b)
   return a.x * b.y - a.y * b.x;
 }
 
+/**
+ * @brief Return a new vec_2d made up of the minimum x- and y-components of two input vec_2d values.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE box_min(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+#ifdef __CUDA_ARCH__
+  return vec_2d<T>{::min(a.x, b.x), ::min(a.y, b.y)};
+#else
+  return vec_2d<T>{std::min(a.x, b.x), std::min(a.y, b.y)};
+#endif
+}
+
+/**
+ * @brief Return a new vec_2d made up of the minimum x- and y-components of two input vec_2d values.
+ */
+template <typename T>
+vec_2d<T> CUSPATIAL_HOST_DEVICE box_max(vec_2d<T> const& a, vec_2d<T> const& b)
+{
+#ifdef __CUDA_ARCH__
+  return vec_2d<T>{::max(a.x, b.x), ::max(a.y, b.y)};
+#else
+  return vec_2d<T>{std::max(a.x, b.x), std::max(a.y, b.y)};
+#endif
+}
+
 /**
  * @} // end of doxygen group
  */

cpp/src/trajectory/derive_trajectories.cu

@@ -14,14 +14,14 @@
  * limitations under the License.
  */
 
-#include <cudf/utilities/type_dispatcher.hpp>
 #include <cuspatial/error.hpp>
 #include <cuspatial/experimental/derive_trajectories.cuh>
 #include <cuspatial/experimental/iterator_factory.cuh>
 
 #include <cudf/column/column_factories.hpp>
 #include <cudf/copying.hpp>
 #include <cudf/table/table.hpp>
+#include <cudf/utilities/type_dispatcher.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>

cpp/src/trajectory/trajectory_bounding_boxes.cu

@@ -15,7 +15,8 @@
  */
 
 #include <cuspatial/error.hpp>
-#include <cuspatial/trajectory.hpp>
+#include <cuspatial/experimental/iterator_factory.cuh>
+#include <cuspatial/experimental/trajectory_bounding_boxes.cuh>
 
 #include <cudf/column/column_device_view.cuh>
 #include <cudf/column/column_factories.hpp>
@@ -26,22 +27,19 @@
 #include <cudf/utilities/type_dispatcher.hpp>
 
 #include <rmm/cuda_stream_view.hpp>
-#include <rmm/exec_policy.hpp>
 
-#include <thrust/fill.h>
-#include <thrust/functional.h>
-#include <thrust/iterator/discard_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
-#include <thrust/reduce.h>
-#include <thrust/tuple.h>
+
+#include <type_traits>
+#include <vector>
 
 namespace cuspatial {
 
 namespace {
 
 struct dispatch_element {
-  template <typename Element>
-  std::enable_if_t<std::is_floating_point<Element>::value, std::unique_ptr<cudf::table>> operator()(
+  template <typename T>
+  std::enable_if_t<std::is_floating_point<T>::value, std::unique_ptr<cudf::table>> operator()(
     cudf::size_type num_trajectories,
     cudf::column_view const& object_id,
     cudf::column_view const& x,
@@ -50,7 +48,7 @@ struct dispatch_element {
     rmm::mr::device_memory_resource* mr)
   {
     // Construct output columns
-    auto type = cudf::data_type{cudf::type_to_id<Element>()};
+    auto type = cudf::data_type{cudf::type_to_id<T>()};
     std::vector<std::unique_ptr<cudf::column>> cols{};
     cols.reserve(4);
     // allocate bbox_x1 output column
@@ -66,40 +64,18 @@ struct dispatch_element {
     cols.push_back(
       cudf::make_numeric_column(type, num_trajectories, cudf::mask_state::UNALLOCATED, stream, mr));
 
-    auto points = thrust::make_zip_iterator(thrust::make_tuple(
-      x.begin<Element>(), y.begin<Element>(), x.begin<Element>(), y.begin<Element>()));
-
-    auto bboxes = thrust::make_zip_iterator(
-      thrust::make_tuple(cols.at(0)->mutable_view().begin<Element>(),  // bbox_x1
-                         cols.at(1)->mutable_view().begin<Element>(),  // bbox_y1
-                         cols.at(2)->mutable_view().begin<Element>(),  // bbox_x2
-                         cols.at(3)->mutable_view().begin<Element>())  // bbox_y2
-    );
-
-    thrust::fill(rmm::exec_policy(stream),
-                 bboxes,
-                 bboxes + num_trajectories,
-                 thrust::make_tuple(std::numeric_limits<Element>::max(),
-                                    std::numeric_limits<Element>::max(),
-                                    std::numeric_limits<Element>::min(),
-                                    std::numeric_limits<Element>::min()));
-
-    thrust::reduce_by_key(
-      rmm::exec_policy(stream),         // execution policy
-      object_id.begin<int32_t>(),       // keys_first
-      object_id.end<int32_t>(),         // keys_last
-      points,                           // values_first
-      thrust::make_discard_iterator(),  // keys_output
-      bboxes,                           // values_output
-      thrust::equal_to<int32_t>(),      // binary_pred
-      [] __device__(auto a, auto b) {   // binary_op
-        Element x1, y1, x2, y2, x3, y3, x4, y4;
-        thrust::tie(x1, y1, x2, y2) = a;
-        thrust::tie(x3, y3, x4, y4) = b;
-        return thrust::make_tuple(
-          min(min(x1, x2), x3), min(min(y1, y2), y3), max(max(x1, x2), x4), max(max(y1, y2), y4));
-      });
+    auto points_begin = cuspatial::make_vec_2d_iterator(x.begin<T>(), y.begin<T>());
+
+    auto bbox_mins  = cuspatial::make_vec_2d_output_iterator(cols.at(0)->mutable_view().begin<T>(),
+                                                            cols.at(1)->mutable_view().begin<T>());
+    auto bbox_maxes = cuspatial::make_vec_2d_output_iterator(cols.at(2)->mutable_view().begin<T>(),
+                                                             cols.at(3)->mutable_view().begin<T>());
 
+    trajectory_bounding_boxes(object_id.begin<cudf::size_type>(),
+                              object_id.end<cudf::size_type>(),
+                              points_begin,
+                              thrust::make_zip_iterator(bbox_mins, bbox_maxes),
+                              stream);
     // check for errors
     CUSPATIAL_CHECK_CUDA(stream.value());
 
@@ -137,7 +113,8 @@ std::unique_ptr<cudf::table> trajectory_bounding_boxes(cudf::size_type num_traje
 {
   CUSPATIAL_EXPECTS(object_id.size() == x.size() && x.size() == y.size(), "Data size mismatch");
   CUSPATIAL_EXPECTS(x.type().id() == y.type().id(), "Data type mismatch");
-  CUSPATIAL_EXPECTS(object_id.type().id() == cudf::type_id::INT32, "Invalid object_id type");
+  CUSPATIAL_EXPECTS(object_id.type().id() == cudf::type_to_id<cudf::size_type>(),
+                    "Invalid object_id type");
   CUSPATIAL_EXPECTS(!(x.has_nulls() || y.has_nulls() || object_id.has_nulls()),
                     "NULL support unimplemented");
 

cpp/tests/CMakeLists.txt

@@ -140,5 +140,8 @@ ConfigureTest(POINT_IN_POLYGON_TEST_EXP
 ConfigureTest(DERIVE_TRAJECTORIES_TEST_EXP
     experimental/trajectory/derive_trajectories_test.cu)
 
+ConfigureTest(TRAJECTORY_BOUNDING_BOXES_TEST_EXP
+    experimental/trajectory/trajectory_bounding_boxes_test.cu)
+
 ConfigureTest(POINT_QUADTREE_TEST_EXP
     experimental/indexing/point_quadtree_test.cu)