Add Header-Only point_polygon_distance, add non-owning polygon objects

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2022, NVIDIA CORPORATION.
+ * Copyright (c) 2022-2023, NVIDIA CORPORATION.
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -13,6 +13,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
+#pragma once
 
 #include <cuspatial/cuda_utils.hpp>
 

cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh → cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh

@@ -16,9 +16,12 @@
 
 #pragma once
 
+#include <cuspatial/experimental/geometry_collection/multipoint_ref.cuh>
 #include <cuspatial/traits.hpp>
+#include <cuspatial/vec_2d.hpp>
 
 #include <cuspatial/detail/utility/floating_point.cuh>
+#include <cuspatial/experimental/geometry/polygon_ref.cuh>
 
 namespace cuspatial {
 namespace detail {
@@ -30,44 +33,34 @@ namespace detail {
  * See "Crossings test" section of http://erich.realtimerendering.com/ptinpoly/
  * The improvement in addenda is also adopted to remove divisions in this kernel.
  *
+ * @tparam T type of coordinate
+ * @tparam PolygonRef polygon_ref type
+ * @param test_point point to test for point in polygon
+ * @param polygon polygon to test for point in polygon
+ * @return boolean to indicate if point is inside the polygon.
+ * `false` if point is on the edge of the polygon.
+ *
  * TODO: the ultimate goal of refactoring this as independent function is to remove
  * src/utility/point_in_polygon.cuh and its usage in quadtree_point_in_polygon.cu. It isn't
  * possible today without further work to refactor quadtree_point_in_polygon into header only
  * API.
  */
-template <class Cart2d,
-          class OffsetType,
-          class OffsetIterator,
-          class Cart2dIt,
-          class OffsetItDiffType = typename std::iterator_traits<OffsetIterator>::difference_type,
-          class Cart2dItDiffType = typename std::iterator_traits<Cart2dIt>::difference_type>
-__device__ inline bool is_point_in_polygon(Cart2d const& test_point,
-                                           OffsetType poly_begin,
-                                           OffsetType poly_end,
-                                           OffsetIterator ring_offsets_first,
-                                           OffsetItDiffType const& num_rings,
-                                           Cart2dIt poly_points_first,
-                                           Cart2dItDiffType const& num_poly_points)
+template <typename T, class PolygonRef>
+__device__ inline bool is_point_in_polygon(vec_2d<T> const& test_point, PolygonRef const& polygon)
 {
-  using T = iterator_vec_base_type<Cart2dIt>;
-
   bool point_is_within = false;
   bool is_colinear     = false;
-  // for each ring
-  for (auto ring_idx = poly_begin; ring_idx < poly_end; ring_idx++) {
-    int32_t ring_idx_next = ring_idx + 1;
-    int32_t ring_begin    = ring_offsets_first[ring_idx];
-    int32_t ring_end =
-      (ring_idx_next < num_rings) ? ring_offsets_first[ring_idx_next] : num_poly_points;
+  for (auto ring : polygon) {
+    auto last_segment = ring.segment(ring.num_segments() - 1);
 
-    Cart2d b     = poly_points_first[ring_end - 1];
+    auto b       = last_segment.v2;
     bool y0_flag = b.y > test_point.y;
     bool y1_flag;
-    // for each line segment, including the segment between the last and first vertex
-    for (auto point_idx = ring_begin; point_idx < ring_end; point_idx++) {
-      Cart2d const a = poly_points_first[point_idx];
-      T run          = b.x - a.x;
-      T rise         = b.y - a.y;
+    auto ring_points = multipoint_ref{ring.point_begin(), ring.point_end()};
+    for (vec_2d<T> a : ring_points) {
+      // for each line segment, including the segment between the last and first vertex
+      T run  = b.x - a.x;
+      T rise = b.y - a.y;
 
       // Points on the line segment are the same, so intersection is impossible.
       // This is possible because we allow closed or unclosed polygons.
@@ -100,5 +93,30 @@ __device__ inline bool is_point_in_polygon(Cart2d const& test_point,
 
   return point_is_within;
 }
+
+/**
+ * @brief Compatibility layer with non-OOP style input
+ */
+template <class Cart2d,
+          class OffsetType,
+          class OffsetIterator,
+          class Cart2dIt,
+          class OffsetItDiffType = typename std::iterator_traits<OffsetIterator>::difference_type,
+          class Cart2dItDiffType = typename std::iterator_traits<Cart2dIt>::difference_type>
+__device__ inline bool is_point_in_polygon(Cart2d const& test_point,
+                                           OffsetType poly_begin,
+                                           OffsetType poly_end,
+                                           OffsetIterator ring_offsets_first,
+                                           OffsetItDiffType const& num_rings,
+                                           Cart2dIt poly_points_first,
+                                           Cart2dItDiffType const& num_poly_points)
+{
+  auto polygon = polygon_ref{thrust::next(ring_offsets_first, poly_begin),
+                             thrust::next(ring_offsets_first, poly_end + 1),
+                             poly_points_first,
+                             thrust::next(poly_points_first, num_poly_points)};
+  return is_point_in_polygon(test_point, polygon);
+}
+
 }  // namespace detail
 }  // namespace cuspatial

cpp/include/cuspatial/experimental/detail/geometry/linestring_ref.cuh

@@ -54,7 +54,7 @@ template <typename VecIterator>
 CUSPATIAL_HOST_DEVICE auto linestring_ref<VecIterator>::num_segments() const
 {
   // The number of segment equals the number of points minus 1. And the number of points
-  // is thrust::distance(_point_begin, _point_end) - 1.
+  // is thrust::distance(_point_begin, _point_end).
   return thrust::distance(_point_begin, _point_end) - 1;
 }
 
@@ -70,6 +70,18 @@ CUSPATIAL_HOST_DEVICE auto linestring_ref<VecIterator>::segment_end() const
   return segment_begin() + num_segments();
 }
 
+template <typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto linestring_ref<VecIterator>::point_begin() const
+{
+  return _point_begin;
+}
+
+template <typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto linestring_ref<VecIterator>::point_end() const
+{
+  return _point_end;
+}
+
 template <typename VecIterator>
 template <typename IndexType>
 CUSPATIAL_HOST_DEVICE auto linestring_ref<VecIterator>::segment(IndexType i) const

cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh

@@ -0,0 +1,78 @@
+/*
+ * Copyright (c) 2023, 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/cuda_utils.hpp>
+#include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/experimental/geometry/segment.cuh>
+#include <cuspatial/experimental/geometry_collection/multilinestring_ref.cuh>
+#include <cuspatial/traits.hpp>
+
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/tuple.h>
+
+namespace cuspatial {
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE polygon_ref<RingIterator, VecIterator>::polygon_ref(RingIterator ring_begin,
+                                                                          RingIterator ring_end,
+                                                                          VecIterator point_begin,
+                                                                          VecIterator point_end)
+  : _ring_begin(ring_begin), _ring_end(ring_end), _point_begin(point_begin), _point_end(point_end)
+{
+  using T = iterator_vec_base_type<VecIterator>;
+  static_assert(is_same<vec_2d<T>, iterator_value_type<VecIterator>>(), "must be vec2d type");
+}
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::num_rings() const
+{
+  return thrust::distance(_ring_begin, _ring_end) - 1;
+}
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::ring_begin() const
+{
+  return detail::make_counting_transform_iterator(0,
+                                                  to_linestring_functor{_ring_begin, _point_begin});
+}
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::ring_end() const
+{
+  return ring_begin() + size();
+}
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::point_begin() const
+{
+  return _point_begin;
+}
+
+template <typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::point_end() const
+{
+  return _point_end;
+}
+
+template <typename RingIterator, typename VecIterator>
+template <typename IndexType>
+CUSPATIAL_HOST_DEVICE auto polygon_ref<RingIterator, VecIterator>::ring(IndexType i) const
+{
+  return *(ring_begin() + i);
+}
+
+}  // namespace cuspatial

cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh

@@ -1,6 +1,22 @@
+/*
+ * Copyright (c) 2023, 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/cuda_utils.hpp"
+#include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/detail/iterator.hpp>
 #include <cuspatial/experimental/geometry/linestring_ref.cuh>
 

cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh

@@ -0,0 +1,117 @@
+#pragma once
+
+#include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/detail/iterator.hpp>
+#include <cuspatial/experimental/geometry/linestring_ref.cuh>
+#include <cuspatial/experimental/geometry/polygon_ref.cuh>
+
+#include <thrust/iterator/transform_iterator.h>
+
+namespace cuspatial {
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+struct to_polygon_functor {
+  using difference_type = typename thrust::iterator_difference<PartIterator>::type;
+  PartIterator part_begin;
+  RingIterator ring_begin;
+  VecIterator point_begin;
+  VecIterator point_end;
+
+  CUSPATIAL_HOST_DEVICE
+  to_polygon_functor(PartIterator part_begin,
+                     RingIterator ring_begin,
+                     VecIterator point_begin,
+                     VecIterator point_end)
+    : part_begin(part_begin), ring_begin(ring_begin), point_begin(point_begin), point_end(point_end)
+  {
+  }
+
+  CUSPATIAL_HOST_DEVICE auto operator()(difference_type i)
+  {
+    return polygon_ref{ring_begin + part_begin[i],
+                       thrust::next(ring_begin + part_begin[i + 1]),
+                       point_begin,
+                       point_end};
+  }
+};
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+class multipolygon_ref;
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE multipolygon_ref<PartIterator, RingIterator, VecIterator>::multipolygon_ref(
+  PartIterator part_begin,
+  PartIterator part_end,
+  RingIterator ring_begin,
+  RingIterator ring_end,
+  VecIterator point_begin,
+  VecIterator point_end)
+  : _part_begin(part_begin),
+    _part_end(part_end),
+    _ring_begin(ring_begin),
+    _ring_end(ring_end),
+    _point_begin(point_begin),
+    _point_end(point_end)
+{
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::num_polygons()
+  const
+{
+  return thrust::distance(_part_begin, _part_end) - 1;
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::part_begin()
+  const
+{
+  return detail::make_counting_transform_iterator(
+    0, to_polygon_functor{_part_begin, _ring_begin, _point_begin, _point_end});
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::part_end()
+  const
+{
+  return part_begin() + num_polygons();
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::ring_begin()
+  const
+{
+  return detail::make_counting_transform_iterator(0,
+                                                  to_linestring_functor{_part_begin, _point_begin});
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::ring_end()
+  const
+{
+  return part_begin() + num_polygons();
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::point_begin()
+  const
+{
+  return _point_begin;
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::point_end()
+  const
+{
+  return _point_end;
+}
+
+template <typename PartIterator, typename RingIterator, typename VecIterator>
+template <typename IndexType>
+CUSPATIAL_HOST_DEVICE auto multipolygon_ref<PartIterator, RingIterator, VecIterator>::operator[](
+  IndexType i) const
+{
+  return *(part_begin() + i);
+}
+
+}  // namespace cuspatial