From 965b80acaa7ce73363086f7b9350b12702957042 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 28 Feb 2023 10:18:10 -0800 Subject: [PATCH 01/34] initial --- .../detail/geometry/polygon_ref.cuh | 66 ++++++ .../geometry_collection/multipolygon_ref.cuh | 110 +++++++++ .../detail/point_polygon_distance.cuh | 87 +++++++ .../detail/ranges/multilinestring_range.cuh | 2 +- .../detail/ranges/multipolygon_range.cuh | 223 ++++++++++++++++++ .../experimental/geometry/polygon_ref.cuh | 64 +++++ .../geometry_collection/multipolygon_ref.cuh | 79 +++++++ .../experimental/point_polygon_distance.cuh | 41 ++++ .../ranges/multipolygon_range.cuh | 136 +++++++++++ 9 files changed, 807 insertions(+), 1 deletion(-) create mode 100644 cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh create mode 100644 cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh create mode 100644 cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh create mode 100644 cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh create mode 100644 cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh create mode 100644 cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh create mode 100644 cpp/include/cuspatial/experimental/point_polygon_distance.cuh create mode 100644 cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh diff --git a/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh new file mode 100644 index 000000000..8b156ea6a --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh @@ -0,0 +1,66 @@ +/* + * 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 "linestring_ref.cuh" + +#include +#include +#include +#include + +#include +#include + +namespace cuspatial { + +template +CUSPATIAL_HOST_DEVICE polygon_ref::polygon_ref(RingIterator ring_begin, + RingIterator ring_end, + VecIterator point_begin, + VecIterator point_end) + : _ring_begin(ring_begin), _ring_end(ring_end), _point_begin(begin), _point_end(end) +{ + using T = iterator_vec_base_type; + static_assert(is_same, iterator_value_type>(), "must be vec2d type"); +} + +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::num_rings() const +{ + return thrust::distance(_point_begin, _point_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::ring_begin() const +{ + return detail::make_counting_transform_iterator(0, to_linestring_ref{_point_begin}); +} + +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::ring_end() const +{ + return ring_begin() + num_rings(); +} + +template +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::ring(IndexType i) const +{ + return *(ring_begin() + i); +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh new file mode 100644 index 000000000..fe4e60cde --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh @@ -0,0 +1,110 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace cuspatial { + +template +struct to_polygon_functor { + using difference_type = typename thrust::iterator_difference::type; + PartIterator part_begin; + RingIterator ring_begin; + VecIterator point_begin; + + CUSPATIAL_HOST_DEVICE + to_polygon_functor(PartIterator part_begin, RingIterator ring_begin, VecIterator point_begin) + : part_begin(part_begin), ring_begin(ring_begin), point_begin(point_begin) + { + } + + CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) + { + return polygon_ref{point_begin + part_begin[i], point_begin + part_begin[i + 1]}; + } +}; + +template +class multipolygon_ref; + +template +CUSPATIAL_HOST_DEVICE multipolygon_ref::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 +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::num_polygons() + const +{ + return thrust::distance(_part_begin, _part_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::part_begin() + const +{ + return detail::make_counting_transform_iterator( + 0, to_polygon_functor{_part_begin, _ring_begin, _point_begin}); +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::part_end() + const +{ + return part_begin() + num_polygons(); +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::ring_begin() + const +{ + return detail::make_counting_transform_iterator(0, + to_linestring_functor{_part_begin, _point_begin}); +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::ring_end() + const +{ + return part_begin() + num_polygons(); +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::point_begin() + const +{ + return _point_begin; +} + +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::point_end() + const +{ + return _point_end; +} + +template +template +CUSPATIAL_HOST_DEVICE auto multipolygon_ref::operator[]( + IndexType i) const +{ + return *(part_begin() + i); +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh new file mode 100644 index 000000000..b51c252aa --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -0,0 +1,87 @@ +/* + * 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 +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace cuspatial { +namespace detail { + +/** + * @brief Kernel to compute the distance between pairs of point and polygon. + */ +template +void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoints, + MultiPolygonRange multipolygons, + OutputIterator distances) +{ + using T = typename MultiPointRange::element_t; + + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multipolygons.num_points(); + idx += gridDim.x * blockDim.x) { + auto geometry_idx = multipolygons.geometry_idx_from_point_idx(idx); + } +} + +} // namespace detail +template +OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, + MultiPolygonRange multipoiygons, + OutputIt distances_first, + rmm::cuda_stream_view stream) +{ + using T = typename MultiPointRange::element_t; + + static_assert(is_same_floating_point(), + "Inputs must have same floating point value type."); + + static_assert( + is_same, typename MultiPointRange::point_t, typename MultiPolygonRange::point_t>(), + "Inputs must be cuspatial::vec_2d"); + + CUSPATIAL_EXPECTS(multipoints.size() == multipolygons.size(), + "Must have the same number of input rows."); + + auto [threads_per_block, n_blocks] = grid_id(multipolygons.num_points()); + + pairwise_point_polygon_distance_kernel<<>>( + multipoints, multipolygons, distances_first); + + CUSPATIAL_CHECK_CUDA(stream.value()); +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh index dd373be3f..851269313 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh @@ -16,7 +16,6 @@ #pragma once -#include #include #include #include @@ -29,6 +28,7 @@ #include #include +#include namespace cuspatial { diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh new file mode 100644 index 000000000..7fe9a9387 --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -0,0 +1,223 @@ +/* + * 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 +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +namespace cuspatial { + +using namespace detail; + +template +struct to_multipolygon_functor { + using difference_type = typename thrust::iterator_difference::type; + GeometryIterator _geometry_begin; + PartIterator _part_begin; + RingIterator _ring_begin; + VecIterator _point_begin; + VecIterator _point_end; + + CUSPATIAL_HOST_DEVICE + to_multipolygon_functor(GeometryIterator geometry_begin, + PartIterator part_begin, + RingIterator ring_begin, + VecIterator point_begin, + VecIterator point_end) + : _geometry_begin(geometry_begin), + _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 multipolygon_ref{_part_begin + _geometry_begin[i], + thrust::next(_part_begin + _geometry_begin[i + 1]), + _point_begin, + _point_end}; + } +}; + +template +class multipolygon_range; + +template +multipolygon_range::multipolygon_range( + GeometryIterator geometry_begin, + GeometryIterator geometry_end, + PartIterator part_begin, + PartIterator part_end, + VecIterator point_begin, + VecIterator point_end) + : _geometry_begin(geometry_begin), + _geometry_end(geometry_end), + _part_begin(part_begin), + _part_end(part_end), + _point_begin(point_begin), + _point_end(point_end) +{ +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::num_multipolygons() +{ + return thrust::distance(_geometry_begin, _geometry_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::num_polygons() +{ + return thrust::distance(_part_begin, _part_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::num_rings() +{ + return thrust::distance(_ring_begin, _ring_end) - 1; +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::num_points() +{ + return thrust::distance(_point_begin, _point_end); +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::multipolygon_begin() +{ + return detail::make_counting_transform_iterator( + 0, + to_multipolygon_functor{_geometry_begin, _part_begin, _ring_begin, _point_begin, _point_end}); +} + +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::multipolygon_end() +{ + return multipolygon_begin() + num_multipolygons(); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::ring_idx_from_point_idx( + IndexType point_idx) +{ + return thrust::distance(_ring_begin, + thrust::prev(thrust::upper_bound(_ring_begin, _ring_end, point_idx))); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::part_idx_from_ring_idx( + IndexType ring_idx) +{ + return thrust::distance(_part_begin, + thrust::prev(thrust::upper_bound(_part_begin, _part_begin, ring_idx))); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::geometry_idx_from_part_idx( + IndexType part_idx) +{ + return thrust::distance( + _geometry_begin, thrust::prev(thrust::upper_bound(_geometry_begin, _geometry_end, part_idx))); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::geometry_idx_from_point_idx( + IndexType point_idx) +{ + return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx_from_part_idx(point_idx))); +} + +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::operator[]( + IndexType multipolygon_idx) +{ + return multipolygon_begin()[multipolygon_idx]; +} + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh b/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh new file mode 100644 index 000000000..42c3ce1a9 --- /dev/null +++ b/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh @@ -0,0 +1,64 @@ +/* + * 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 +#include + +namespace cuspatial { + +/** + * @brief Represent a reference to a polygon stored in a structure of arrays. + * + * @tparam VecIterator type of iterator to the underlying point array. + */ +template +class polygon_ref { + public: + CUSPATIAL_HOST_DEVICE polygon_ref(RingIterator ring_begin, + RingIterator ring_end, + VecIterator point_begin, + VecIterator point_end); + + /// Return the number of rings in the polygon + CUSPATIAL_HOST_DEVICE auto num_rings() const; + + /// Return iterator to the first ring of the polygon + CUSPATIAL_HOST_DEVICE auto ring_begin() const; + /// Return iterator to one past the last ring + CUSPATIAL_HOST_DEVICE auto ring_end() const; + + /// Return iterator to the first ring of the polygon + CUSPATIAL_HOST_DEVICE auto begin() const { return ring_begin(); } + /// Return iterator to one past the last ring + CUSPATIAL_HOST_DEVICE auto end() const { return ring_end(); } + + /// Return an enumerated range to the rings. + CUSPATIAL_HOST_DEVICE auto enumerate() { return detail::enumerate_range{begin(), end()}; } + + /// Return the `ring_idx`th ring in the polygon. + template + CUSPATIAL_HOST_DEVICE auto ring(IndexType ring_idx) const; + + protected: + RingIterator _ring_begin; + RingIterator _ring_end; + VecIterator _point_begin; + VecIterator _point_end; +}; + +} // namespace cuspatial +#include diff --git a/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh new file mode 100644 index 000000000..a6aab7cc4 --- /dev/null +++ b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh @@ -0,0 +1,79 @@ +/* + * 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 +#include + +namespace cuspatial { + +/** + * @brief Represent a reference to a multipolygon stored in a structure of arrays. + * + * @tparam PartIterator type of iterator to the part offset array. + * @tparam VecIterator type of iterator to the underlying point array. + */ +template +class multipolygon_ref { + public: + CUSPATIAL_HOST_DEVICE multipolygon_ref(PartIterator part_begin, + PartIterator part_end, + VecIterator point_begin, + VecIterator point_end); + /// Return the number of polygons in the multipolygon. + CUSPATIAL_HOST_DEVICE auto num_polygons() const; + /// Return the number of polygons in the multipolygon. + CUSPATIAL_HOST_DEVICE auto size() const { return num_polygons(); } + + /// Return iterator to the first polygon. + CUSPATIAL_HOST_DEVICE auto part_begin() const; + /// Return iterator to one past the last polygon. + CUSPATIAL_HOST_DEVICE auto part_end() const; + + /// Return iterator to the first polygon. + CUSPATIAL_HOST_DEVICE auto ring_begin() const; + /// Return iterator to one past the last polygon. + CUSPATIAL_HOST_DEVICE auto ring_begin() const; + + /// Return iterator to the first point of the multipolygon. + CUSPATIAL_HOST_DEVICE auto point_begin() const; + /// Return iterator to one past the last point of the multipolygon. + CUSPATIAL_HOST_DEVICE auto point_end() const; + + /// Return iterator to the first polygon of the multipolygon. + CUSPATIAL_HOST_DEVICE auto begin() const { return part_begin(); } + /// Return iterator to one past the last polygon of the multipolygon. + CUSPATIAL_HOST_DEVICE auto end() const { return part_end(); } + + /// Return an enumerated range to the polygons. + CUSPATIAL_HOST_DEVICE auto enumerate() const { return detail::enumerate_range{begin(), end()}; } + + /// Return `polygon_idx`th polygon in the multipolygon. + template + CUSPATIAL_HOST_DEVICE auto operator[](IndexType polygon_idx) const; + + protected: + PartIterator _part_begin; + PartIterator _part_end; + RingIterator _ring_begin; + RingIterator _ring_end; + VecIterator _point_begin; + VecIterator _point_end; +}; + +} // namespace cuspatial + +#include diff --git a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh new file mode 100644 index 000000000..c5394e902 --- /dev/null +++ b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh @@ -0,0 +1,41 @@ +/* + * 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 + +namespace cuspatial { + +/** + * @ingroup distance + * @copybrief cuspatial::pairwise_point_polygon_distance + * + * @tparam MultiPointRange An instance of template type `MultiPointRange` + * @tparam MultiPolygonRangeB An instance of template type `MultiPolygonRange` + * + * @param multipoints Range of multipoints in each distance pair. + * @param multipolygons Range of multipolygons in each distance pair. + * @return Iterator past the last distance computed + */ +template +OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, + MultiPolygonRangeB multipoiygons, + OutputIt distances_first, + rmm::cuda_stream_view stream = rmm::cuda_stream_default); +} // namespace cuspatial + +#include diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh new file mode 100644 index 000000000..331302348 --- /dev/null +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -0,0 +1,136 @@ +/* + * 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 + +#include +#include +#include +#include + +namespace cuspatial { + +/** + * @brief Non-owning range-based interface to multipolygon data + * @ingroup ranges + * + * Provides a range-based interface to contiguous storage of multipolygon data, to make it easier + * to access and iterate over multipolygons, polygons, rings and points. + * + * Conforms to GeoArrow's specification of multipolygon: + * https://github.com/geopandas/geo-arrow-spec/blob/main/format.md + * + * @tparam GeometryIterator iterator type for the geometry offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam PartIterator iterator type for the part offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam RingIterator iterator type for the ring offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam VecIterator iterator type for the point array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * + * @note Though this object is host/device compatible, + * The underlying iterator must be device accessible if used in device kernel. + * + * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator + * "LegacyRandomAccessIterator" + */ +template +class multipolygon_range { + public: + using geometry_it_t = GeometryIterator; + using part_it_t = PartIterator; + using ring_it_t = RingIterator; + using point_it_t = VecIterator; + using point_t = iterator_value_type; + using element_t = iterator_vec_base_type; + + multipolygon_range(GeometryIterator geometry_begin, + GeometryIterator geometry_end, + PartIterator part_begin, + PartIterator part_end, + RingIterator ring_begin, + RingIterator ring_end, + VecIterator points_begin, + VecIterator points_end); + + /// Return the number of multipolygons in the array. + CUSPATIAL_HOST_DEVICE auto size() { return num_multipolygons(); } + + /// Return the number of multipolygons in the array. + CUSPATIAL_HOST_DEVICE auto num_multipolygons(); + + /// Return the total number of polygons in the array. + CUSPATIAL_HOST_DEVICE auto num_polygons(); + + /// Return the total number of rings in the array. + CUSPATIAL_HOST_DEVICE auto num_rings(); + + /// Return the total number of points in the array. + CUSPATIAL_HOST_DEVICE auto num_points(); + + /// Return the iterator to the first multipolygon in the range. + CUSPATIAL_HOST_DEVICE auto multipolygon_begin(); + + /// Return the iterator to the one past the last multipolygon in the range. + CUSPATIAL_HOST_DEVICE auto multipolygon_end(); + + /// Return the iterator to the first multipolygon in the range. + CUSPATIAL_HOST_DEVICE auto begin() { return multipolygon_begin(); } + + /// Return the iterator to the one past the last multipolygon in the range. + CUSPATIAL_HOST_DEVICE auto end() { return multipolygon_end(); } + + /// Given the index of a point, return the ring index where the point locates. + template + CUSPATIAL_HOST_DEVICE auto ring_idx_from_point_idx(IndexType point_idx); + + /// Given the index of a ring, return the part (polygon) index + /// where the polygon locates. + template + CUSPATIAL_HOST_DEVICE auto part_idx_from_ring_idx(IndexType ring_idx); + + /// Given the index of a part (polygon), return the geometry (multipolygon) index + /// where the polygon locates. + template + CUSPATIAL_HOST_DEVICE auto geometry_idx_from_part_idx(IndexType part_idx); + + /// Given the index of a point, return the geometry (multipolygon) index where the + /// point locates. + template + CUSPATIAL_HOST_DEVICE auto geometry_idx_from_point_idx(IndexType point_idx); + + /// Returns the `multipolygon_idx`th multipolygon in the range. + template + CUSPATIAL_HOST_DEVICE auto operator[](IndexType multipolygon_idx); + + protected: + GeometryIterator _geometry_begin; + GeometryIterator _geometry_end; + PartIterator _part_begin; + PartIterator _part_end; + VecIterator _point_begin; + VecIterator _point_end; +}; + +} // namespace cuspatial + +#include From 02f61feca933bca3de3f545d7331b1f2d8fb678b Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:01:25 -0800 Subject: [PATCH 02/34] add pragma once for floating_point.cuh --- cpp/include/cuspatial/detail/utility/floating_point.cuh | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/include/cuspatial/detail/utility/floating_point.cuh b/cpp/include/cuspatial/detail/utility/floating_point.cuh index 744c7cb88..7a44aeb73 100644 --- a/cpp/include/cuspatial/detail/utility/floating_point.cuh +++ b/cpp/include/cuspatial/detail/utility/floating_point.cuh @@ -13,6 +13,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#pragma once #include From 7bd797fc23864c9cc7a83e0e4ac743a0b64ebdb4 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:02:20 -0800 Subject: [PATCH 03/34] add polygon_ref structure --- .../detail/geometry/polygon_ref.cuh | 26 ++++++++++++++----- .../experimental/geometry/polygon_ref.cuh | 12 ++++++--- 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh index 8b156ea6a..02a917cd4 100644 --- a/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh @@ -15,11 +15,10 @@ */ #pragma once -#include "linestring_ref.cuh" - #include #include #include +#include #include #include @@ -32,28 +31,41 @@ CUSPATIAL_HOST_DEVICE polygon_ref::polygon_ref(RingIt RingIterator ring_end, VecIterator point_begin, VecIterator point_end) - : _ring_begin(ring_begin), _ring_end(ring_end), _point_begin(begin), _point_end(end) + : _ring_begin(ring_begin), _ring_end(ring_end), _point_begin(point_begin), _point_end(point_end) { - using T = iterator_vec_base_type; + using T = iterator_vec_base_type; static_assert(is_same, iterator_value_type>(), "must be vec2d type"); } template CUSPATIAL_HOST_DEVICE auto polygon_ref::num_rings() const { - return thrust::distance(_point_begin, _point_end) - 1; + return thrust::distance(_ring_begin, _ring_end) - 1; } template CUSPATIAL_HOST_DEVICE auto polygon_ref::ring_begin() const { - return detail::make_counting_transform_iterator(0, to_linestring_ref{_point_begin}); + return detail::make_counting_transform_iterator(0, + to_linestring_functor{_ring_begin, _point_begin}); } template CUSPATIAL_HOST_DEVICE auto polygon_ref::ring_end() const { - return ring_begin() + num_rings(); + return ring_begin() + size(); +} + +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::point_begin() const +{ + return _point_begin; +} + +template +CUSPATIAL_HOST_DEVICE auto polygon_ref::point_end() const +{ + return _point_end; } template diff --git a/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh b/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh index 42c3ce1a9..5905eab2c 100644 --- a/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry/polygon_ref.cuh @@ -16,7 +16,6 @@ #pragma once #include -#include namespace cuspatial { @@ -36,19 +35,24 @@ class polygon_ref { /// Return the number of rings in the polygon CUSPATIAL_HOST_DEVICE auto num_rings() const; + /// Return the number of rings in the polygon + CUSPATIAL_HOST_DEVICE auto size() const { return num_rings(); } + /// Return iterator to the first ring of the polygon CUSPATIAL_HOST_DEVICE auto ring_begin() const; /// Return iterator to one past the last ring CUSPATIAL_HOST_DEVICE auto ring_end() const; + /// Return iterator to the first point of the polygon + CUSPATIAL_HOST_DEVICE auto point_begin() const; + /// Return iterator to one past the last point + CUSPATIAL_HOST_DEVICE auto point_end() const; + /// Return iterator to the first ring of the polygon CUSPATIAL_HOST_DEVICE auto begin() const { return ring_begin(); } /// Return iterator to one past the last ring CUSPATIAL_HOST_DEVICE auto end() const { return ring_end(); } - /// Return an enumerated range to the rings. - CUSPATIAL_HOST_DEVICE auto enumerate() { return detail::enumerate_range{begin(), end()}; } - /// Return the `ring_idx`th ring in the polygon. template CUSPATIAL_HOST_DEVICE auto ring(IndexType ring_idx) const; From 0968c15b72d3632700b284bacd4dbc9f2481bfa5 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:04:23 -0800 Subject: [PATCH 04/34] add multipolygon_ref class --- .../geometry_collection/multipolygon_ref.cuh | 15 +++++++++++---- .../geometry_collection/multipolygon_ref.cuh | 10 ++++++---- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh index fe4e60cde..940c6ef83 100644 --- a/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/detail/geometry_collection/multipolygon_ref.cuh @@ -15,16 +15,23 @@ struct to_polygon_functor { 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) - : part_begin(part_begin), ring_begin(ring_begin), point_begin(point_begin) + 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{point_begin + part_begin[i], point_begin + part_begin[i + 1]}; + return polygon_ref{ring_begin + part_begin[i], + thrust::next(ring_begin + part_begin[i + 1]), + point_begin, + point_end}; } }; @@ -60,7 +67,7 @@ CUSPATIAL_HOST_DEVICE auto multipolygon_ref diff --git a/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh index a6aab7cc4..64da05797 100644 --- a/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh @@ -26,11 +26,13 @@ namespace cuspatial { * @tparam PartIterator type of iterator to the part offset array. * @tparam VecIterator type of iterator to the underlying point array. */ -template +template class multipolygon_ref { public: CUSPATIAL_HOST_DEVICE multipolygon_ref(PartIterator part_begin, PartIterator part_end, + RingIterator ring_begin, + RingIterator ring_end, VecIterator point_begin, VecIterator point_end); /// Return the number of polygons in the multipolygon. @@ -43,10 +45,10 @@ class multipolygon_ref { /// Return iterator to one past the last polygon. CUSPATIAL_HOST_DEVICE auto part_end() const; - /// Return iterator to the first polygon. - CUSPATIAL_HOST_DEVICE auto ring_begin() const; - /// Return iterator to one past the last polygon. + /// Return iterator to the first ring. CUSPATIAL_HOST_DEVICE auto ring_begin() const; + /// Return iterator to one past the last ring. + CUSPATIAL_HOST_DEVICE auto ring_end() const; /// Return iterator to the first point of the multipolygon. CUSPATIAL_HOST_DEVICE auto point_begin() const; From a659eab1df056666af3fec2ac382c1221f6c5d91 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:05:06 -0800 Subject: [PATCH 05/34] update multipolygon_range class --- .../detail/ranges/multipolygon_range.cuh | 104 ++++++++++++++---- .../ranges/multipolygon_range.cuh | 51 ++++++--- 2 files changed, 114 insertions(+), 41 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index 7fe9a9387..b44420532 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -43,6 +44,7 @@ struct to_multipolygon_functor { GeometryIterator _geometry_begin; PartIterator _part_begin; RingIterator _ring_begin; + RingIterator _ring_end; VecIterator _point_begin; VecIterator _point_end; @@ -50,11 +52,13 @@ struct to_multipolygon_functor { to_multipolygon_functor(GeometryIterator geometry_begin, PartIterator part_begin, RingIterator ring_begin, + RingIterator ring_end, VecIterator point_begin, VecIterator point_end) : _geometry_begin(geometry_begin), _part_begin(part_begin), _ring_begin(ring_begin), + _ring_end(ring_end), _point_begin(point_begin), _point_end(point_end) { @@ -62,8 +66,22 @@ struct to_multipolygon_functor { CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) { + auto new_part_begin = _part_begin + _geometry_begin[i]; + auto new_part_end = thrust::next(_part_begin, _geometry_begin[i + 1] + 1); + + printf( + "In to_multipolygon_functor: %d %d %d\n\t new_part_begin_val: %d, " + "new_part_end_dist_from_begin: %d\n", + static_cast(i), + static_cast(_geometry_begin[i]), + static_cast(_geometry_begin[i + 1]), + static_cast(*new_part_begin), + static_cast(thrust::distance(new_part_end, new_part_begin))); + return multipolygon_ref{_part_begin + _geometry_begin[i], - thrust::next(_part_begin + _geometry_begin[i + 1]), + thrust::next(_part_begin, _geometry_begin[i + 1] + 1), + _ring_begin, + _ring_end, _point_begin, _point_end}; } @@ -79,17 +97,21 @@ template -multipolygon_range::multipolygon_range( +multipolygon_range::multipolygon_range( GeometryIterator geometry_begin, GeometryIterator geometry_end, PartIterator part_begin, PartIterator part_end, + RingIterator ring_begin, + RingIterator ring_end, VecIterator point_begin, VecIterator point_end) : _geometry_begin(geometry_begin), _geometry_end(geometry_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) { @@ -100,7 +122,7 @@ template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_multipolygons() +multipolygon_range::num_multipolygons() { return thrust::distance(_geometry_begin, _geometry_end) - 1; } @@ -110,7 +132,7 @@ template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_polygons() +multipolygon_range::num_polygons() { return thrust::distance(_part_begin, _part_end) - 1; } @@ -120,7 +142,7 @@ template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_rings() +multipolygon_range::num_rings() { return thrust::distance(_ring_begin, _ring_end) - 1; } @@ -130,7 +152,7 @@ template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_points() +multipolygon_range::num_points() { return thrust::distance(_point_begin, _point_end); } @@ -140,11 +162,12 @@ template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::multipolygon_begin() +multipolygon_range::multipolygon_begin() { return detail::make_counting_transform_iterator( 0, - to_multipolygon_functor{_geometry_begin, _part_begin, _ring_begin, _point_begin, _point_end}); + to_multipolygon_functor{ + _geometry_begin, _part_begin, _ring_begin, _ring_end, _point_begin, _point_end}); } template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::multipolygon_end() +multipolygon_range::multipolygon_end() { return multipolygon_begin() + num_multipolygons(); } @@ -163,11 +186,11 @@ template template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::ring_idx_from_point_idx( - IndexType point_idx) +multipolygon_range:: + ring_idx_from_point_idx(IndexType point_idx) { - return thrust::distance(_ring_begin, - thrust::prev(thrust::upper_bound(_ring_begin, _ring_end, point_idx))); + return thrust::distance( + _ring_begin, thrust::prev(thrust::upper_bound(thrust::seq, _ring_begin, _ring_end, point_idx))); } template template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::part_idx_from_ring_idx( - IndexType ring_idx) +multipolygon_range:: + part_idx_from_ring_idx(IndexType ring_idx) { - return thrust::distance(_part_begin, - thrust::prev(thrust::upper_bound(_part_begin, _part_begin, ring_idx))); + return thrust::distance( + _part_begin, + thrust::prev(thrust::upper_bound(thrust::seq, _part_begin, _part_begin, ring_idx))); } template template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::geometry_idx_from_part_idx( - IndexType part_idx) +multipolygon_range:: + geometry_idx_from_part_idx(IndexType part_idx) { return thrust::distance( - _geometry_begin, thrust::prev(thrust::upper_bound(_geometry_begin, _geometry_end, part_idx))); + _geometry_begin, + thrust::prev(thrust::upper_bound(thrust::seq, _geometry_begin, _geometry_end, part_idx))); } template template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::geometry_idx_from_point_idx( - IndexType point_idx) +multipolygon_range:: + geometry_idx_from_segment_idx(IndexType segment_idx) +{ + auto ring_idx = ring_idx_from_point_idx(segment_idx); + if (!is_valid_segment_id(segment_idx, ring_idx)) + return multipolygon_range:: + INVALID_INDEX; + return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx)); +} + +template +template +CUSPATIAL_HOST_DEVICE bool +multipolygon_range::is_valid_segment_id( + IndexType1 point_idx, IndexType2 ring_idx) { - return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx_from_part_idx(point_idx))); + if constexpr (std::is_signed_v) + return point_idx >= 0 && point_idx < (_ring_begin[ring_idx + 1] - 1); + else + return point_idx < (_ring_begin[ring_idx + 1] - 1); } template template CUSPATIAL_HOST_DEVICE auto -multipolygon_range::operator[]( +multipolygon_range::operator[]( IndexType multipolygon_idx) { return multipolygon_begin()[multipolygon_idx]; } +template +template +CUSPATIAL_HOST_DEVICE auto +multipolygon_range::get_segment( + IndexType segment_idx) +{ + return segment{_point_begin[segment_idx], _point_begin[segment_idx + 1]}; +} + } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index 331302348..b9efd64e5 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -63,6 +63,8 @@ class multipolygon_range { using point_t = iterator_value_type; using element_t = iterator_vec_base_type; + int64_t static constexpr INVALID_INDEX = -1; + multipolygon_range(GeometryIterator geometry_begin, GeometryIterator geometry_end, PartIterator part_begin, @@ -98,37 +100,52 @@ class multipolygon_range { /// Return the iterator to the one past the last multipolygon in the range. CUSPATIAL_HOST_DEVICE auto end() { return multipolygon_end(); } - - /// Given the index of a point, return the ring index where the point locates. + /// Given the index of a segment, return the geometry (multipolygon) index where the + /// segment locates. + /// Segment index is the index to the starting point of the segment. If the + /// index is the last point of the ring, then it is not a valid index. + /// This function returns multipolygon_range::INVALID_INDEX if the index is invalid. template - CUSPATIAL_HOST_DEVICE auto ring_idx_from_point_idx(IndexType point_idx); + CUSPATIAL_HOST_DEVICE auto geometry_idx_from_segment_idx(IndexType segment_idx); - /// Given the index of a ring, return the part (polygon) index - /// where the polygon locates. + /// Returns the `multipolygon_idx`th multipolygon in the range. template - CUSPATIAL_HOST_DEVICE auto part_idx_from_ring_idx(IndexType ring_idx); + CUSPATIAL_HOST_DEVICE auto operator[](IndexType multipolygon_idx); - /// Given the index of a part (polygon), return the geometry (multipolygon) index - /// where the polygon locates. - template - CUSPATIAL_HOST_DEVICE auto geometry_idx_from_part_idx(IndexType part_idx); + // template + // CUSPATIAL_HOST_DEVICE auto get_point(IndexType point_idx); - /// Given the index of a point, return the geometry (multipolygon) index where the - /// point locates. template - CUSPATIAL_HOST_DEVICE auto geometry_idx_from_point_idx(IndexType point_idx); - - /// Returns the `multipolygon_idx`th multipolygon in the range. - template - CUSPATIAL_HOST_DEVICE auto operator[](IndexType multipolygon_idx); + CUSPATIAL_HOST_DEVICE auto get_segment(IndexType segment_idx); protected: GeometryIterator _geometry_begin; GeometryIterator _geometry_end; PartIterator _part_begin; PartIterator _part_end; + RingIterator _ring_begin; + RingIterator _ring_end; VecIterator _point_begin; VecIterator _point_end; + + private: + /// Given the index of a point, return the ring index + /// where the point locates. + template + CUSPATIAL_HOST_DEVICE auto ring_idx_from_point_idx(IndexType point_idx); + + /// Given the index of a ring, return the part (polygon) index + /// where the ring locates. + template + CUSPATIAL_HOST_DEVICE auto part_idx_from_ring_idx(IndexType ring_idx); + + /// Given the index of a part (polygon), return the geometry (multipolygon) index + /// where the polygon locates. + template + CUSPATIAL_HOST_DEVICE auto geometry_idx_from_part_idx(IndexType part_idx); + + template + CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 ring_idx); }; } // namespace cuspatial From c2740702c2497774d05a377f91a58ee0bec1bee3 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:05:53 -0800 Subject: [PATCH 06/34] update multipoint_range class --- .../experimental/detail/ranges/multipoint_range.cuh | 7 +++++++ .../experimental/ranges/multipoint_range.cuh | 12 ++++++++++++ 2 files changed, 19 insertions(+) diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh index 08b3d5d5e..d3d8926d5 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh @@ -133,4 +133,11 @@ multipoint_range::geometry_idx_from_point_idx(Ind thrust::prev(thrust::upper_bound(thrust::seq, _geometry_begin, _geometry_end, idx))); } +template +template +CUSPATIAL_HOST_DEVICE auto multipoint_range::point(IndexType idx) +{ + return _points_begin[idx]; +} + } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh index 366f9c18e..1ee08b5cc 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh @@ -51,6 +51,8 @@ class multipoint_range { using point_t = iterator_value_type; using element_t = iterator_vec_base_type; + int32_t INVALID_IDX = -1; + /** * @brief Construct a new multipoint array object */ @@ -129,6 +131,16 @@ class multipoint_range { template CUSPATIAL_HOST_DEVICE auto operator[](IndexType idx); + /** + * @brief Returns the `idx`th point in the array. + * + * @tparam IndexType type of the index + * @param idx the index to the point + * @return a vec_2d object + */ + template + CUSPATIAL_HOST_DEVICE auto point(IndexType idx); + protected: /// Iterator to the start of the index array of start positions to each multipoint. GeometryIterator _geometry_begin; From 12ffa53f6a917a18bf103fc02c8eb11884b84464 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:06:36 -0800 Subject: [PATCH 07/34] update is_point_in_polygon usage with polygon_ref --- .../detail/is_point_in_polygon.cuh | 104 ------------------ .../detail/pairwise_point_in_polygon.cuh | 4 +- .../experimental/detail/point_in_polygon.cuh | 4 +- 3 files changed, 4 insertions(+), 108 deletions(-) delete mode 100644 cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh diff --git a/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh deleted file mode 100644 index b9c7d5ac6..000000000 --- a/cpp/include/cuspatial/experimental/detail/is_point_in_polygon.cuh +++ /dev/null @@ -1,104 +0,0 @@ -/* - * 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 - -#include - -namespace cuspatial { -namespace detail { - -/** - * @brief Kernel to test if a point is inside a polygon. - * - * Implemented based on Eric Haines's crossings-multiply algorithm: - * See "Crossings test" section of http://erich.realtimerendering.com/ptinpoly/ - * The improvement in addenda is also addopted to remove divisions in this kernel. - * - * 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 ::difference_type, - class Cart2dItDiffType = typename std::iterator_traits::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) -{ - using T = iterator_vec_base_type; - - 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; - - Cart2d b = poly_points_first[ring_end - 1]; - 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; - - // Points on the line segment are the same, so intersection is impossible. - // This is possible because we allow closed or unclosed polygons. - T constexpr zero = 0.0; - if (float_equal(run, zero) && float_equal(rise, zero)) continue; - - T rise_to_point = test_point.y - a.y; - - // colinearity test - T run_to_point = test_point.x - a.x; - is_colinear = float_equal(run * rise_to_point, run_to_point * rise); - if (is_colinear) { break; } - - y1_flag = a.y > test_point.y; - if (y1_flag != y0_flag) { - // Transform the following inequality to avoid division - // test_point.x < (run / rise) * rise_to_point + a.x - auto lhs = (test_point.x - a.x) * rise; - auto rhs = run * rise_to_point; - if (lhs < rhs != y1_flag) { point_is_within = not point_is_within; } - } - b = a; - y0_flag = y1_flag; - } - if (is_colinear) { - point_is_within = false; - break; - } - } - - return point_is_within; -} -} // namespace detail -} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/pairwise_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/pairwise_point_in_polygon.cuh index 8753367f6..f4659477a 100644 --- a/cpp/include/cuspatial/experimental/detail/pairwise_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/pairwise_point_in_polygon.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. @@ -18,7 +18,7 @@ #include #include -#include +#include #include #include diff --git a/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh index 5232cae1f..58737fe61 100644 --- a/cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_in_polygon.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. @@ -17,7 +17,7 @@ #pragma once #include -#include +#include #include #include From 749033340a81e816f30f95ce3a0b401d17ee8eb1 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:07:04 -0800 Subject: [PATCH 08/34] update multilinestring_range --- .../multilinestring_ref.cuh | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh index d1041818f..b8768d18f 100644 --- a/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh +++ b/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 #include #include @@ -22,6 +38,9 @@ struct to_linestring_functor { CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) { + printf("In to_linestring_functor: %d %d\n", + static_cast(part_begin[i]), + static_cast(part_begin[i + 1])); return linestring_ref{point_begin + part_begin[i], point_begin + part_begin[i + 1]}; } }; From f66528713dcb4eab82984648050c55c1e50d2253 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:07:22 -0800 Subject: [PATCH 09/34] add point to polygon kernel --- .../detail/point_polygon_distance.cuh | 118 +++++++++++++++--- .../experimental/point_polygon_distance.cuh | 2 +- 2 files changed, 102 insertions(+), 18 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index b51c252aa..7ab760b0b 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -17,50 +17,108 @@ #pragma once #include +#include #include #include +#include #include -#include -#include +#include +#include #include +#include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include -#include -#include -#include +#include #include namespace cuspatial { namespace detail { +template +struct point_in_multipolygon_test_functor { + MultiPointRange multipoints; + MultiPolygonRange multipolygons; + + point_in_multipolygon_test_functor(MultiPointRange multipoints, MultiPolygonRange multipolygons) + : multipoints(multipoints), multipolygons(multipolygons) + { + } + + template + uint8_t __device__ operator()(IndexType pidx) + { + printf("%d\n", static_cast(pidx)); + + auto point = thrust::raw_reference_cast(multipoints.point(pidx)); + + printf("%f, %f\n", point.x, point.y); + + auto geometry_idx = multipoints.geometry_idx_from_point_idx(pidx); + + printf("%d\n", static_cast(geometry_idx)); + + bool intersects = false; + for (auto polygon : multipolygons[geometry_idx]) { + printf("here\n"); + intersects = intersects || is_point_in_polygon(point, polygon); + } + + return static_cast(intersects); + } +}; + /** * @brief Kernel to compute the distance between pairs of point and polygon. */ -template +template void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoints, MultiPolygonRange multipolygons, + IntersectionRange intersects, OutputIterator distances) { using T = typename MultiPointRange::element_t; + T dist_squared = std::numeric_limits::max(); for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multipolygons.num_points(); idx += gridDim.x * blockDim.x) { - auto geometry_idx = multipolygons.geometry_idx_from_point_idx(idx); + auto geometry_idx = multipolygons.geometry_idx_from_segment_idx(idx); + if (geometry_idx == MultiPolygonRange::INVALID_INDEX) continue; + + if (intersects[geometry_idx]) { + // TODO: only the leading thread of the pair need to store the result, atomics is not needed. + atomicMin(&distances[geometry_idx], T{0.0}); + continue; + } + + printf("In distance kernel: %d %d %d", + static_cast(idx), + static_cast(geometry_idx), + static_cast(intersects.size())); + + auto [a, b] = multipolygons.get_segment(idx); + for (vec_2d point : multipoints[geometry_idx]) { + dist_squared = min(dist_squared, point_to_segment_distance_squared(point, a, b)); + } + + atomicMin(&distances[geometry_idx], sqrt(dist_squared)); } } } // namespace detail + template OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, - MultiPolygonRange multipoiygons, + MultiPolygonRange multipolygons, OutputIt distances_first, rmm::cuda_stream_view stream) { @@ -76,12 +134,38 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, CUSPATIAL_EXPECTS(multipoints.size() == multipolygons.size(), "Must have the same number of input rows."); - auto [threads_per_block, n_blocks] = grid_id(multipolygons.num_points()); + auto multipoint_intersects = [&]() { + rmm::device_uvector point_intersects(multipoints.num_points(), stream); + + thrust::tabulate(rmm::exec_policy(stream), + point_intersects.begin(), + point_intersects.end(), + detail::point_in_multipolygon_test_functor{multipoints, multipolygons}); + + rmm::device_uvector multipoint_intersects(multipoints.num_multipoints(), stream); + auto offset_as_key_it = detail::make_counting_transform_iterator( + 0, offsets_to_keys_functor{multipoints.offsets_begin(), multipoints.offsets_end()}); - pairwise_point_polygon_distance_kernel<<>>( - multipoints, multipolygons, distances_first); + thrust::reduce_by_key(rmm::exec_policy(stream), + offset_as_key_it, + offset_as_key_it + multipoints.num_points(), + point_intersects.begin(), + thrust::make_discard_iterator(), + multipoint_intersects.begin(), + thrust::logical_or()); + + return multipoint_intersects; + }(); + + auto [threads_per_block, n_blocks] = grid_1d(multipolygons.num_points()); + + detail:: + pairwise_point_polygon_distance_kernel<<>>( + multipoints, multipolygons, multipoint_intersects.begin(), distances_first); CUSPATIAL_CHECK_CUDA(stream.value()); + + return distances_first + multipoints.size(); } } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh index c5394e902..95c3a7d2c 100644 --- a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh @@ -33,7 +33,7 @@ namespace cuspatial { */ template OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, - MultiPolygonRangeB multipoiygons, + MultiPolygonRange multipoiygons, OutputIt distances_first, rmm::cuda_stream_view stream = rmm::cuda_stream_default); } // namespace cuspatial From 291f6e679eea9e87fb87646a09593d07bf863326 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:08:10 -0800 Subject: [PATCH 10/34] add segment deduction guide --- cpp/include/cuspatial/experimental/geometry/segment.cuh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cpp/include/cuspatial/experimental/geometry/segment.cuh b/cpp/include/cuspatial/experimental/geometry/segment.cuh index 30d9240d4..8667e174c 100644 --- a/cpp/include/cuspatial/experimental/geometry/segment.cuh +++ b/cpp/include/cuspatial/experimental/geometry/segment.cuh @@ -19,6 +19,8 @@ #include #include +#include + #include namespace cuspatial { @@ -61,4 +63,7 @@ class alignas(sizeof(Vertex)) segment { template segment(vec_2d a, vec_2d b) -> segment>; +template +segment(thrust::device_reference> a, thrust::device_reference> b) + -> segment>; } // namespace cuspatial From efa68830493511998064a00fba9f3b293528865b Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:08:30 -0800 Subject: [PATCH 11/34] add owning object type to vector factories --- .../cuspatial_test/vector_factories.cuh | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh index 00e4c5e07..7747dfe33 100644 --- a/cpp/include/cuspatial_test/vector_factories.cuh +++ b/cpp/include/cuspatial_test/vector_factories.cuh @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -50,6 +51,63 @@ auto make_device_uvector(std::initializer_list inl, return res; } +/** + * @brief Owning object of a multipolygon array following geoarrow layout. + * + * @tparam GeometryArray Array type of geometry offset array + * @tparam PartArray Array type of part offset array + * @tparam RingArray Array type of ring offset array + * @tparam CoordinateArray Array type of coordinate array + */ +template +class multipolygon_array { + public: + multipolygon_array(GeometryArray geometry_offsets_array, + PartArray part_offsets_array, + RingArray ring_offsets_array, + CoordinateArray coordinate_offsets_array) + : _geometry_offsets_array(geometry_offsets_array), + _part_offsets_array(part_offsets_array), + _ring_offsets_array(ring_offsets_array), + _coordinate_offsets_array(coordinate_offsets_array) + { + } + + /// Return the number of multilinestrings + auto size() { return _geometry_offsets_array.size() - 1; } + + /// Return range object of the multilinestring array + auto range() + { + return multipolygon_range(_geometry_offsets_array.begin(), + _geometry_offsets_array.end(), + _part_offsets_array.begin(), + _part_offsets_array.end(), + _ring_offsets_array.begin(), + _ring_offsets_array.end(), + _coordinate_offsets_array.begin(), + _coordinate_offsets_array.end()); + } + + protected: + GeometryArray _geometry_offsets_array; + PartArray _part_offsets_array; + RingArray _ring_offsets_array; + CoordinateArray _coordinate_offsets_array; +}; + +template +auto make_multipolygon_array(std::initializer_list geometry_inl, + std::initializer_list part_inl, + std::initializer_list ring_inl, + std::initializer_list> coord_inl) +{ + return multipolygon_array{make_device_vector(geometry_inl), + make_device_vector(part_inl), + make_device_vector(ring_inl), + make_device_vector(coord_inl)}; +} + /** * @brief Owning object of a multilinestring array following geoarrow layout. * From 23146ef32dbac357085464c11b47c4a73af56b8b Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:08:58 -0800 Subject: [PATCH 12/34] add tests --- cpp/tests/CMakeLists.txt | 3 + .../spatial/point_polygon_distance_test.cu | 90 +++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 cpp/tests/experimental/spatial/point_polygon_distance_test.cu diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 6182e1238..821c7881b 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -133,6 +133,9 @@ ConfigureTest(POINT_DISTANCE_TEST_EXP ConfigureTest(POINT_LINESTRING_DISTANCE_TEST_EXP experimental/spatial/point_linestring_distance_test.cu) +ConfigureTest(POINT_POLYGON_DISTANCE_TEST_EXP + experimental/spatial/point_polygon_distance_test.cu) + ConfigureTest(HAUSDORFF_TEST_EXP experimental/spatial/hausdorff_test.cu) diff --git a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu new file mode 100644 index 000000000..987cefa04 --- /dev/null +++ b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu @@ -0,0 +1,90 @@ +/* + * 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. + */ + +#include +#include + +#include +#include + +#include +#include + +#include + +using namespace cuspatial; +using namespace cuspatial::test; + +template +struct PairwisePointPolygonDistanceTest : public ::testing::Test { + rmm::cuda_stream_view stream() { return rmm::cuda_stream_default; } + rmm::mr::device_memory_resource* mr() { return rmm::mr::get_current_device_resource(); } + + void run_single(std::initializer_list>> multipoints, + std::initializer_list multipolygon_geometry_offsets, + std::initializer_list multipolygon_part_offsets, + std::initializer_list multipolygon_ring_offsets, + std::initializer_list> multipolygon_coordinates, + std::initializer_list expected) + { + auto d_multipoints = make_multipoints_array(multipoints); + auto d_multipolygons = make_multipolygon_array(multipolygon_geometry_offsets, + multipolygon_part_offsets, + multipolygon_ring_offsets, + multipolygon_coordinates); + + auto got = rmm::device_uvector(d_multipoints.size(), stream()); + + auto ret = pairwise_point_polygon_distance( + d_multipoints.range(), d_multipolygons.range(), got.begin(), stream()); + + auto d_expected = make_device_vector(expected); + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); + EXPECT_EQ(ret, got.end()); + } +}; + +using TestTypes = ::testing::Types; + +TYPED_TEST_CASE(PairwisePointPolygonDistanceTest, TestTypes); + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{0, 0}}}, + {0, 1}, + {0, 1}, + {0, 5}, + {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}}, + {0.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing2) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{0, 2}}}, + {0, 1}, + {0, 1}, + {0, 5}, + {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}}, + {1.0}); +} From ead160a14c111f683c9c550944b95b7db7bfc12a Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 7 Mar 2023 11:09:21 -0800 Subject: [PATCH 13/34] add helper files --- .../detail/utility/offset_to_keys.cuh | 51 ++++++++ .../detail/algorithm/is_point_in_polygon.cuh | 112 ++++++++++++++++++ 2 files changed, 163 insertions(+) create mode 100644 cpp/include/cuspatial/detail/utility/offset_to_keys.cuh create mode 100644 cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh diff --git a/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh b/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh new file mode 100644 index 000000000..70665b4db --- /dev/null +++ b/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh @@ -0,0 +1,51 @@ +/* + * 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 +#include + +/** @brief Given list offset and row `i`, return a unique key that represent the list of `i`. + * + * The key is computed by performing a `upper_bound` search with `i` in the offset array. + * Then subtracts the position with the start of offset array. + * + * Example: + * offset: 0 0 0 1 3 4 4 4 + * i: 0 1 2 3 + * key: 3 4 4 5 + * + * Note that the values of `key`, {offset[3], offset[4], offset[5]} denotes the ending + * position of the first 3 non-empty list. + */ +template +struct offsets_to_keys_functor { + Iterator _offsets_begin; + Iterator _offsets_end; + + offsets_to_keys_functor(Iterator offset_begin, Iterator offset_end) + : _offsets_begin(offset_begin), _offsets_end(offset_end) + { + } + + template + IndexType __device__ operator()(IndexType i) + { + return thrust::distance(_offsets_begin, + thrust::upper_bound(thrust::seq, _offsets_begin, _offsets_end, i)); + } +}; diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh new file mode 100644 index 000000000..e489c22dd --- /dev/null +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -0,0 +1,112 @@ +/* + * 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 +#include + +#include +#include + +namespace cuspatial { +namespace detail { + +/** + * @brief Kernel to test if a point is inside a polygon. + * + * Implemented based on Eric Haines's crossings-multiply algorithm: + * See "Crossings test" section of http://erich.realtimerendering.com/ptinpoly/ + * The improvement in addenda is also addopted to remove divisions in this kernel. + * + * 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 +__device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonRef const& polygon) +{ + bool point_is_within = false; + bool is_colinear = false; + printf("Polygon num rings: %d\n", static_cast(polygon.num_rings())); + for (auto ring : polygon) { + printf("here in ring\n"); + bool y0_flag = ring.segment(ring.num_segments() - 1).v2.y > test_point.y; + bool y1_flag; + for (auto [a, b] : ring) { + printf("here in segment (%f, %f) -> (%f, %f)\n", a.x, a.y, b.x, b.y); + // 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. + T constexpr zero = 0.0; + if (float_equal(run, zero) && float_equal(rise, zero)) continue; + + T rise_to_point = test_point.y - a.y; + + // colinearity test + T run_to_point = test_point.x - a.x; + is_colinear = float_equal(run * rise_to_point, run_to_point * rise); + if (is_colinear) { break; } + + y1_flag = a.y > test_point.y; + if (y1_flag != y0_flag) { + // Transform the following inequality to avoid division + // test_point.x < (run / rise) * rise_to_point + a.x + auto lhs = (test_point.x - a.x) * rise; + auto rhs = run * rise_to_point; + if (lhs < rhs != y1_flag) { point_is_within = not point_is_within; } + } + b = a; + y0_flag = y1_flag; + } + if (is_colinear) { + point_is_within = false; + break; + } + } + + printf("Exiting pip.\n"); + + return point_is_within; +} + +template ::difference_type, + class Cart2dItDiffType = typename std::iterator_traits::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{ring_offsets_first + poly_begin, + ring_offsets_first + poly_end, + poly_points_first, + poly_points_first + num_poly_points}; + return is_point_in_polygon(test_point, polygon); +} + +} // namespace detail +} // namespace cuspatial From 09bd35f2bd8ead3a6fa525bb28ccfb4cf5981ecd Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 13:15:27 -0800 Subject: [PATCH 14/34] add more tests --- .../spatial/point_polygon_distance_test.cu | 337 ++++++++++++++++++ 1 file changed, 337 insertions(+) diff --git a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu index 987cefa04..519ed6819 100644 --- a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu +++ b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -61,6 +62,20 @@ using TestTypes = ::testing::Types; TYPED_TEST_CASE(PairwisePointPolygonDistanceTest, TestTypes); +TYPED_TEST(PairwisePointPolygonDistanceTest, ZeroPairs) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + std::initializer_list>{}, + {0}, + {0}, + {0}, + std::initializer_list

