Header-only Refactor of point_in_polygon

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

@@ -23,5 +23,44 @@ constexpr bool is_floating_point()
   return std::conjunction_v<std::is_floating_point<Ts>...>;
 }
 
+/**
+ * @internal
+ * @brief returns true if all types are floating point types.
+ */
+template <typename... Ts>
+constexpr bool is_integral()
+{
+  return std::conjunction_v<std::is_integral<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>...>;
+}
+
+/**
+ * @internal
+ * @brief Get the value type of @p Iterator type
+ *
+ * @tparam Iterator The iterator type to get from
+ */
+template <typename Iterator>
+using iterator_value_type = typename std::iterator_traits<Iterator>::value_type;
+
+/**
+ * @internal
+ * @brief Get the value type of the underlying vec_2d type from @p Iterator type
+ *
+ * @tparam Iterator The value type to get from, must point to a cuspatial::vec_2d
+ */
+template <typename Iterator>
+using iterator_vec_base_type = typename iterator_value_type<Iterator>::value_type;
+
 }  // namespace detail
 }  // namespace cuspatial

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

@@ -0,0 +1,174 @@
+/*
+ * 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 <thrust/memory.h>
+
+#include <iterator>
+#include <type_traits>
+
+namespace cuspatial {
+namespace detail {
+
+/**
+ * @brief Kernel to test if a point is inside all polygons.
+ *
+ * The algorithm is based on testing if the point is on one side of the segments in a polygon ring.
+ * Each point is tested against all segments in the polygon ring. If the point is on the the "same
+ * side" of a segment, it will flip the flag of `point_is_within`. Starting with `point_is_within`
+ * as `False`, a point is in the polygon if the flag is flipped odd number of times. Note that for a
+ * polygon ring with n vertices, the algorithm tests `n` segments (not `n-1`), including the segment
+ * between the last and first vertex.
+ */
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OutputIt>
+__global__ void point_in_polygon_kernel(Cart2dItA test_points_first,
+                                        int32_t const num_test_points,
+                                        OffsetIteratorA poly_offsets_first,
+                                        int32_t const num_polys,
+                                        OffsetIteratorB ring_offsets_first,
+                                        int32_t const num_rings,
+                                        Cart2dItB poly_points_first,
+                                        int32_t const num_poly_points,
+                                        OutputIt result)
+{
+  using T = iterator_vec_base_type<Cart2dItA>;
+
+  auto idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if (idx > num_test_points) { return; }
+
+  int32_t hit_mask = 0;
+
+  auto const test_point = thrust::raw_reference_cast(test_points_first[idx]);
+
+  // for each polygon
+  for (auto poly_idx = 0; poly_idx < num_polys; poly_idx++) {
+    auto poly_idx_next = poly_idx + 1;
+    auto poly_begin    = poly_offsets_first[poly_idx];
+    auto poly_end = (poly_idx_next < num_polys) ? poly_offsets_first[poly_idx_next] : num_rings;
+
+    bool point_is_within = false;
+
+    // for each ring
+    for (auto ring_idx = poly_begin; ring_idx < poly_end; ring_idx++) {
+      auto ring_idx_next = ring_idx + 1;
+      auto ring_begin    = ring_offsets_first[ring_idx];
+      auto ring_end =
+        (ring_idx_next < num_rings) ? ring_offsets_first[ring_idx_next] : num_poly_points;
+      auto ring_len = ring_end - ring_begin;
+
+      // for each line segment, including the segment between the last and first vertex
+      for (auto point_idx = 0; point_idx < ring_len; point_idx++) {
+        auto const a =
+          thrust::raw_reference_cast(poly_points_first[ring_begin + ((point_idx + 0) % ring_len)]);
+        auto const b =
+          thrust::raw_reference_cast(poly_points_first[ring_begin + ((point_idx + 1) % ring_len)]);
+
+        bool y_between_ay_by =
+          a.y <= test_point.y && test_point.y < b.y;  // is y in range [ay, by) when ay < by?
+        bool y_between_by_ay =
+          b.y <= test_point.y && test_point.y < a.y;  // is y in range [by, ay) when by < ay?
+        bool y_in_bounds = y_between_ay_by || y_between_by_ay;  // is y in range [by, ay]?
+        T run            = b.x - a.x;
+        T rise           = b.y - a.y;
+        T rise_to_point  = test_point.y - a.y;
+
+        if (y_in_bounds && test_point.x < (run / rise) * rise_to_point + a.x) {
+          point_is_within = not point_is_within;
+        }
+      }
+    }
+
+    hit_mask |= point_is_within << poly_idx;
+  }
+  result[idx] = hit_mask;
+}
+
+}  // namespace detail
+
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OutputIt>
+OutputIt point_in_polygon(Cart2dItA test_points_first,
+                          Cart2dItA test_points_last,
+                          OffsetIteratorA polygon_offsets_first,
+                          OffsetIteratorA polygon_offsets_last,
+                          OffsetIteratorB poly_ring_offsets_first,
+                          OffsetIteratorB poly_ring_offsets_last,
+                          Cart2dItB polygon_points_first,
+                          Cart2dItB polygon_points_last,
+                          OutputIt output,
+                          rmm::cuda_stream_view stream)
+{
+  using T = detail::iterator_vec_base_type<Cart2dItA>;
+
+  auto const num_test_points = std::distance(test_points_first, test_points_last);
+  auto const num_polys       = std::distance(polygon_offsets_first, polygon_offsets_last);
+  auto const num_rings       = std::distance(poly_ring_offsets_first, poly_ring_offsets_last);
+  auto const num_poly_points = std::distance(polygon_points_first, polygon_points_last);
+
+  static_assert(detail::is_same_floating_point<T, detail::iterator_vec_base_type<Cart2dItB>>(),
+                "Underlying type of Cart2dItA and Cart2dItB must be the same floating point type");
+  static_assert(detail::is_same<cartesian_2d<T>,
+                                detail::iterator_value_type<Cart2dItA>,
+                                detail::iterator_value_type<Cart2dItB>>(),
+                "Inputs must be cuspatial::cartesian_2d");
+
+  static_assert(detail::is_integral<detail::iterator_value_type<OffsetIteratorA>,
+                                    detail::iterator_value_type<OffsetIteratorB>>(),
+                "OffsetIterators must point to integral type.");
+
+  static_assert(std::is_same_v<detail::iterator_value_type<OutputIt>, int32_t>,
+                "OutputIt must point to 32 bit integer type.");
+
+  CUSPATIAL_EXPECTS(num_polys <= std::numeric_limits<int32_t>::digits,
+                    "Number of polygons cannot exceed 31");
+
+  CUSPATIAL_EXPECTS(num_rings >= num_polys, "Each polygon must have at least one ring");
+  CUSPATIAL_EXPECTS(num_poly_points >= num_polys * 4, "Each ring must have at least four vertices");
+
+  auto constexpr block_size = 256;
+  auto const num_blocks     = (num_test_points + block_size - 1) / block_size;
+
+  detail::point_in_polygon_kernel<<<num_blocks, block_size, 0, stream.value()>>>(
+    test_points_first,
+    num_test_points,
+    polygon_offsets_first,
+    num_polys,
+    poly_ring_offsets_first,
+    num_rings,
+    polygon_points_first,
+    num_poly_points,
+    output);
+  CUSPATIAL_CUDA_TRY(cudaGetLastError());
+
+  return output + num_test_points;
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/point_in_polygon.cuh

@@ -0,0 +1,114 @@
+/*
+ * 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 {
+
+/**
+ * @ingroup spatial_relationship
+ *
+ * @brief Tests whether the specified points are inside any of the specified polygons.
+ *
+ * Tests whether points are inside at most 31 polygons. Polygons are a collection of one or more
+ * rings. Rings are a collection of three or more vertices.
+ *
+ * Each input point will map to one `int32_t` element in the output. Each bit (except the sign bit)
+ * represents a hit or miss for each of the input polygons in least-significant-bit order. i.e.
+ * `output[3] & 0b0010` indicates a hit or miss for the 3rd point against the 2nd polygon.
+ *
+ *
+ * @tparam Cart2dItA iterator type for point array. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam Cart2dItB iterator type for point array. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorA iterator type for offset array. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OffsetIteratorB iterator type for offset array. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible.
+ * @tparam OutputIt iterator type for output array. Must meet
+ * the requirements of [LegacyRandomAccessIterator][LinkLRAI], be device-accessible, mutable and
+ * iterate on `int32_t` type.
+ *
+ * @param test_points_first begin of range of test points
+ * @param test_points_last end of range of test points
+ * @param polygon_offsets_first begin of range of indices to the first ring in each polygon
+ * @param polygon_offsets_last end of range of indices to the first ring in each polygon
+ * @param ring_offsets_first begin of range of indices to the first point in each ring
+ * @param ring_offsets_last end of range of indices to the first point in each ring
+ * @param polygon_points_first begin of range of polygon points
+ * @param polygon_points_last end of range of polygon points
+ * @param output begin iterator to the output buffer
+ * @param stream The CUDA stream to use for kernel launches.
+ * @return iterator to one past the last element in the output buffer
+ *
+ * @note Limit 31 polygons per call. Polygons may contain multiple rings.
+ * @note Direction of rings does not matter.
+ * @note This algorithm supports the ESRI shapefile format, but assumes all polygons are "clean" (as
+ * defined by the format), and does _not_ verify whether the input adheres to the shapefile format.
+ * @note Overlapping rings negate each other. This behavior is not limited to a single negation,
+ * allowing for "islands" within the same polygon.
+ * @note `poly_ring_offsets` must contain only the rings that make up the polygons indexed by
+ * `poly_offsets`. If there are rings in `poly_ring_offsets` that are not part of the polygons in
+ * `poly_offsets`, results are likely to be incorrect and behavior is undefined.
+ *
+ * ```
+ *   poly w/two rings         poly w/four rings
+ * +-----------+          +------------------------+
+ * :███████████:          :████████████████████████:
+ * :███████████:          :██+------------------+██:
+ * :██████+----:------+   :██:  +----+  +----+  :██:
+ * :██████:    :██████:   :██:  :████:  :████:  :██:
+ * +------;----+██████:   :██:  :----:  :----:  :██:
+ *        :███████████:   :██+------------------+██:
+ *        :███████████:   :████████████████████████:
+ *        +-----------+   +------------------------+
+ * ```
+ *
+ * @pre All point iterators must have the same `cartesian_2d` type, with  the same underlying
+ * floating-point coordinate type (e.g. `cuspatial::cartesian_2d<float>`).
+ * @pre All offset iterators must have the same `cartesian_2d` type, with  the same underlying
+ * integeral type.
+ * @pre Output iterator must be mutable and iterate on int32_t type.
+ *
+ * @throw cuspatial::logic_error if the number of polygons or rings exceeds 31.
+ * @throw cuspatial::logic_error polygon has less than 1 ring.
+ * @throw cuspatial::logic_error polygon has less than 4 vertices.
+ *
+ * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator
+ * "LegacyRandomAccessIterator"
+ */
+template <class Cart2dItA,
+          class Cart2dItB,
+          class OffsetIteratorA,
+          class OffsetIteratorB,
+          class OutputIt>
+OutputIt point_in_polygon(Cart2dItA test_points_first,
+                          Cart2dItA test_points_last,
+                          OffsetIteratorA polygon_offsets_first,
+                          OffsetIteratorA polygon_offsets_last,
+                          OffsetIteratorB poly_ring_offsets_first,
+                          OffsetIteratorB poly_ring_offsets_last,
+                          Cart2dItB polygon_points_first,
+                          Cart2dItB polygon_points_last,
+                          OutputIt output,
+                          rmm::cuda_stream_view stream = rmm::cuda_stream_default);
+
+}  // namespace cuspatial
+
+#include <cuspatial/experimental/detail/point_in_polygon.cuh>

cpp/include/doxygen_groups.h

@@ -58,6 +58,7 @@
  *          @brief APIs related to spatial relationship
  *
  *          @file point_in_polygon.hpp
+ *          @file experimental/point_in_polygon.cuh
  *          @file polygon_bounding_box.hpp
  *          @file polyline_bounding_box.hpp
  *          @file spatial_window.hpp