{}, + std::initializer_list{}); +} + TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing) { using T = TypeParam; @@ -88,3 +103,325 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing2) {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}}, {1.0}); } + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{0, 0}}}, + {0, 1}, + {0, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{2, -2}, + P{2, 2}, + P{-2, 2}, + P{-2, -2}, + P{-1, -1}, + P{1, -1}, + P{1, 1}, + P{-1, 1}, + P{-1, -1}, + }, + {1.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings2) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{1.5, 0}}}, + {0, 1}, + {0, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{2, -2}, + P{2, 2}, + P{-2, 2}, + P{-2, -2}, + P{-1, -1}, + P{1, -1}, + P{1, 1}, + P{-1, 1}, + P{-1, -1}, + }, + {0.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings3) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{3, 0}}}, + {0, 1}, + {0, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{2, -2}, + P{2, 2}, + P{-2, 2}, + P{-2, -2}, + P{-1, -1}, + P{1, -1}, + P{1, 1}, + P{-1, 1}, + P{-1, -1}, + }, + {1.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{1, 1}}}, + {0, 2}, + {0, 1, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{0, -2}, + P{0, 0}, + P{-2, 0}, + P{-2, -2}, + P{0, 0}, + P{2, 0}, + P{2, 2}, + P{0, 2}, + P{0, 0}, + }, + {0.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing2) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{-1, -1}}}, + {0, 2}, + {0, 1, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{0, -2}, + P{0, 0}, + P{-2, 0}, + P{-2, -2}, + P{0, 0}, + P{2, 0}, + P{2, 2}, + P{0, 2}, + P{0, 0}, + }, + {0.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing3) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{-1, 0.5}}}, + {0, 2}, + {0, 1, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{0, -2}, + P{0, 0}, + P{-2, 0}, + P{-2, -2}, + P{0, 0}, + P{2, 0}, + P{2, 2}, + P{0, 2}, + P{0, 0}, + }, + {0.5}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing4) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{-0.3, 1}}}, + {0, 2}, + {0, 1, 2}, + {0, 5, 10}, + { + P{-2, -2}, + P{0, -2}, + P{0, 0}, + P{-2, 0}, + P{-2, -2}, + P{0, 0}, + P{2, 0}, + P{2, 2}, + P{0, 2}, + P{0, 0}, + }, + {0.3}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairOnePolygonOneRing) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{-0.6, -0.6}}, {P{0, 0}}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + { + P{-1, -1}, + P{0, 0}, + P{0, 1}, + P{-1, -1}, + P{1, 1}, + P{1, 0}, + P{2, 2}, + P{1, 1}, + }, + {0.0, 1.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairTwoPolygonTwoRing) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{2.5, 3}}, {P{-2.5, -1.5}}}, + {0, 1, 2}, + {0, 2, 4}, + {0, 5, 10, 14, 18}, + { + P{0, 0}, + P{3, 0}, + P{3, 3}, + P{0, 3}, + P{0, 0}, + P{1, 1}, + P{2, 1}, + P{2, 2}, + P{1, 2}, + P{1, 1}, + + P{0, 0}, + P{-3, 0}, + P{-3, -3}, + P{0, 0}, + P{-1, -1}, + P{-2, -1}, + P{-2, -2}, + P{-1, -1}, + }, + {0.0, 0.5}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, ThreePolygons) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{1, 1}}, {P{2, 2}}, {P{1.5, 0}}}, + {0, 1, 2, 3}, + {0, 2, 3, 6}, + {0, 4, 8, 12, 17, 22, 27}, + {// POLYGON ((0 1, -1 -1, 1 -1, 0 1), (0 0.5, 0.5 -0.5, -0.5 -0.5, 0 0.5)) + P{0, 1}, + P{-1, -1}, + P{1, -1}, + P{0, 1}, + P{0, 0.5}, + P{0.5, -0.5}, + P{-0.5, -0.5}, + P{0, 0.5}, + // POLYGON ((1 1, 1 2, 2 1, 1 1)) + P{1, 1}, + P{1, 2}, + P{2, 1}, + P{1, 1}, + // POLYGON ( + // (-3 -3, 3 -3, 3 3, -3 3, -3 -3), + // (-2 -2, -1 -2, -1 2, -2 2, -2 -2), + // (2 2, 2 -2, 1 -2, 1 2, 2 2) + // ) + P{-3, -3}, + P{3, -3}, + P{3, 3}, + P{-3, 3}, + P{-3, -3}, + P{-2, -2}, + P{-1, -2}, + P{-1, 2}, + P{-2, 2}, + P{-2, -2}, + P{2, 2}, + P{2, -2}, + P{1, -2}, + P{1, 2}, + P{2, 2}}, + {0.894427190999916, 0.7071067811865476, 0.5}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{0, 3}, P{2, 0}}}, + {0, 1}, + {0, 1}, + {0, 5}, + {P{0, 1}, P{-1, -1}, P{1, -1}, P{0, 1}}, + {1.3416407864998738}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon2) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + {{P{0, 3}, P{0, 0}}}, + {0, 1}, + {0, 1}, + {0, 5}, + {P{0, 1}, P{-1, -1}, P{1, -1}, P{0, 1}}, + {0.0}); +} + +TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairMultiPointOnePolygon2) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST( + this->run_single, + {{P{0, 2}, P{0, 0}}, {P{1, 1}, P{-1, -1}}}, + {0, 1, 2}, + {0, 1, 2}, + {0, 5, 9}, + {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}, P{-1, 1}, P{1, 1}, P{0, -1}, P{-1, 1}}, + {0.0, 0.0}); +} From 92760d13d330f63ca5f4d999220edfa3233e4427 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 13:23:03 -0800 Subject: [PATCH 15/34] bug fixes --- .../detail/algorithm/is_point_in_polygon.cuh | 23 ++++++++++-- .../detail/geometry/linestring_ref.cuh | 14 +++++++- .../detail/point_polygon_distance.cuh | 35 +++++++++++++++---- .../detail/ranges/multipolygon_range.cuh | 8 +++-- .../experimental/geometry/linestring_ref.cuh | 5 +++ .../ranges/multipolygon_range.cuh | 3 -- 6 files changed, 73 insertions(+), 15 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh index e489c22dd..61b80b365 100644 --- a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -45,9 +45,18 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR printf("Polygon num rings: %d\n", static_cast(polygon.num_rings())); for (auto ring : polygon) { printf("here in ring\n"); - bool y0_flag = ring.segment(ring.num_segments() - 1).v2.y > test_point.y; + auto last_segment = ring.segment(ring.num_segments() - 1); + printf("Last segment is: (%f %f)->(%f %f)\n", + last_segment.v1.x, + last_segment.v1.y, + last_segment.v2.x, + last_segment.v2.y); + + auto b = last_segment.v2; + bool y0_flag = b.y > test_point.y; bool y1_flag; - for (auto [a, b] : ring) { + for (auto it = ring.point_begin(); it != ring.point_end(); ++it) { + vec_2d a = *it; printf("here in segment (%f, %f) -> (%f, %f)\n", a.x, a.y, b.x, b.y); // for each line segment, including the segment between the last and first vertex T run = b.x - a.x; @@ -65,7 +74,12 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR is_colinear = float_equal(run * rise_to_point, run_to_point * rise); if (is_colinear) { break; } + // y0_flag = a.y > test_point.y; y1_flag = a.y > test_point.y; + printf("\t y0_flag: %d, y1_flag: %d, point_is_within: %d\n", + static_cast(y0_flag), + static_cast(y1_flag), + static_cast(point_is_within)); if (y1_flag != y0_flag) { // Transform the following inequality to avoid division // test_point.x < (run / rise) * rise_to_point + a.x @@ -82,11 +96,14 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR } } - printf("Exiting pip.\n"); + printf("Exiting pip. Result: %d\n", static_cast(point_is_within)); return point_is_within; } +/** + * @brief Compatibility layer with non-OOP style input + */ template CUSPATIAL_HOST_DEVICE auto linestring_ref::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::segment_end() const return segment_begin() + num_segments(); } +template +CUSPATIAL_HOST_DEVICE auto linestring_ref::point_begin() const +{ + return _point_begin; +} + +template +CUSPATIAL_HOST_DEVICE auto linestring_ref::point_end() const +{ + return _point_end; +} + template template CUSPATIAL_HOST_DEVICE auto linestring_ref::segment(IndexType i) const diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index 7ab760b0b..efd825cb3 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -16,11 +16,14 @@ #pragma once +#include + #include #include #include #include #include +#include #include #include #include @@ -36,6 +39,7 @@ #include #include +#include #include namespace cuspatial { @@ -70,6 +74,7 @@ struct point_in_multipolygon_test_functor { intersects = intersects || is_point_in_polygon(point, polygon); } + printf("Intersect: %d\n", static_cast(intersects)); return static_cast(intersects); } }; @@ -88,25 +93,32 @@ void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoin { using T = typename MultiPointRange::element_t; - T dist_squared = std::numeric_limits::max(); for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < multipolygons.num_points(); idx += gridDim.x * blockDim.x) { auto geometry_idx = multipolygons.geometry_idx_from_segment_idx(idx); if (geometry_idx == MultiPolygonRange::INVALID_INDEX) continue; + printf("Intesects? %d Geometry_idx: %d\n", + static_cast(intersects[geometry_idx]), + static_cast(geometry_idx)); + if (intersects[geometry_idx]) { // TODO: only the leading thread of the pair need to store the result, atomics is not needed. + printf("In intersects, idx: %d\n", static_cast(idx)); atomicMin(&distances[geometry_idx], T{0.0}); continue; } - printf("In distance kernel: %d %d %d", + printf("In distance kernel: point_idx: %d segment_idx: %d\n", static_cast(idx), - static_cast(geometry_idx), - static_cast(intersects.size())); + static_cast(geometry_idx)); - auto [a, b] = multipolygons.get_segment(idx); + T dist_squared = std::numeric_limits::max(); + auto [a, b] = multipolygons.get_segment(idx); for (vec_2d point : multipoints[geometry_idx]) { + printf("point: %f %f\n", point.x, point.y); + printf("segment: (%f %f) -> (%f %f)\n", a.x, a.y, b.x, b.y); + printf("dist: %f\n", point_to_segment_distance_squared(point, a, b)); dist_squared = min(dist_squared, point_to_segment_distance_squared(point, a, b)); } @@ -134,6 +146,8 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, CUSPATIAL_EXPECTS(multipoints.size() == multipolygons.size(), "Must have the same number of input rows."); + if (multipoints.size() == 0) return distances_first; + auto multipoint_intersects = [&]() { rmm::device_uvector point_intersects(multipoints.num_points(), stream); @@ -142,7 +156,10 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, point_intersects.end(), detail::point_in_multipolygon_test_functor{multipoints, multipolygons}); + // TODO: optimize when input is not a multipolygon rmm::device_uvector multipoint_intersects(multipoints.num_multipoints(), stream); + detail::zero_data_async(multipoint_intersects.begin(), multipoint_intersects.end(), stream); + auto offset_as_key_it = detail::make_counting_transform_iterator( 0, offsets_to_keys_functor{multipoints.offsets_begin(), multipoints.offsets_end()}); @@ -157,10 +174,16 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, return multipoint_intersects; }(); + thrust::fill(rmm::exec_policy(stream), + distances_first, + distances_first + multipoints.size(), + std::numeric_limits::max()); auto [threads_per_block, n_blocks] = grid_1d(multipolygons.num_points()); + std::cout << "Size of multipoint intersects: " << multipoint_intersects.size() << std::endl; + detail:: - pairwise_point_polygon_distance_kernel<<>>( + pairwise_point_polygon_distance_kernel<<>>( multipoints, multipolygons, multipoint_intersects.begin(), distances_first); CUSPATIAL_CHECK_CUDA(stream.value()); diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index b44420532..33b7a512e 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -203,8 +203,7 @@ multipolygon_range:: part_idx_from_ring_idx(IndexType ring_idx) { return thrust::distance( - _part_begin, - thrust::prev(thrust::upper_bound(thrust::seq, _part_begin, _part_begin, ring_idx))); + _part_begin, thrust::prev(thrust::upper_bound(thrust::seq, _part_begin, _part_end, ring_idx))); } template :: geometry_idx_from_segment_idx(IndexType segment_idx) { + printf("segment_idx: %d\n", static_cast(segment_idx)); auto ring_idx = ring_idx_from_point_idx(segment_idx); + printf("ring_idx: %d\n", static_cast(ring_idx)); if (!is_valid_segment_id(segment_idx, ring_idx)) return multipolygon_range:: INVALID_INDEX; + + auto part_idx = part_idx_from_ring_idx(ring_idx); + printf("part_idx: %d\n", static_cast(part_idx)); return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx)); } diff --git a/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh b/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh index 0ebad40f4..707b9e0cc 100644 --- a/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh @@ -38,6 +38,11 @@ class linestring_ref { /// Return iterator to one past the last segment CUSPATIAL_HOST_DEVICE auto segment_end() const; + /// Return iterator to the first point of the linestring + CUSPATIAL_HOST_DEVICE auto point_begin() const; + /// Return iterator to one past the last point + CUSPATIAL_HOST_DEVICE auto point_end() const; + /// Return iterator to the first segment of the linestring CUSPATIAL_HOST_DEVICE auto begin() const { return segment_begin(); } /// Return iterator to one past the last segment diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index b9efd64e5..8a7300213 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -112,9 +112,6 @@ class multipolygon_range { template CUSPATIAL_HOST_DEVICE auto operator[](IndexType multipolygon_idx); - // template - // CUSPATIAL_HOST_DEVICE auto get_point(IndexType point_idx); - template CUSPATIAL_HOST_DEVICE auto get_segment(IndexType segment_idx); From 8acb5dc48127664d92b5f89370616aa1dfa25b6e Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 13:47:29 -0800 Subject: [PATCH 16/34] cleanups --- .../detail/algorithm/is_point_in_polygon.cuh | 14 --------- .../multilinestring_ref.cuh | 3 -- .../detail/point_polygon_distance.cuh | 30 ++----------------- .../detail/ranges/multipolygon_range.cuh | 27 ++++++++--------- .../ranges/multipolygon_range.cuh | 7 +++++ 5 files changed, 22 insertions(+), 59 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh index 61b80b365..135011122 100644 --- a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -42,22 +42,14 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR { bool point_is_within = false; bool is_colinear = false; - printf("Polygon num rings: %d\n", static_cast(polygon.num_rings())); for (auto ring : polygon) { - printf("here in ring\n"); auto last_segment = ring.segment(ring.num_segments() - 1); - printf("Last segment is: (%f %f)->(%f %f)\n", - last_segment.v1.x, - last_segment.v1.y, - last_segment.v2.x, - last_segment.v2.y); auto b = last_segment.v2; bool y0_flag = b.y > test_point.y; bool y1_flag; for (auto it = ring.point_begin(); it != ring.point_end(); ++it) { vec_2d a = *it; - printf("here in segment (%f, %f) -> (%f, %f)\n", a.x, a.y, b.x, b.y); // 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; @@ -76,10 +68,6 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR // y0_flag = a.y > test_point.y; y1_flag = a.y > test_point.y; - printf("\t y0_flag: %d, y1_flag: %d, point_is_within: %d\n", - static_cast(y0_flag), - static_cast(y1_flag), - static_cast(point_is_within)); if (y1_flag != y0_flag) { // Transform the following inequality to avoid division // test_point.x < (run / rise) * rise_to_point + a.x @@ -96,8 +84,6 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR } } - printf("Exiting pip. Result: %d\n", static_cast(point_is_within)); - return point_is_within; } diff --git a/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh index b8768d18f..b4c76ec32 100644 --- a/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh +++ b/cpp/include/cuspatial/experimental/detail/geometry_collection/multilinestring_ref.cuh @@ -38,9 +38,6 @@ struct to_linestring_functor { CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) { - printf("In to_linestring_functor: %d %d\n", - static_cast(part_begin[i]), - static_cast(part_begin[i + 1])); return linestring_ref{point_begin + part_begin[i], point_begin + part_begin[i + 1]}; } }; diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index efd825cb3..21bc5a332 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -58,23 +58,13 @@ struct point_in_multipolygon_test_functor { template uint8_t __device__ operator()(IndexType pidx) { - printf("%d\n", static_cast(pidx)); - - auto point = thrust::raw_reference_cast(multipoints.point(pidx)); - - printf("%f, %f\n", point.x, point.y); - + auto point = thrust::raw_reference_cast(multipoints.point(pidx)); auto geometry_idx = multipoints.geometry_idx_from_point_idx(pidx); - printf("%d\n", static_cast(geometry_idx)); - bool intersects = false; for (auto polygon : multipolygons[geometry_idx]) { - printf("here\n"); intersects = intersects || is_point_in_polygon(point, polygon); } - - printf("Intersect: %d\n", static_cast(intersects)); return static_cast(intersects); } }; @@ -98,27 +88,15 @@ void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoin auto geometry_idx = multipolygons.geometry_idx_from_segment_idx(idx); if (geometry_idx == MultiPolygonRange::INVALID_INDEX) continue; - printf("Intesects? %d Geometry_idx: %d\n", - static_cast(intersects[geometry_idx]), - static_cast(geometry_idx)); - if (intersects[geometry_idx]) { - // TODO: only the leading thread of the pair need to store the result, atomics is not needed. - printf("In intersects, idx: %d\n", static_cast(idx)); - atomicMin(&distances[geometry_idx], T{0.0}); + if (multipolygons.is_first_point_of_multipolygon(idx, geometry_idx)) + distances[geometry_idx] = T{0.0}; continue; } - printf("In distance kernel: point_idx: %d segment_idx: %d\n", - static_cast(idx), - static_cast(geometry_idx)); - T dist_squared = std::numeric_limits::max(); auto [a, b] = multipolygons.get_segment(idx); for (vec_2d point : multipoints[geometry_idx]) { - printf("point: %f %f\n", point.x, point.y); - printf("segment: (%f %f) -> (%f %f)\n", a.x, a.y, b.x, b.y); - printf("dist: %f\n", point_to_segment_distance_squared(point, a, b)); dist_squared = min(dist_squared, point_to_segment_distance_squared(point, a, b)); } @@ -180,8 +158,6 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, std::numeric_limits::max()); auto [threads_per_block, n_blocks] = grid_1d(multipolygons.num_points()); - std::cout << "Size of multipoint intersects: " << multipoint_intersects.size() << std::endl; - detail:: pairwise_point_polygon_distance_kernel<<>>( multipoints, multipolygons, multipoint_intersects.begin(), distances_first); diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index 33b7a512e..0228d1de8 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -66,18 +66,6 @@ struct to_multipolygon_functor { CUSPATIAL_HOST_DEVICE auto operator()(difference_type i) { - auto new_part_begin = _part_begin + _geometry_begin[i]; - auto new_part_end = thrust::next(_part_begin, _geometry_begin[i + 1] + 1); - - printf( - "In to_multipolygon_functor: %d %d %d\n\t new_part_begin_val: %d, " - "new_part_end_dist_from_begin: %d\n", - static_cast(i), - static_cast(_geometry_begin[i]), - static_cast(_geometry_begin[i + 1]), - static_cast(*new_part_begin), - static_cast(thrust::distance(new_part_end, new_part_begin))); - return multipolygon_ref{_part_begin + _geometry_begin[i], thrust::next(_part_begin, _geometry_begin[i + 1] + 1), _ring_begin, @@ -229,15 +217,12 @@ CUSPATIAL_HOST_DEVICE auto multipolygon_range:: geometry_idx_from_segment_idx(IndexType segment_idx) { - printf("segment_idx: %d\n", static_cast(segment_idx)); auto ring_idx = ring_idx_from_point_idx(segment_idx); - printf("ring_idx: %d\n", static_cast(ring_idx)); if (!is_valid_segment_id(segment_idx, ring_idx)) return multipolygon_range:: INVALID_INDEX; auto part_idx = part_idx_from_ring_idx(ring_idx); - printf("part_idx: %d\n", static_cast(part_idx)); return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx)); } @@ -280,4 +265,16 @@ multipolygon_range::g return segment{_point_begin[segment_idx], _point_begin[segment_idx + 1]}; } +template +template +CUSPATIAL_HOST_DEVICE bool +multipolygon_range:: + is_first_point_of_multipolygon(IndexType1 point_idx, IndexType2 geometry_idx) +{ + return point_idx == _ring_begin[_part_begin[_geometry_begin[geometry_idx]]]; +} + } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index 8a7300213..775c7dba4 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -112,9 +112,16 @@ class multipolygon_range { template CUSPATIAL_HOST_DEVICE auto operator[](IndexType multipolygon_idx); + /// Returns the `segment_idx`th segment in the multipolygon range. template CUSPATIAL_HOST_DEVICE auto get_segment(IndexType segment_idx); + /// Returns `true` if `point_idx`th point is the first point of its + /// multipolygon + template + CUSPATIAL_HOST_DEVICE bool is_first_point_of_multipolygon(IndexType1 point_idx, + IndexType2 geometry_idx); + protected: GeometryIterator _geometry_begin; GeometryIterator _geometry_end; From a2b94fe7410d0e4a32940ee52a13db59955e0e27 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 14:38:32 -0800 Subject: [PATCH 17/34] fix tests --- cpp/tests/experimental/spatial/point_polygon_distance_test.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu index 519ed6819..ce5f78b23 100644 --- a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu +++ b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu @@ -308,7 +308,7 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairTwoPolygonTwoRing) using P = vec_2d; CUSPATIAL_RUN_TEST(this->run_single, - {{P{2.5, 3}}, {P{-2.5, -1.5}}}, + {{P{2.5, 3}}, {P{-1.75, -1.5}}}, {0, 1, 2}, {0, 2, 4}, {0, 5, 10, 14, 18}, @@ -333,7 +333,7 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairTwoPolygonTwoRing) P{-2, -2}, P{-1, -1}, }, - {0.0, 0.5}); + {0.0, 0.17677669529663687}); } TYPED_TEST(PairwisePointPolygonDistanceTest, ThreePolygons) From 46a67fe5612967a44d7a1708b0b742e0997a7e02 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 14:39:07 -0800 Subject: [PATCH 18/34] optimize single point range input --- .../experimental/detail/point_polygon_distance.cuh | 7 ++++++- .../experimental/detail/ranges/multipoint_range.cuh | 6 ++++++ .../cuspatial/experimental/ranges/multipoint_range.cuh | 5 +++++ 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index 21bc5a332..b217c5f35 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -126,6 +126,9 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, if (multipoints.size() == 0) return distances_first; + // Compute whether each multipoint intersects with the corresponding multipolygon. + // First, compute the point-multipolygon intersection. Then use reduce-by-key to + // compute the multipoint-multipolygon intersection. auto multipoint_intersects = [&]() { rmm::device_uvector point_intersects(multipoints.num_points(), stream); @@ -134,7 +137,9 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, point_intersects.end(), detail::point_in_multipolygon_test_functor{multipoints, multipolygons}); - // TODO: optimize when input is not a multipolygon + // `multipoints` contains only single points, no need to reduce. + if (multipoints.is_single_point_range()) return point_intersects; + rmm::device_uvector multipoint_intersects(multipoints.num_multipoints(), stream); detail::zero_data_async(multipoint_intersects.begin(), multipoint_intersects.end(), stream); diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh index d3d8926d5..7b7ce5243 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh @@ -140,4 +140,10 @@ CUSPATIAL_HOST_DEVICE auto multipoint_range::poin return _points_begin[idx]; } +template +CUSPATIAL_HOST_DEVICE bool multipoint_range::is_single_point_range() +{ + return num_multipoints() == num_points(); +} + } // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh index 1ee08b5cc..1e2ac7c42 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh @@ -141,6 +141,11 @@ class multipoint_range { template CUSPATIAL_HOST_DEVICE auto point(IndexType idx); + /** + * @brief Returns `true` if the range contains only single points + */ + CUSPATIAL_HOST_DEVICE bool is_single_point_range(); + protected: /// Iterator to the start of the index array of start positions to each multipoint. GeometryIterator _geometry_begin; From b725b52f70067c3171ee67081de9f9b0af6b5280 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 14:51:02 -0800 Subject: [PATCH 19/34] docs, type checks in range ctor --- .../experimental/detail/point_polygon_distance.cuh | 9 --------- .../detail/ranges/multipoint_range.cuh | 2 ++ .../detail/ranges/multipolygon_range.cuh | 2 ++ .../experimental/point_polygon_distance.cuh | 14 +++++++++----- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index b217c5f35..b39ba6dec 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -16,8 +16,6 @@ #pragma once -#include - #include #include #include @@ -114,13 +112,6 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, { using T = typename MultiPointRange::element_t; - static_assert(is_same_floating_point(), - "Inputs must have same floating point value type."); - - static_assert( - is_same, typename MultiPointRange::point_t, typename MultiPolygonRange::point_t>(), - "Inputs must be cuspatial::vec_2d"); - CUSPATIAL_EXPECTS(multipoints.size() == multipolygons.size(), "Must have the same number of input rows."); diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh index 7b7ce5243..f28ed5ef2 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh @@ -63,6 +63,8 @@ multipoint_range::multipoint_range(GeometryIterat _points_begin(points_begin), _points_end(points_end) { + static_assert(is_vec_2d>(), + "Coordinate range must be constructed with iterators to vec_2d."); } template diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index 0228d1de8..7e0e868bf 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -103,6 +103,8 @@ multipolygon_range::m _point_begin(point_begin), _point_end(point_end) { + static_assert(is_vec_2d>(), + "Coordinate range must be constructed with iterators to vec_2d."); } template OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, From ab59e7d1352ed3b3c843abbfc95f09b06e8ad57f Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 15:13:15 -0800 Subject: [PATCH 20/34] use range based for loop in is_point_in_polygon --- .../experimental/detail/algorithm/is_point_in_polygon.cuh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh index 135011122..5ddf208d7 100644 --- a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -16,6 +16,7 @@ #pragma once +#include "cuspatial/experimental/geometry_collection/multipoint_ref.cuh" #include #include @@ -48,8 +49,8 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR auto b = last_segment.v2; bool y0_flag = b.y > test_point.y; bool y1_flag; - for (auto it = ring.point_begin(); it != ring.point_end(); ++it) { - vec_2d a = *it; + auto ring_points = multipoint_ref{ring.point_begin(), ring.point_end()}; + for (vec_2d 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; From b136c0b3f1707e870af3c4fae2f2dbfebc3f89a5 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 17:49:05 -0800 Subject: [PATCH 21/34] Apply suggestions from code review --- .../experimental/detail/algorithm/is_point_in_polygon.cuh | 3 +-- .../cuspatial/experimental/detail/geometry/polygon_ref.cuh | 2 +- .../cuspatial/experimental/detail/point_polygon_distance.cuh | 1 + .../experimental/detail/ranges/multipolygon_range.cuh | 3 +-- .../cuspatial/experimental/ranges/multipolygon_range.cuh | 5 +++-- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh index 5ddf208d7..e9f43818d 100644 --- a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -16,7 +16,7 @@ #pragma once -#include "cuspatial/experimental/geometry_collection/multipoint_ref.cuh" +#include #include #include @@ -67,7 +67,6 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR is_colinear = float_equal(run * rise_to_point, run_to_point * rise); if (is_colinear) { break; } - // y0_flag = a.y > test_point.y; y1_flag = a.y > test_point.y; if (y1_flag != y0_flag) { // Transform the following inequality to avoid division diff --git a/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh index 02a917cd4..3c1dfeb68 100644 --- a/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/detail/geometry/polygon_ref.cuh @@ -1,5 +1,5 @@ /* - * Copyright (c) 2022, NVIDIA CORPORATION. + * 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. diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index b39ba6dec..110a16087 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -87,6 +87,7 @@ void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoin if (geometry_idx == MultiPolygonRange::INVALID_INDEX) continue; if (intersects[geometry_idx]) { + // Leading thread of the pair writes to the output if (multipolygons.is_first_point_of_multipolygon(idx, geometry_idx)) distances[geometry_idx] = T{0.0}; continue; diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh index 7e0e868bf..452d1be2e 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipolygon_range.cuh @@ -1,5 +1,5 @@ /* - * Copyright (c) 2022, NVIDIA CORPORATION. + * 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. @@ -224,7 +224,6 @@ multipolygon_range:: return multipolygon_range:: INVALID_INDEX; - auto part_idx = part_idx_from_ring_idx(ring_idx); return geometry_idx_from_part_idx(part_idx_from_ring_idx(ring_idx)); } diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index 775c7dba4..7850a48db 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -1,5 +1,5 @@ /* - * Copyright (c) 2022, NVIDIA CORPORATION. + * 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. @@ -100,6 +100,7 @@ class multipolygon_range { /// Return the iterator to the one past the last multipolygon in the range. CUSPATIAL_HOST_DEVICE auto end() { return multipolygon_end(); } + /// Given the index of a segment, return the geometry (multipolygon) index where the /// segment locates. /// Segment index is the index to the starting point of the segment. If the @@ -116,7 +117,7 @@ class multipolygon_range { template CUSPATIAL_HOST_DEVICE auto get_segment(IndexType segment_idx); - /// Returns `true` if `point_idx`th point is the first point of its + /// Returns `true` if `point_idx`th point is the first point of `geometry_idx`th /// multipolygon template CUSPATIAL_HOST_DEVICE bool is_first_point_of_multipolygon(IndexType1 point_idx, From 744f32f2e0febf59131ad274bd1282a3998a7b3f Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 17:52:43 -0800 Subject: [PATCH 22/34] add docs --- .../cuspatial/experimental/detail/point_polygon_distance.cuh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index 110a16087..a3395de97 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -43,6 +43,9 @@ namespace cuspatial { namespace detail { +/** + * @brief For each point in the multipoint, compute point-in-multipolygon in corresponding pair. + */ template struct point_in_multipolygon_test_functor { MultiPointRange multipoints; From bb6c63785a12cb42fd9124f1e9cbce64b97e9341 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 8 Mar 2023 17:53:19 -0800 Subject: [PATCH 23/34] style --- .../cuspatial/experimental/ranges/multipolygon_range.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index 7850a48db..b9cf666ea 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -100,7 +100,7 @@ class multipolygon_range { /// Return the iterator to the one past the last multipolygon in the range. CUSPATIAL_HOST_DEVICE auto end() { return multipolygon_end(); } - + /// Given the index of a segment, return the geometry (multipolygon) index where the /// segment locates. /// Segment index is the index to the starting point of the segment. If the From 8319e7c0f387ee86c16a5968d64d25976a0f2e02 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Fri, 10 Mar 2023 09:55:35 -0800 Subject: [PATCH 24/34] fix bug in PiP tests --- .../detail/algorithm/is_point_in_polygon.cuh | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh index e9f43818d..c59322803 100644 --- a/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/experimental/detail/algorithm/is_point_in_polygon.cuh @@ -33,6 +33,13 @@ namespace detail { * See "Crossings test" section of http://erich.realtimerendering.com/ptinpoly/ * The improvement in addenda is also addopted 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 @@ -104,10 +111,10 @@ __device__ inline bool is_point_in_polygon(Cart2d const& test_point, Cart2dIt poly_points_first, Cart2dItDiffType const& num_poly_points) { - auto polygon = polygon_ref{ring_offsets_first + poly_begin, - ring_offsets_first + poly_end, + auto polygon = polygon_ref{thrust::next(ring_offsets_first, poly_begin), + thrust::next(ring_offsets_first, poly_end + 1), poly_points_first, - poly_points_first + num_poly_points}; + thrust::next(poly_points_first, num_poly_points)}; return is_point_in_polygon(test_point, polygon); } From 86d36d9b36fa936d192f772ea88674f81812cf3c Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 10:22:55 -0700 Subject: [PATCH 25/34] updates with test docs and API docs, address review --- .../geometry_collection/multipolygon_ref.cuh | 1 + .../experimental/point_polygon_distance.cuh | 1 + .../cuspatial_test/vector_factories.cuh | 4 +- .../spatial/point_polygon_distance_test.cu | 60 ++++++++++++++++++- 4 files changed, 62 insertions(+), 4 deletions(-) diff --git a/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh index 64da05797..9affb3fba 100644 --- a/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry_collection/multipolygon_ref.cuh @@ -24,6 +24,7 @@ namespace cuspatial { * @brief Represent a reference to a multipolygon stored in a structure of arrays. * * @tparam PartIterator type of iterator to the part offset array. + * @tparam RingIterator type of iterator to the ring offset array. * @tparam VecIterator type of iterator to the underlying point array. */ template diff --git a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh index a1a0093d9..6398531ef 100644 --- a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh @@ -27,6 +27,7 @@ namespace cuspatial { * @tparam MultiPointRange An instance of template type `MultiPointRange` * @tparam MultiPolygonRange An instance of template type `MultiPolygonRange` * @tparam OutputIt iterator type for output array. Must meet the requirements of [LRAI](LinkLRAI). + * Must be an iterator to type convertible from floating points. * * @param multipoints Range of multipoints, one per computed distance pair. * @param multipolygons Range of multilinestrings, one per computed distance pair. diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh index 7747dfe33..9053ed2c0 100644 --- a/cpp/include/cuspatial_test/vector_factories.cuh +++ b/cpp/include/cuspatial_test/vector_factories.cuh @@ -73,10 +73,10 @@ class multipolygon_array { { } - /// Return the number of multilinestrings + /// Return the number of multipolygons auto size() { return _geometry_offsets_array.size() - 1; } - /// Return range object of the multilinestring array + /// Return range object of the multipolygon array auto range() { return multipolygon_range(_geometry_offsets_array.begin(), diff --git a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu index ce5f78b23..64219b55f 100644 --- a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu +++ b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu @@ -62,6 +62,7 @@ using TestTypes = ::testing::Types; TYPED_TEST_CASE(PairwisePointPolygonDistanceTest, TestTypes); +// Inputs are empty columns TYPED_TEST(PairwisePointPolygonDistanceTest, ZeroPairs) { using T = TypeParam; @@ -76,6 +77,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, ZeroPairs) std::initializer_list{}); } +// Point in 1 ring polygon. +// POINT (0 0) +// POLYGON ((-1 -1, 1, -1, 1 1, -1 1, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing) { using T = TypeParam; @@ -90,6 +94,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing) {0.0}); } +// Point outside 1 ring polygon. +// POINT (0 2) +// POLYGON ((-1 -1, 1 -1, 1 1, -1 1, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing2) { using T = TypeParam; @@ -104,6 +111,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonOneRing2) {1.0}); } +// Point in the hole. Polygon has two rings. Point in the hole. +// POINT (0 0) +// POLYGON ((-2 -2, 2 -2, 2 2, -2 2, -2 -2), (-1 -1, 1 -1, 1 1, -1 1, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings) { using T = TypeParam; @@ -129,6 +139,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings) {1.0}); } +// Point in polygon. Polygon has two rings. Point outside of polygon. +// POINT (1.5 0) +// POLYGON ((-2 -2, 2 -2, 2 2, -2 2, -2 -2), (-1 -1, 1 -1, 1 1, -1 1, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings2) { using T = TypeParam; @@ -154,6 +167,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings2) {0.0}); } +// Point outside of polygon. Polygon has two rings. Point outside of polygon. +// POINT (3 0) +// POLYGON ((-2 -2, 2 -2, 2 2, -2 2, -2 -2), (-1 -1, 1 -1, 1 1, -1 1, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings3) { using T = TypeParam; @@ -179,6 +195,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairOnePolygonTwoRings3) {1.0}); } +// 1 Multipolygon with 2 Polygons. Point intersects with second polygon +// POINT (1 1) +// MULTIPOLYGON (((-2 -2, 0 -2, 0 0, -2 0, -2 -2)), ((0 0, 2 0, 2 2, 0 2, 0 0))) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing) { using T = TypeParam; @@ -204,6 +223,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing) {0.0}); } +// 1 Multipolygon with 2 Polygons. Point intersects with first polygon. +// POINT (-1 -1) +// MULTIPOLYGON (((-2 -2, 0 -2, 0 0, -2 0, -2 -2)), ((0 0, 2 0, 2 2, 0 2, 0 0))) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing2) { using T = TypeParam; @@ -229,6 +251,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing2) {0.0}); } +// 1 Multipolygon with 2 Polygons. Point does not intersect. Closer to first polygon. +// POINT (-1 0.5) +// MULTIPOLYGON (((-2 -2, 0 -2, 0 0, -2 0, -2 -2)), ((0 0, 2 0, 2 2, 0 2, 0 0))) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing3) { using T = TypeParam; @@ -254,6 +279,9 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing3) {0.5}); } +// 1 Multipolygon with 2 Polygons. Point does not intersect. Closer to second polygon. +// POINT (-0.3, 1) +// MULTIPOLYGON (((-2 -2, 0 -2, 0 0, -2 0, -2 -2)), ((0 0, 2 0, 2 2, 0 2, 0 0))) TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing4) { using T = TypeParam; @@ -279,6 +307,12 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairTwoPolygonOneRing4) {0.3}); } +// Two Pairs. +// POINT (-0.6 -0.6) +// POLYGON ((-1 -1, 0 0, 0 1, -1 -1)) +// +// POINT (0 0) +// POLYGON ((1 1, 1 0, 2 2, 1 1)) TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairOnePolygonOneRing) { using T = TypeParam; @@ -302,6 +336,12 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairOnePolygonOneRing) {0.0, 1.0}); } +// Two Pairs, each polygon has two rings. +// POINT (2.5, 3) +// POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1)) +// +// POINT (-1.75, -1.5) +// POLYGON ((0 0, -3 0, -3 -3, 0 0), (-1 -1, -2 -1, -2 -2, -1 -1)) TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairTwoPolygonTwoRing) { using T = TypeParam; @@ -336,6 +376,19 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairTwoPolygonTwoRing) {0.0, 0.17677669529663687}); } +// Three Polygons +// POINT (1 1) +// POLYGON ((0 1, -1 -1, 1 -1, 0 1), (0 0.5, 0.5 -0.5, -0.5 -0.5, 0 0.5)) +// +// POINT (2 2) +// POLYGON ((1 1, 1 2, 2 1, 1 1)) +// +// POINT (1.5 0) +// POLYGON ( +// (-3 -3, 3 -3, 3 3, -3 3, -3 -3), +// (-2 -2, -1 -2, -1 2, -2 2, -2 -2), +// (2 2, 2 -2, 1 -2, 1 2, 2 2) +// ) TYPED_TEST(PairwisePointPolygonDistanceTest, ThreePolygons) { using T = TypeParam; @@ -383,6 +436,7 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, ThreePolygons) {0.894427190999916, 0.7071067811865476, 0.5}); } +// Multipoint tests: 1 multipoint - 1 polygon. No Intersection. TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon) { using T = TypeParam; @@ -397,6 +451,7 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon) {1.3416407864998738}); } +// Multipoint tests: 1 multipoint - 1 polygon. Intesects. TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon2) { using T = TypeParam; @@ -411,6 +466,7 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, OnePairMultiPointOnePolygon2) {0.0}); } +// Multipoint tests: 2 multipoints - 2 polygons. TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairMultiPointOnePolygon2) { using T = TypeParam; @@ -418,10 +474,10 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairMultiPointOnePolygon2) CUSPATIAL_RUN_TEST( this->run_single, - {{P{0, 2}, P{0, 0}}, {P{1, 1}, P{-1, -1}}}, + {{P{0, 2}, P{3, 0}}, {P{1, 1}, P{-1, -1}}}, {0, 1, 2}, {0, 1, 2}, {0, 5, 9}, {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}, P{-1, 1}, P{1, 1}, P{0, -1}, P{-1, 1}}, - {0.0, 0.0}); + {1.0, 0.0}); } From 7bf8cbfa634ccfc1e83fcf76b79842cdfe01a13c Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 10:48:41 -0700 Subject: [PATCH 26/34] update offsets_to_keys --- .../detail/utility/offset_to_keys.cuh | 51 ------------------- ...inestring_intersection_with_duplicates.cuh | 36 ++----------- .../detail/point_polygon_distance.cuh | 4 +- .../spatial/linestring_intersection_test.cu | 5 +- 4 files changed, 8 insertions(+), 88 deletions(-) delete mode 100644 cpp/include/cuspatial/detail/utility/offset_to_keys.cuh diff --git a/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh b/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh deleted file mode 100644 index 70665b4db..000000000 --- a/cpp/include/cuspatial/detail/utility/offset_to_keys.cuh +++ /dev/null @@ -1,51 +0,0 @@ -/* - * 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 -#include - -/** @brief Given list offset and row `i`, return a unique key that represent the list of `i`. - * - * The key is computed by performing a `upper_bound` search with `i` in the offset array. - * Then subtracts the position with the start of offset array. - * - * Example: - * offset: 0 0 0 1 3 4 4 4 - * i: 0 1 2 3 - * key: 3 4 4 5 - * - * Note that the values of `key`, {offset[3], offset[4], offset[5]} denotes the ending - * position of the first 3 non-empty list. - */ -template -struct offsets_to_keys_functor { - Iterator _offsets_begin; - Iterator _offsets_end; - - offsets_to_keys_functor(Iterator offset_begin, Iterator offset_end) - : _offsets_begin(offset_begin), _offsets_end(offset_end) - { - } - - template - IndexType __device__ operator()(IndexType i) - { - return thrust::distance(_offsets_begin, - thrust::upper_bound(thrust::seq, _offsets_begin, _offsets_end, i)); - } -}; diff --git a/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh b/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh index dc3f4c5fd..a01ba0a63 100644 --- a/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh +++ b/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -123,37 +124,6 @@ struct offsets_update_functor { } }; -/** @brief Given list offset and row `i`, return a unique key that represent the list of `i`. - * - * The key is computed by performing a `upper_bound` search with `i` in the offset array. - * Then subtracts the position with the start of offset array. - * - * Example: - * offset: 0 0 0 1 3 4 4 4 - * i: 0 1 2 3 - * key: 3 4 4 5 - * - * Note that the values of `key`, {offset[3], offset[4], offset[5]} denotes the ending - * position of the first 3 non-empty list. - */ -template -struct offsets_to_keys_functor { - Iterator _offsets_begin; - Iterator _offsets_end; - - offsets_to_keys_functor(Iterator offset_begin, Iterator offset_end) - : _offsets_begin(offset_begin), _offsets_end(offset_end) - { - } - - template - IndexType __device__ operator()(IndexType i) - { - return thrust::distance(_offsets_begin, - thrust::upper_bound(thrust::seq, _offsets_begin, _offsets_end, i)); - } -}; - } // namespace intersection_functors /// Internal structure to provide convenient access to the intersection intermediate id arrays. @@ -318,7 +288,7 @@ struct linestring_intersection_intermediates { rmm::device_uvector reduced_keys(num_pairs(), stream); rmm::device_uvector reduced_flags(num_pairs(), stream); auto keys_begin = make_counting_transform_iterator( - 0, intersection_functors::offsets_to_keys_functor{offsets->begin(), offsets->end()}); + 0, upper_bound_index_functor{offsets->begin(), offsets->end()}); auto [keys_end, flags_end] = thrust::reduce_by_key(rmm::exec_policy(stream), @@ -387,7 +357,7 @@ struct linestring_intersection_intermediates { auto keys_begin() { return make_counting_transform_iterator( - 0, intersection_functors::offsets_to_keys_functor{offsets->begin(), offsets->end()}); + 0, upper_bound_index_functor{offsets->begin(), offsets->end()}); } /// Return the number of pairs in the intermediates diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index a3395de97..7207e3c19 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -20,7 +20,7 @@ #include #include #include -#include +#include #include #include #include @@ -139,7 +139,7 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, detail::zero_data_async(multipoint_intersects.begin(), multipoint_intersects.end(), stream); auto offset_as_key_it = detail::make_counting_transform_iterator( - 0, offsets_to_keys_functor{multipoints.offsets_begin(), multipoints.offsets_end()}); + 0, detail::upper_bound_index_functor{multipoints.offsets_begin(), multipoints.offsets_end()}); thrust::reduce_by_key(rmm::exec_policy(stream), offset_as_key_it, diff --git a/cpp/tests/experimental/spatial/linestring_intersection_test.cu b/cpp/tests/experimental/spatial/linestring_intersection_test.cu index 762a88c64..02769dc33 100644 --- a/cpp/tests/experimental/spatial/linestring_intersection_test.cu +++ b/cpp/tests/experimental/spatial/linestring_intersection_test.cu @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -83,8 +84,8 @@ linestring_intersection_result segment_sort_intersection_result( rmm::device_uvector geometry_collection_keys(num_geoms, stream); auto geometry_collection_keys_begin = detail::make_counting_transform_iterator( 0, - detail::intersection_functors::offsets_to_keys_functor{ - result.geometry_collection_offset->begin(), result.geometry_collection_offset->end()}); + detail::upper_bound_index_functor{result.geometry_collection_offset->begin(), + result.geometry_collection_offset->end()}); thrust::copy(rmm::exec_policy(stream), geometry_collection_keys_begin, geometry_collection_keys_begin + num_geoms, From ff48dc7d2b9a5b3bbeabe95e8756cc9cfbbf3ee7 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 10:49:30 -0700 Subject: [PATCH 27/34] update floating_point docs --- cpp/include/cuspatial/detail/utility/floating_point.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/include/cuspatial/detail/utility/floating_point.cuh b/cpp/include/cuspatial/detail/utility/floating_point.cuh index 7a44aeb73..362994e8b 100644 --- a/cpp/include/cuspatial/detail/utility/floating_point.cuh +++ b/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. From b4420a73a651e2e9b29e3a73082458607c143be0 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 11:08:49 -0700 Subject: [PATCH 28/34] use any_of --- .../detail/point_polygon_distance.cuh | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index 7207e3c19..e4b71cb76 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -33,6 +34,7 @@ #include #include #include +#include #include #include @@ -48,6 +50,8 @@ namespace detail { */ template struct point_in_multipolygon_test_functor { + using T = typename MultiPointRange::element_t; + MultiPointRange multipoints; MultiPolygonRange multipolygons; @@ -59,13 +63,15 @@ struct point_in_multipolygon_test_functor { template uint8_t __device__ operator()(IndexType pidx) { - auto point = thrust::raw_reference_cast(multipoints.point(pidx)); - auto geometry_idx = multipoints.geometry_idx_from_point_idx(pidx); + vec_2d const& point = multipoints.point(pidx); + auto geometry_idx = multipoints.geometry_idx_from_point_idx(pidx); + + auto const& polys = multipolygons[geometry_idx]; + bool intersects = + thrust::any_of(thrust::seq, polys.begin(), polys.end(), [&point] __device__(auto poly) { + return is_point_in_polygon(point, poly); + }); - bool intersects = false; - for (auto polygon : multipolygons[geometry_idx]) { - intersects = intersects || is_point_in_polygon(point, polygon); - } return static_cast(intersects); } }; From 44e2843be533e117696c773c301fd7bb363a1166 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 13:40:44 -0700 Subject: [PATCH 29/34] add large (multigrid) tests, address reviews --- .../detail/point_polygon_distance.cuh | 8 +-- .../detail/ranges/multilinestring_range.cuh | 2 +- .../detail/ranges/multipoint_range.cuh | 2 +- .../experimental/geometry/linestring_ref.cuh | 2 +- .../experimental/geometry/segment.cuh | 2 +- .../experimental/point_polygon_distance.cuh | 5 +- .../experimental/ranges/multipoint_range.cuh | 2 - .../ranges/multipolygon_range.cuh | 9 ++-- .../cuspatial_test/vector_factories.cuh | 19 ++++--- .../spatial/point_polygon_distance_test.cu | 50 +++++++++++++++++-- 10 files changed, 74 insertions(+), 27 deletions(-) diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index e4b71cb76..d5d12d3a2 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -67,6 +67,7 @@ struct point_in_multipolygon_test_functor { auto geometry_idx = multipoints.geometry_idx_from_point_idx(pidx); auto const& polys = multipolygons[geometry_idx]; + // TODO: benchmark against range based for loop bool intersects = thrust::any_of(thrust::seq, polys.begin(), polys.end(), [&point] __device__(auto poly) { return is_point_in_polygon(point, poly); @@ -96,14 +97,13 @@ void __global__ pairwise_point_polygon_distance_kernel(MultiPointRange multipoin if (geometry_idx == MultiPolygonRange::INVALID_INDEX) continue; if (intersects[geometry_idx]) { - // Leading thread of the pair writes to the output - if (multipolygons.is_first_point_of_multipolygon(idx, geometry_idx)) - distances[geometry_idx] = T{0.0}; + distances[geometry_idx] = T{0.0}; continue; } + auto [a, b] = multipolygons.get_segment(idx); + T dist_squared = std::numeric_limits::max(); - auto [a, b] = multipolygons.get_segment(idx); for (vec_2d point : multipoints[geometry_idx]) { dist_squared = min(dist_squared, point_to_segment_distance_squared(point, a, b)); } diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh index 851269313..5add3a769 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multilinestring_range.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. diff --git a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh index f28ed5ef2..fa7d4002b 100644 --- a/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/detail/ranges/multipoint_range.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. diff --git a/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh b/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh index 707b9e0cc..4d0092420 100644 --- a/cpp/include/cuspatial/experimental/geometry/linestring_ref.cuh +++ b/cpp/include/cuspatial/experimental/geometry/linestring_ref.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. diff --git a/cpp/include/cuspatial/experimental/geometry/segment.cuh b/cpp/include/cuspatial/experimental/geometry/segment.cuh index 8667e174c..214819694 100644 --- a/cpp/include/cuspatial/experimental/geometry/segment.cuh +++ b/cpp/include/cuspatial/experimental/geometry/segment.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. diff --git a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh index 6398531ef..b2538a0ae 100644 --- a/cpp/include/cuspatial/experimental/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/point_polygon_distance.cuh @@ -24,13 +24,14 @@ namespace cuspatial { * @ingroup distance * @brief Computes pairwise multipoint to multipolygon distance * - * @tparam MultiPointRange An instance of template type `MultiPointRange` - * @tparam MultiPolygonRange An instance of template type `MultiPolygonRange` + * @tparam MultiPointRange An instance of template type `multipoint_range` + * @tparam MultiPolygonRange An instance of template type `multipolygon_range` * @tparam OutputIt iterator type for output array. Must meet the requirements of [LRAI](LinkLRAI). * Must be an iterator to type convertible from floating points. * * @param multipoints Range of multipoints, one per computed distance pair. * @param multipolygons Range of multilinestrings, one per computed distance pair. + * @param stream The CUDA stream on which to perform computations * @return Output Iterator past the last distance computed * * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator diff --git a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh index 1e2ac7c42..5cc06f62b 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipoint_range.cuh @@ -51,8 +51,6 @@ class multipoint_range { using point_t = iterator_value_type; using element_t = iterator_vec_base_type; - int32_t INVALID_IDX = -1; - /** * @brief Construct a new multipoint array object */ diff --git a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh index b9cf666ea..2e05066ed 100644 --- a/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh +++ b/cpp/include/cuspatial/experimental/ranges/multipolygon_range.cuh @@ -101,11 +101,10 @@ class multipolygon_range { /// Return the iterator to the one past the last multipolygon in the range. CUSPATIAL_HOST_DEVICE auto end() { return multipolygon_end(); } - /// Given the index of a segment, return the geometry (multipolygon) index where the - /// segment locates. - /// Segment index is the index to the starting point of the segment. If the - /// index is the last point of the ring, then it is not a valid index. - /// This function returns multipolygon_range::INVALID_INDEX if the index is invalid. + /// Given the index of a segment, return the index of the geometry (multipolygon) that contains + /// the segment. Segment index is the index to the starting point of the segment. If the index is + /// the last point of the ring, then it is not a valid index. This function returns + /// multipolygon_range::INVALID_INDEX if the index is invalid. template CUSPATIAL_HOST_DEVICE auto geometry_idx_from_segment_idx(IndexType segment_idx); diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh index 9053ed2c0..d4887e6e4 100644 --- a/cpp/include/cuspatial_test/vector_factories.cuh +++ b/cpp/include/cuspatial_test/vector_factories.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. @@ -30,6 +30,13 @@ namespace cuspatial { namespace test { +template +auto make_device_vector(Range rng) +{ + using T = typename Range::value_type; + return rmm::device_vector(rng.begin(), rng.end()); +} + template auto make_device_vector(std::initializer_list inl) { @@ -96,11 +103,11 @@ class multipolygon_array { CoordinateArray _coordinate_offsets_array; }; -template -auto make_multipolygon_array(std::initializer_list geometry_inl, - std::initializer_list part_inl, - std::initializer_list ring_inl, - std::initializer_list> coord_inl) +template +auto make_multipolygon_array(IndexRange geometry_inl, + IndexRange part_inl, + IndexRange ring_inl, + CoordRange coord_inl) { return multipolygon_array{make_device_vector(geometry_inl), make_device_vector(part_inl), diff --git a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu index 64219b55f..bec557c4f 100644 --- a/cpp/tests/experimental/spatial/point_polygon_distance_test.cu +++ b/cpp/tests/experimental/spatial/point_polygon_distance_test.cu @@ -17,7 +17,9 @@ #include #include +#include #include +#include #include #include @@ -29,6 +31,8 @@ using namespace cuspatial; using namespace cuspatial::test; +double constexpr PI = 3.14159265358979323846; + template struct PairwisePointPolygonDistanceTest : public ::testing::Test { rmm::cuda_stream_view stream() { return rmm::cuda_stream_default; } @@ -40,12 +44,29 @@ struct PairwisePointPolygonDistanceTest : public ::testing::Test { std::initializer_list multipolygon_ring_offsets, std::initializer_list> multipolygon_coordinates, std::initializer_list expected) + { + std::vector> multipolygon_coordinates_vec(multipolygon_coordinates); + return this->run_single(multipoints, + multipolygon_geometry_offsets, + multipolygon_part_offsets, + multipolygon_ring_offsets, + multipolygon_coordinates_vec, + expected); + } + + void run_single(std::initializer_list>> multipoints, + std::initializer_list multipolygon_geometry_offsets, + std::initializer_list multipolygon_part_offsets, + std::initializer_list multipolygon_ring_offsets, + std::vector> const& multipolygon_coordinates, + std::initializer_list expected) { auto d_multipoints = make_multipoints_array(multipoints); - auto d_multipolygons = make_multipolygon_array(multipolygon_geometry_offsets, - multipolygon_part_offsets, - multipolygon_ring_offsets, - multipolygon_coordinates); + auto d_multipolygons = make_multipolygon_array( + range{multipolygon_geometry_offsets.begin(), multipolygon_geometry_offsets.end()}, + range{multipolygon_part_offsets.begin(), multipolygon_part_offsets.end()}, + range{multipolygon_ring_offsets.begin(), multipolygon_ring_offsets.end()}, + range{multipolygon_coordinates.begin(), multipolygon_coordinates.end()}); auto got = rmm::device_uvector(d_multipoints.size(), stream()); @@ -481,3 +502,24 @@ TYPED_TEST(PairwisePointPolygonDistanceTest, TwoPairMultiPointOnePolygon2) {P{-1, -1}, P{1, -1}, P{1, 1}, P{-1, 1}, P{-1, -1}, P{-1, 1}, P{1, 1}, P{0, -1}, P{-1, 1}}, {1.0, 0.0}); } + +// Large distance test +TYPED_TEST(PairwisePointPolygonDistanceTest, DistanceTestManyVertex) +{ + using T = TypeParam; + using P = vec_2d; + + std::size_t num_vertex = 2000; + P centroid{0.0, 0.0}; + T radius = 1.0; + + std::vector

polygon; + auto it = detail::make_counting_transform_iterator(0, [](auto i) { + T theta = i / (2 * PI); + return P{cos(theta), sin(theta)}; + }); + std::copy(it, it + num_vertex, std::back_inserter(polygon)); + + CUSPATIAL_RUN_TEST( + this->run_single, {{P{0.0, 0.0}}}, {0, 1}, {0, 1}, {0, num_vertex + 1}, polygon, {0.0}); +} From 49972d4d213fa295387cfe482b78392aacb3d2d5 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 18:04:18 -0700 Subject: [PATCH 30/34] add upper_bound_index.cuh --- .../detail/utility/upper_bound_index.cuh | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 cpp/include/cuspatial/detail/utility/upper_bound_index.cuh diff --git a/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh b/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh new file mode 100644 index 000000000..048b76edf --- /dev/null +++ b/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh @@ -0,0 +1,58 @@ +/* + * 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 +#include + +namespace cuspatial { +namespace detail { + +/** @brief Returns the index of the first element in the range `[first, last)` such that + * `value < element` is true (i.e. strictly greater), or `distance(first, last)` if no such element + * is found. + * + * A common use of the functor is to apply `*_by_key` algorithm to list values. The keys + * iterator can be constructed with `upper_bound_index_functor{offsets.begin(), offsets.end()}`. + * + * Example: + * ``` + * offset: 0 0 0 1 3 4 4 4 + * i: 0 1 2 3 + * key: 3 4 4 5 + * ``` + */ +template +struct upper_bound_index_functor { + Iterator _offsets_begin; + Iterator _offsets_end; + + upper_bound_index_functor(Iterator offset_begin, Iterator offset_end) + : _offsets_begin(offset_begin), _offsets_end(offset_end) + { + } + + template + IndexType __device__ operator()(IndexType i) + { + return thrust::distance(_offsets_begin, + thrust::upper_bound(thrust::seq, _offsets_begin, _offsets_end, i)); + } +}; + +} // namespace detail +} // namespace cuspatial From f000e00585be9bf00c0f831160c94d49d13ca379 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 14 Mar 2023 22:52:16 -0700 Subject: [PATCH 31/34] remove upper_bound_index --- .../detail/utility/upper_bound_index.cuh | 58 ------------------- ...inestring_intersection_with_duplicates.cuh | 11 +--- .../detail/point_polygon_distance.cuh | 9 +-- .../spatial/linestring_intersection_test.cu | 8 +-- 4 files changed, 11 insertions(+), 75 deletions(-) delete mode 100644 cpp/include/cuspatial/detail/utility/upper_bound_index.cuh diff --git a/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh b/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh deleted file mode 100644 index 048b76edf..000000000 --- a/cpp/include/cuspatial/detail/utility/upper_bound_index.cuh +++ /dev/null @@ -1,58 +0,0 @@ -/* - * 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 -#include - -namespace cuspatial { -namespace detail { - -/** @brief Returns the index of the first element in the range `[first, last)` such that - * `value < element` is true (i.e. strictly greater), or `distance(first, last)` if no such element - * is found. - * - * A common use of the functor is to apply `*_by_key` algorithm to list values. The keys - * iterator can be constructed with `upper_bound_index_functor{offsets.begin(), offsets.end()}`. - * - * Example: - * ``` - * offset: 0 0 0 1 3 4 4 4 - * i: 0 1 2 3 - * key: 3 4 4 5 - * ``` - */ -template -struct upper_bound_index_functor { - Iterator _offsets_begin; - Iterator _offsets_end; - - upper_bound_index_functor(Iterator offset_begin, Iterator offset_end) - : _offsets_begin(offset_begin), _offsets_end(offset_end) - { - } - - template - IndexType __device__ operator()(IndexType i) - { - return thrust::distance(_offsets_begin, - thrust::upper_bound(thrust::seq, _offsets_begin, _offsets_end, i)); - } -}; - -} // namespace detail -} // namespace cuspatial diff --git a/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh b/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh index f917927c5..c56c4bfbe 100644 --- a/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh +++ b/cpp/include/cuspatial/experimental/detail/linestring_intersection_with_duplicates.cuh @@ -16,11 +16,11 @@ #include #include -#include #include #include #include #include +#include #include #include @@ -287,8 +287,7 @@ struct linestring_intersection_intermediates { // Use `reduce_by_key` to compute the number of removed geometry per list. rmm::device_uvector reduced_keys(num_pairs(), stream); rmm::device_uvector reduced_flags(num_pairs(), stream); - auto keys_begin = make_counting_transform_iterator( - 0, upper_bound_index_functor{offsets->begin(), offsets->end()}); + auto keys_begin = make_geometry_id_iterator(offsets->begin(), offsets->end()); auto [keys_end, flags_end] = thrust::reduce_by_key(rmm::exec_policy(stream), @@ -354,11 +353,7 @@ struct linestring_intersection_intermediates { } /// Return list-id corresponding to the geometry - auto keys_begin() - { - return make_counting_transform_iterator( - 0, upper_bound_index_functor{offsets->begin(), offsets->end()}); - } + auto keys_begin() { return make_geometry_id_iterator(offsets->begin(), offsets->end()); } /// Return the number of pairs in the intermediates auto num_pairs() { return offsets->size() - 1; } diff --git a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh index d5d12d3a2..704f67a6b 100644 --- a/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh +++ b/cpp/include/cuspatial/experimental/detail/point_polygon_distance.cuh @@ -20,10 +20,10 @@ #include #include #include -#include #include #include #include +#include #include #include @@ -120,7 +120,8 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, OutputIt distances_first, rmm::cuda_stream_view stream) { - using T = typename MultiPointRange::element_t; + using T = typename MultiPointRange::element_t; + using index_t = typename MultiPointRange::index_t; CUSPATIAL_EXPECTS(multipoints.size() == multipolygons.size(), "Must have the same number of input rows."); @@ -144,8 +145,8 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, rmm::device_uvector multipoint_intersects(multipoints.num_multipoints(), stream); detail::zero_data_async(multipoint_intersects.begin(), multipoint_intersects.end(), stream); - auto offset_as_key_it = detail::make_counting_transform_iterator( - 0, detail::upper_bound_index_functor{multipoints.offsets_begin(), multipoints.offsets_end()}); + auto offset_as_key_it = + make_geometry_id_iterator(multipoints.offsets_begin(), multipoints.offsets_end()); thrust::reduce_by_key(rmm::exec_policy(stream), offset_as_key_it, diff --git a/cpp/tests/experimental/spatial/linestring_intersection_test.cu b/cpp/tests/experimental/spatial/linestring_intersection_test.cu index a2956fe76..0ec93957f 100644 --- a/cpp/tests/experimental/spatial/linestring_intersection_test.cu +++ b/cpp/tests/experimental/spatial/linestring_intersection_test.cu @@ -20,7 +20,6 @@ #include #include -#include #include #include #include @@ -82,10 +81,9 @@ linestring_intersection_result segment_sort_intersection_result( // Compute keys for each row in the union column. Rows of the same list // are assigned the same label. rmm::device_uvector geometry_collection_keys(num_geoms, stream); - auto geometry_collection_keys_begin = detail::make_counting_transform_iterator( - 0, - detail::upper_bound_index_functor{result.geometry_collection_offset->begin(), - result.geometry_collection_offset->end()}); + auto geometry_collection_keys_begin = make_geometry_id_iterator( + result.geometry_collection_offset->begin(), result.geometry_collection_offset->end()); + thrust::copy(rmm::exec_policy(stream), geometry_collection_keys_begin, geometry_collection_keys_begin + num_geoms, From 57504e6e4226b7b66d46dd60647f48f65329680d Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 15 Mar 2023 08:39:26 -0700 Subject: [PATCH 32/34] update geometry id factory --- cpp/include/cuspatial/experimental/iterator_factory.cuh | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cpp/include/cuspatial/experimental/iterator_factory.cuh b/cpp/include/cuspatial/experimental/iterator_factory.cuh index 4f1280e9a..67c35e3fc 100644 --- a/cpp/include/cuspatial/experimental/iterator_factory.cuh +++ b/cpp/include/cuspatial/experimental/iterator_factory.cuh @@ -164,9 +164,8 @@ struct index_to_geometry_id { CUSPATIAL_HOST_DEVICE auto operator()(IndexT idx) { - return thrust::distance( - geometry_begin, - thrust::prev(thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx))); + return thrust::distance(geometry_begin, + thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx)); } }; From 48659221fa7d8ff5eb3d68133d63cd53802a7b55 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 15 Mar 2023 09:21:52 -0700 Subject: [PATCH 33/34] add get gtest scripts in cmake --- cpp/CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 8bd8958bc..f9bf49ad9 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -103,7 +103,10 @@ include(cmake/Modules/ConfigureCUDA.cmake) rapids_cpm_init() # find or add cuDF include(cmake/thirdparty/CUSPATIAL_GetCUDF.cmake) - +# find or install GoogleTest +if (CUSPATIAL_BUILD_TESTS) + include(cmake/thirdparty/get_gtests.cmake) +endif() ################################################################################################### # - library targets ------------------------------------------------------------------------------- From cb627fc6982cc7738a998c5a30f93572e2edc9b4 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Wed, 15 Mar 2023 16:52:40 -0700 Subject: [PATCH 34/34] Revert "add get gtest scripts in cmake" This reverts commit 48659221fa7d8ff5eb3d68133d63cd53802a7b55. --- cpp/CMakeLists.txt | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index f9bf49ad9..8bd8958bc 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -103,10 +103,7 @@ include(cmake/Modules/ConfigureCUDA.cmake) rapids_cpm_init() # find or add cuDF include(cmake/thirdparty/CUSPATIAL_GetCUDF.cmake) -# find or install GoogleTest -if (CUSPATIAL_BUILD_TESTS) - include(cmake/thirdparty/get_gtests.cmake) -endif() + ################################################################################################### # - library targets -------------------------------------------------------------------------------