diff --git a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh index 4c988612d..d8a4c0357 100644 --- a/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/algorithm/linestring_distance.cuh @@ -21,8 +21,13 @@ #include +#include #include +#include + +#include + namespace cuspatial { namespace detail { @@ -72,5 +77,84 @@ __global__ void linestring_distance(MultiLinestringRange1 multilinestrings1, } } +template +auto __device__ compute_geometry_id(BoundsIterator bounds_begin, BoundsIterator bounds_end, IndexType idx) +{ + auto it = thrust::prev( + thrust::upper_bound(thrust::seq, bounds_begin, bounds_end, idx)); + auto const geometry_id = thrust::distance(bounds_begin, it); + return thrust::make_pair(it, geometry_id); +} + +/** + * @internal + * @brief The kernel to compute (multi)linestring to (multi)linestring distance + * + * Load balanced kernel to compute distances between one pair of segments from the multilinestring + * and multilinestring. + * + * `intersects` is an optional pointer to a boolean range where the `i`th element indicates the + * `i`th output should be set to 0 and bypass distance computation. This argument is optional, if + * set to nullopt, no distance computation will be bypassed. + */ +template +__global__ void linestring_distance_load_balanced(SegmentRange1 multilinestrings1, + SegmentRange2 multilinestrings2, + IndexRange thread_bounds, + thrust::optional intersects, + OutputIt distances) +{ + using index_t = typename IndexRange::value_type; + + auto num_segment_pairs = thread_bounds[thread_bounds.size() - 1]; + for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_segment_pairs; + idx += gridDim.x * blockDim.x) { + + auto [it, geometry_id] = compute_geometry_id(thread_bounds.begin(), thread_bounds.end(), idx); + + auto first_thread_of_this_block = blockDim.x * blockIdx.x; + auto last_thread_of_block = blockDim.x * (blockIdx.x + 1) - 1; + auto first_thread_of_next_geometry = thread_bounds[geometry_id + 1]; + + bool split_block = first_thread_of_this_block < first_thread_of_next_geometry && first_thread_of_next_geometry <= last_thread_of_block; + + if (intersects.has_value() && intersects.value()[geometry_id]) { + distances[geometry_id] = 0.0f; + continue; + } + + auto const local_idx = idx - *it; + // Retrieve the number of segments in multilinestrings[geometry_id] + auto num_segment_this_multilinestring = + multilinestrings1.multigeometry_count_begin()[geometry_id]; + // The segment id from multilinestring1 this thread is computing (local_id + global_offset) + auto multilinestring1_segment_id = local_idx % num_segment_this_multilinestring + + multilinestrings1.multigeometry_offset_begin()[geometry_id]; + + // The segment id from multilinestring2 this thread is computing (local_id + global_offset) + auto multilinestring2_segment_id = local_idx / num_segment_this_multilinestring + + multilinestrings2.multigeometry_offset_begin()[geometry_id]; + + auto [a, b] = multilinestrings1.begin()[multilinestring1_segment_id]; + auto [c, d] = multilinestrings2.begin()[multilinestring2_segment_id]; + + auto partial = sqrt(squared_segment_distance(a, b, c, d)); + + if (split_block) + atomicMin(&distances[geometry_id], partial); + else + { + // block reduce + typedef cub::BlockReduce BlockReduce; + __shared__ typename BlockReduce::TempStorage temp_storage; + auto aggregate = BlockReduce(temp_storage).Reduce(partial, cub::Min()); + + // atmomic with leading thread + if (cooperative_groups::this_thread_block().thread_rank() == 0) + atomicMin(&distances[geometry_id], aggregate); + } + } +} + } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/distance/distance_utils.cuh b/cpp/include/cuspatial/detail/distance/distance_utils.cuh index 71ee9d5fd..13860f600 100644 --- a/cpp/include/cuspatial/detail/distance/distance_utils.cuh +++ b/cpp/include/cuspatial/detail/distance/distance_utils.cuh @@ -29,6 +29,7 @@ #include #include #include +#include namespace cuspatial { namespace detail { @@ -115,5 +116,45 @@ rmm::device_uvector point_polygon_intersects(MultiPointRange multipoint return multipoint_intersects; } +/** + * @brief Compute the thread bound between two ranges of partitions + * + * @tparam CountIterator1 + * @tparam CountIterator2 + * @param lhs + * @param rhs + * @param stream + * @return rmm::device_uvector + */ +template > +rmm::device_uvector compute_segment_thread_bounds(CountIterator1 lhs_begin, + CountIterator1 lhs_end, + CountIterator2 rhs_begin, + rmm::cuda_stream_view stream) +{ + auto size = thrust::distance(lhs_begin, lhs_end) + 1; + + // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings + // times the number of polygons in a multilinestring-multipolygon pair. + auto segment_count_product_it = + thrust::make_transform_iterator(thrust::make_zip_iterator(lhs_begin, rhs_begin), + thrust::make_zip_function(thrust::multiplies{})); + + // Computes the "thread boundary" of each pair. This array partitions the thread range by + // geometries. E.g. threadIdx within [thread_bounds[i], thread_bounds[i+1]) computes distances of + // the ith pair. + auto thread_bounds = rmm::device_uvector(size, stream); + detail::zero_data_async(thread_bounds.begin(), thread_bounds.end(), stream); + + thrust::inclusive_scan(rmm::exec_policy(stream), + segment_count_product_it, + segment_count_product_it + thread_bounds.size() - 1, + thrust::next(thread_bounds.begin())); + + return thread_bounds; +} + } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh index 19e70b51b..dc901e78b 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_distance.cuh @@ -17,9 +17,11 @@ #pragma once #include +#include #include #include #include +#include #include #include @@ -33,8 +35,8 @@ namespace cuspatial { template -OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, - MultiLinestringRange2 multilinestrings2, +OutputIt pairwise_linestring_distance(MultiLinestringRange1 lhs, + MultiLinestringRange2 rhs, OutputIt distances_first, rmm::cuda_stream_view stream) { @@ -48,25 +50,43 @@ OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, typename MultiLinestringRange2::point_t>(), "All input types must be cuspatial::vec_2d with the same value type"); - CUSPATIAL_EXPECTS(multilinestrings1.size() == multilinestrings2.size(), - "Inputs must have the same number of rows."); + CUSPATIAL_EXPECTS(lhs.size() == rhs.size(), "Inputs must have the same number of rows."); - if (multilinestrings1.size() == 0) return distances_first; + if (lhs.size() == 0) return distances_first; + // Make views to the segments in the multilinestring + auto lhs_segments = lhs._segments(stream); + auto lhs_segments_range = lhs_segments.view(); + + // Make views to the segments in the multilinestring + auto rhs_segments = rhs._segments(stream); + auto rhs_segments_range = rhs_segments.view(); + + auto thread_bounds = + detail::compute_segment_thread_bounds(lhs_segments_range.multigeometry_count_begin(), + lhs_segments_range.multigeometry_count_end(), + rhs_segments_range.multigeometry_count_begin(), + stream); + // Initialize the output range thrust::fill(rmm::exec_policy(stream), distances_first, - distances_first + multilinestrings1.size(), + distances_first + lhs.size(), std::numeric_limits::max()); std::size_t constexpr threads_per_block = 256; - std::size_t const num_blocks = - (multilinestrings1.num_points() + threads_per_block - 1) / threads_per_block; + std::size_t num_threads = 1e8; + std::size_t const num_blocks = (num_threads + threads_per_block - 1) / threads_per_block; - detail::linestring_distance<<>>( - multilinestrings1, multilinestrings2, thrust::nullopt, distances_first); + detail::linestring_distance_load_balanced + <<>>( + lhs_segments_range, + rhs_segments_range, + range(thread_bounds.begin(), thread_bounds.end()), + thrust::nullopt, + distances_first); CUSPATIAL_CUDA_TRY(cudaGetLastError()); - return distances_first + multilinestrings1.size(); + return distances_first + lhs.size(); } } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh index aa0c7cc11..345f6d76e 100644 --- a/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh +++ b/cpp/include/cuspatial/detail/distance/linestring_polygon_distance.cuh @@ -19,8 +19,6 @@ #include "distance_utils.cuh" #include -#include -#include #include #include #include @@ -30,10 +28,10 @@ #include #include -#include -#include +#include #include #include +#include #include #include @@ -45,36 +43,44 @@ namespace cuspatial { namespace detail { +/// Device functor that returns true if any of the geometry is empty. +struct any_input_is_empty { + template + bool __device__ operator()(LhsType lhs, RhsType rhs) + { + return lhs.is_empty() || rhs.is_empty(); + } +}; + /** * @brief Computes distances between the multilinestring and multipolygons * - * @param multilinestrings Range to the multilinestring - * @param multipolygons Range to the multipolygon + * This is a load balanced distance compute kernel. Each thread compute exactly 1 pair of segments + * between the multilinestring and multipolygon. + * + * @tparam T type of the underlying coordinates + * @tparam index_t type of underlying offsets + * + * @param multilinestring_segments Range to the segments of the multilinestring + * @param multipolygon_segments Range to the segments of the multipolygon * @param thread_bounds Range to the boundary of thread partitions - * @param multilinestrings_segment_offsets Range to the indices where the first segment of each - * multilinestring begins - * @param multipolygons_segment_offsets Range to the indices where the first segment of each - * multipolygon begins * @param intersects A uint8_t array that indicates if the corresponding pair of multipoint and * multipolygon intersects * @param distances Output range of distances, pre-filled with std::numerical_limits::max() */ -template void __global__ -pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestrings, - MultiPolygonRange multipolygons, +pairwise_linestring_polygon_distance_kernel(MultiLinestringSegmentRange multilinestring_segments, + MultiPolygonSegmentRange multipolygon_segments, IndexRange thread_bounds, - IndexRange multilinestrings_segment_offsets, - IndexRange multipolygons_segment_offsets, uint8_t* intersects, OutputIt* distances) { - using T = typename MultiLinestringRange::element_t; - using index_t = iterator_value_type; - auto num_threads = thread_bounds[thread_bounds.size() - 1]; for (auto idx = blockDim.x * blockIdx.x + threadIdx.x; idx < num_threads; idx += blockDim.x * gridDim.x) { @@ -90,16 +96,17 @@ pairwise_linestring_polygon_distance_kernel(MultiLinestringRange multilinestring // Retrieve the number of segments in multilinestrings[geometry_id] auto num_segment_this_multilinestring = - multilinestrings.multilinestring_segment_count_begin()[geometry_id]; + multilinestring_segments.multigeometry_count_begin()[geometry_id]; // The segment id from the multilinestring this thread is computing (local_id + global_offset) auto multilinestring_segment_id = - local_idx % num_segment_this_multilinestring + multilinestrings_segment_offsets[geometry_id]; + local_idx % num_segment_this_multilinestring + + multilinestring_segments.multigeometry_offset_begin()[geometry_id]; // The segment id from the multipolygon this thread is computing (local_id + global_offset) - auto multipolygon_segment_id = - local_idx / num_segment_this_multilinestring + multipolygons_segment_offsets[geometry_id]; + auto multipolygon_segment_id = local_idx / num_segment_this_multilinestring + + multipolygon_segments.multigeometry_offset_begin()[geometry_id]; - auto [a, b] = multilinestrings.segment_begin()[multilinestring_segment_id]; - auto [c, d] = multipolygons.segment_begin()[multipolygon_segment_id]; + auto [a, b] = multilinestring_segments.begin()[multilinestring_segment_id]; + auto [c, d] = multipolygon_segments.begin()[multipolygon_segment_id]; atomicMin(&distances[geometry_id], sqrt(squared_segment_distance(a, b, c, d))); } @@ -119,18 +126,31 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri CUSPATIAL_EXPECTS(multilinestrings.size() == multipolygons.size(), "Must have the same number of input rows."); - if (multilinestrings.size() == 0) return distances_first; + auto size = multilinestrings.size(); + + if (size == 0) return distances_first; // Create a multipoint range from multilinestrings, computes intersection auto multipoints = multilinestrings.as_multipoint_range(); auto multipoint_intersects = point_polygon_intersects(multipoints, multipolygons, stream); + // Make views to the segments in the multilinestring + auto multilinestring_segments = multilinestrings._segments(stream); + auto multilinestring_segments_range = multilinestring_segments.view(); + auto multilinestring_segment_count_begin = + multilinestring_segments_range.multigeometry_count_begin(); + + // Make views to the segments in the multilinestring + auto multipolygon_segments = multipolygons._segments(stream); + auto multipolygon_segments_range = multipolygon_segments.view(); + auto multipolygon_segment_count_begin = multipolygon_segments_range.multigeometry_count_begin(); + // Compute the "boundary" of threads. Threads are partitioned based on the number of linestrings - // times the number of polygons in a multipoint-multipolygon pair. - auto segment_count_product_it = thrust::make_transform_iterator( - thrust::make_zip_iterator(multilinestrings.multilinestring_segment_count_begin(), - multipolygons.multipolygon_segment_count_begin()), - thrust::make_zip_function(thrust::multiplies{})); + // times the number of polygons in a multilinestring-multipolygon pair. + auto segment_count_product_it = + thrust::make_transform_iterator(thrust::make_zip_iterator(multilinestring_segment_count_begin, + multipolygon_segment_count_begin), + thrust::make_zip_function(thrust::multiplies{})); // Computes the "thread boundary" of each pair. This array partitions the thread range by // geometries. E.g. threadIdx within [thread_bounds[i], thread_bounds[i+1]) computes distances of @@ -143,48 +163,33 @@ OutputIt pairwise_linestring_polygon_distance(MultiLinestringRange multilinestri segment_count_product_it + thread_bounds.size() - 1, thrust::next(thread_bounds.begin())); - // Compute offsets to the first segment of each multilinestring and multipolygon - auto multilinestring_segment_offsets = - rmm::device_uvector(multilinestrings.num_multilinestrings() + 1, stream); - detail::zero_data_async( - multilinestring_segment_offsets.begin(), multilinestring_segment_offsets.end(), stream); - - auto multipolygon_segment_offsets = - rmm::device_uvector(multipolygons.num_multipolygons() + 1, stream); - detail::zero_data_async( - multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end(), stream); - - thrust::inclusive_scan(rmm::exec_policy(stream), - multilinestrings.multilinestring_segment_count_begin(), - multilinestrings.multilinestring_segment_count_begin() + - multilinestrings.num_multilinestrings(), - thrust::next(multilinestring_segment_offsets.begin())); - - thrust::inclusive_scan( - rmm::exec_policy(stream), - multipolygons.multipolygon_segment_count_begin(), - multipolygons.multipolygon_segment_count_begin() + multipolygons.num_multipolygons(), - thrust::next(multipolygon_segment_offsets.begin())); - // Initialize output range thrust::fill(rmm::exec_policy(stream), distances_first, - distances_first + multilinestrings.num_multilinestrings(), + distances_first + size, std::numeric_limits::max()); + // If any input multigeometries is empty, result is nan. + auto nan_it = thrust::make_constant_iterator(std::numeric_limits::quiet_NaN()); + thrust::transform_if(rmm::exec_policy(stream), + nan_it, + nan_it + size, + thrust::make_zip_iterator(multilinestrings.begin(), multipolygons.begin()), + distances_first, + thrust::identity{}, + thrust::make_zip_function(detail::any_input_is_empty{})); + auto num_threads = thread_bounds.back_element(stream); auto [tpb, num_blocks] = grid_1d(num_threads); - detail::pairwise_linestring_polygon_distance_kernel<<>>( - multilinestrings, - multipolygons, - range{thread_bounds.begin(), thread_bounds.end()}, - range{multilinestring_segment_offsets.begin(), multilinestring_segment_offsets.end()}, - range{multipolygon_segment_offsets.begin(), multipolygon_segment_offsets.end()}, - multipoint_intersects.begin(), - distances_first); + detail::pairwise_linestring_polygon_distance_kernel + <<>>(multilinestring_segments_range, + multipolygon_segments_range, + range{thread_bounds.begin(), thread_bounds.end()}, + multipoint_intersects.begin(), + distances_first); - return distances_first + multilinestrings.num_multilinestrings(); + return distances_first + size; } } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/functors.cuh b/cpp/include/cuspatial/detail/functors.cuh deleted file mode 100644 index 1f390753d..000000000 --- a/cpp/include/cuspatial/detail/functors.cuh +++ /dev/null @@ -1,142 +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 -#include - -#include - -namespace cuspatial { -namespace detail { - -/** - * @brief Given iterator a pair of offsets, return the number of elements between the offsets. - * - * Used to create iterator to geometry counts, such as `multi*_point_count_begin`, - * `multi*_segment_count_begin`. - * - * Example: - * pair of offsets: (0, 3), (3, 5), (5, 8) - * number of elements between offsets: 3, 2, 3 - * - * @tparam OffsetPairIterator Must be iterator type to thrust::pair of indices. - * @param p Iterator of thrust::pair of indices. - */ -struct offset_pair_to_count_functor { - template - CUSPATIAL_HOST_DEVICE auto operator()(OffsetPairIterator p) - { - return thrust::get<1>(p) - thrust::get<0>(p); - } -}; - -/** - * @brief Convert counts of points to counts of segments in a linestring. - * - * A Multilinestring is composed of a series of Linestrings. Each Linestring is composed of a - * segments. The number of segments in a multilinestring is the number of points in the - * multilinestring minus the number of linestrings. - * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. - * - * Used to create segment count iterators, such as `multi*_segment_count_begin`. - * - * @tparam IndexPair Must be iterator to a pair of counts - * @param n_point_linestring_pair A pair of counts, the first is the number of points, the second is - * the number of linestrings. - */ -struct point_count_to_segment_count_functor { - template - CUSPATIAL_HOST_DEVICE auto operator()(IndexPair n_point_linestring_pair) - { - auto nPoints = thrust::get<0>(n_point_linestring_pair); - auto nLinestrings = thrust::get<1>(n_point_linestring_pair); - return nPoints - nLinestrings; - } -}; - -/** - * @brief Given an offset iterator it, returns an iterator of the distance between it and an input - * index i - * - * @tparam OffsetIterator Iterator type to the offset - * - * Caveats: This has a strong assumption that the Multilinestring does not contain empty - * linestrings. While each non-empty linestring in the multilinestring represents 1 extra segment, - * an empty multilinestring does not introduce any extra segments since it does not contain any - * points. - * - * Used to create iterator to segment offsets, such as `segment_offset_begin`. - */ -template -struct to_distance_iterator { - OffsetIterator begin; - - template - CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) - { - return begin[i] - i; - } -}; - -/// Deduction guide for to_distance_iterator -template -to_distance_iterator(OffsetIterator) -> to_distance_iterator; - -/** - * @brief Return a segment from the a partitioned range of points - * - * Used in a counting transform iterator. Given an index of the segment, offset it by the number of - * skipped segments preceding i in the partitioned range of points. Dereference the corresponding - * point and the point following to make a segment. - * - * Used to create iterator to segments, such as `segment_begin`. - * - * @tparam OffsetIterator the iterator type indicating partitions of the point range. - * @tparam CoordinateIterator the iterator type to the point range. - */ -template -struct to_valid_segment_functor { - using element_t = iterator_vec_base_type; - - OffsetIterator begin; - OffsetIterator end; - CoordinateIterator point_begin; - - template - CUSPATIAL_HOST_DEVICE segment operator()(IndexType i) - { - auto kit = thrust::upper_bound(thrust::seq, begin, end, i); - auto k = thrust::distance(begin, kit); - auto pid = i + k - 1; - - return segment{point_begin[pid], point_begin[pid + 1]}; - } -}; - -/// Deduction guide for to_valid_segment_functor -template -to_valid_segment_functor(OffsetIterator, OffsetIterator, CoordinateIterator) - -> to_valid_segment_functor; - -} // namespace detail -} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/iterator_factory.cuh b/cpp/include/cuspatial/detail/iterator_factory.cuh new file mode 100644 index 000000000..8f66de0ed --- /dev/null +++ b/cpp/include/cuspatial/detail/iterator_factory.cuh @@ -0,0 +1,203 @@ +/* + * 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 + +#include + +namespace cuspatial { +namespace detail { + +/** + * @internal + * @brief Helper to create a `transform_iterator` that transforms sequential values. + */ +template +inline CUSPATIAL_HOST_DEVICE auto make_counting_transform_iterator(IndexType start, UnaryFunction f) +{ + return thrust::make_transform_iterator(thrust::make_counting_iterator(start), f); +} + +/** + * @internal + * @brief Helper to convert a tuple of elements into a `vec_2d` + */ +template > +struct tuple_to_vec_2d { + __device__ VectorType operator()(thrust::tuple const& pos) + { + return VectorType{thrust::get<0>(pos), thrust::get<1>(pos)}; + } +}; + +/** + * @internal + * @brief Helper to convert a `vec_2d` into a tuple of elements + */ +template > +struct vec_2d_to_tuple { + __device__ thrust::tuple operator()(VectorType const& xy) + { + return thrust::make_tuple(xy.x, xy.y); + } +}; + +/** + * @internal + * @brief Helper to convert a `box` into a tuple of elements + */ +template > +struct box_to_tuple { + __device__ thrust::tuple operator()(Box const& box) + { + return thrust::make_tuple(box.v1.x, box.v1.y, box.v2.x, box.v2.y); + } +}; + +/** + * @internal + * @brief Helper to convert a tuple of vec_2d into a `box` + */ +template > +struct vec_2d_tuple_to_box { + __device__ box operator()(thrust::tuple pair) + { + auto v1 = thrust::get<0>(pair); + auto v2 = thrust::get<1>(pair); + return {v1, v2}; + } +}; + +/** + * @internal + * @brief Generic to convert any iterator pointing to interleaved xy range into + * iterator of vec_2d. + * + * This generic version does not use of vectorized load. + * + * @pre `Iter` has operator[] defined. + * @pre std::iterator_traits::value_type is convertible to `T`. + */ +template +struct interleaved_to_vec_2d { + using element_t = typename std::iterator_traits::value_type; + using value_type = vec_2d; + Iter it; + constexpr interleaved_to_vec_2d(Iter it) : it{it} {} + + CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) + { + return vec_2d{it[2 * i], it[2 * i + 1]}; + } +}; + +/** + * @brief Specialization for thrust iterators conforming to `contiguous_iterator`. (including raw + * pointer) + * + * This iterator specific version uses vectorized load. + * + * @throw cuspatial::logic_error if `Iter` is not aligned to type `vec_2d` + * @pre `Iter` is a `contiguous_iterator` (including raw pointer). + */ +template +struct interleaved_to_vec_2d>> { + using element_t = typename std::iterator_traits::value_type; + using value_type = vec_2d; + + element_t const* ptr; + + constexpr interleaved_to_vec_2d(Iter it) : ptr{&thrust::raw_reference_cast(*it)} + { + CUSPATIAL_EXPECTS(!((intptr_t)ptr % alignof(vec_2d)), + "Misaligned interleaved data."); + } + + CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) + { + auto const aligned = + static_cast(__builtin_assume_aligned(ptr + 2 * i, 2 * sizeof(element_t))); + return vec_2d{aligned[0], aligned[1]}; + } +}; + +/** + * @internal + * @brief Functor to transform an index to strided index. + */ +template +struct strided_functor { + auto __device__ operator()(std::size_t i) { return i * stride; } +}; + +/** + * @internal + * @brief Functor to transform an index into a geometry ID determined by a range of offsets. + */ +template +struct index_to_geometry_id { + GeometryIter geometry_begin; + GeometryIter geometry_end; + + index_to_geometry_id(GeometryIter begin, GeometryIter end) + : geometry_begin(begin), geometry_end(end) + { + } + + CUSPATIAL_HOST_DEVICE auto operator()(IndexT idx) + { + return thrust::distance(geometry_begin, + thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx)); + } +}; + +/** + * @brief Given iterator a pair of offsets, return the number of elements between the offsets. + * + * Used to create iterator to geometry counts, such as `multi*_point_count_begin`, + * `multi*_segment_count_begin`. + * + * Example: + * pair of offsets: (0, 3), (3, 5), (5, 8) + * number of elements between offsets: 3, 2, 3 + * + * @tparam OffsetPairIterator Must be iterator type to thrust::pair of indices. + * @param p Iterator of thrust::pair of indices. + */ +struct offset_pair_to_count_functor { + template + CUSPATIAL_HOST_DEVICE auto operator()(OffsetPairIterator p) + { + return thrust::get<1>(p) - thrust::get<0>(p); + } +}; + +} // namespace detail +} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/multilinestring_segment.cuh b/cpp/include/cuspatial/detail/multilinestring_segment.cuh new file mode 100644 index 000000000..d57f125dc --- /dev/null +++ b/cpp/include/cuspatial/detail/multilinestring_segment.cuh @@ -0,0 +1,127 @@ +/* + * 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 + +namespace cuspatial { +namespace detail { + +/** + * @internal + * @brief Functor that returns true if current value is greater than 0. + */ +template +struct greater_than_zero_functor { + __device__ IndexType operator()(IndexType x) const { return x > 0; } +}; + +/** + * @internal + * @brief Owning class to provide iterators to segments in a multilinestring range + * + * The owned memory in this struct is `_non_empty_linestring_prefix_sum`, which equals the + * number of linestrings in the multilinestring range plus 1. This vector holds the number + * of non empty linestrings that precedes the current linestring. + * + * This class is only meant for tracking the life time of the owned memory. To access the + * segment iterators, call `view()` function to create a non-owning object of this class. + * + * For detailed explanation on the implementation of the segment iterators, see documentation + * of `multilinestring_segment_range`. + * + * @note To use this class with a multipolygon range, cast the multipolygon range as a + * multilinestring range. + * + * TODO: Optimization: for range that does not contain any empty linestrings, + * `_non_empty_linestring_prefix_sum` can be substituted with `counting_iterator`. + * + * @tparam MultilinestringRange The multilinestring range to initialize this class with. + */ +template +class multilinestring_segment { + using index_t = iterator_value_type; + + public: + /** + * @brief Construct a new multilinestring segment object + * + * @note multilinestring_segment is always internal use, thus memory consumed is always + * temporary, therefore always use default device memory resource. + * + * @param parent The parent multilinestring object to construct from + * @param stream The stream to perform computation on + */ + multilinestring_segment(MultilinestringRange parent, rmm::cuda_stream_view stream) + : _parent(parent), _non_empty_linestring_prefix_sum(parent.num_linestrings() + 1, stream) + { + auto offset_range = ::cuspatial::range{_parent.part_offset_begin(), _parent.part_offset_end()}; + auto count_begin = thrust::make_transform_iterator( + thrust::make_zip_iterator(offset_range.begin(), thrust::next(offset_range.begin())), + offset_pair_to_count_functor{}); + + auto count_greater_than_zero = + thrust::make_transform_iterator(count_begin, greater_than_zero_functor{}); + + zero_data_async( + _non_empty_linestring_prefix_sum.begin(), _non_empty_linestring_prefix_sum.end(), stream); + + thrust::inclusive_scan(rmm::exec_policy(stream), + count_greater_than_zero, + count_greater_than_zero + _parent.num_linestrings(), + thrust::next(_non_empty_linestring_prefix_sum.begin())); + + _num_segments = _parent.num_points() - _non_empty_linestring_prefix_sum.element( + _non_empty_linestring_prefix_sum.size() - 1, stream); + } + + /** + * @brief Return a non-owning `multilinestring_segment_range` object from this class + * + * @return multilinestring_segment_range + */ + auto view() + { + auto index_range = ::cuspatial::range{_non_empty_linestring_prefix_sum.begin(), + _non_empty_linestring_prefix_sum.end()}; + return multilinestring_segment_range{ + _parent, index_range, _num_segments}; + } + + private: + MultilinestringRange _parent; + index_t _num_segments; + rmm::device_uvector _non_empty_linestring_prefix_sum; +}; + +} // namespace detail + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh index 5983efc00..3b0aa597f 100644 --- a/cpp/include/cuspatial/detail/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/detail/range/multilinestring_range.cuh @@ -16,14 +16,8 @@ #pragma once -#include -#include -#include -#include -#include - #include -#include +#include #include #include #include @@ -31,7 +25,13 @@ #include #include +#include +#include +#include #include +#include +#include +#include #include #include @@ -115,13 +115,6 @@ multilinestring_range::num_points() return thrust::distance(_point_begin, _point_end); } -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::num_segments() -{ - return num_points() - num_linestrings(); -} - template CUSPATIAL_HOST_DEVICE auto multilinestring_range::multilinestring_begin() @@ -219,9 +212,7 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range @@ -231,29 +222,11 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range -CUSPATIAL_HOST_DEVICE auto multilinestring_range:: - multilinestring_segment_count_begin() -{ - auto n_point_linestring_pair_it = thrust::make_zip_iterator( - multilinestring_point_count_begin(), multilinestring_linestring_count_begin()); - return thrust::make_transform_iterator(n_point_linestring_pair_it, - detail::point_count_to_segment_count_functor{}); -} - -template -CUSPATIAL_HOST_DEVICE auto multilinestring_range:: - multilinestring_segment_count_end() -{ - return multilinestring_segment_count_begin() + num_multilinestrings(); -} - template CUSPATIAL_HOST_DEVICE auto multilinestring_range:: multilinestring_linestring_count_begin() { - auto paired_it = thrust::make_zip_iterator(_geometry_begin, thrust::next(_geometry_begin)); - return thrust::make_transform_iterator(paired_it, detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(_geometry_begin); } template @@ -264,34 +237,10 @@ CUSPATIAL_HOST_DEVICE auto multilinestring_range -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_begin() -{ - return detail::make_counting_transform_iterator( - 0, - detail::to_valid_segment_functor{ - this->segment_offset_begin(), this->segment_offset_end(), _point_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_end() -{ - return segment_begin() + num_segments(); -} - -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_offset_begin() -{ - return detail::make_counting_transform_iterator(0, detail::to_distance_iterator{_part_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multilinestring_range::segment_offset_end() +auto multilinestring_range::_segments( + rmm::cuda_stream_view stream) { - return segment_offset_begin() + thrust::distance(_part_begin, _part_end); + return multilinestring_segment{*this, stream}; } template diff --git a/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh new file mode 100644 index 000000000..58817c46e --- /dev/null +++ b/cpp/include/cuspatial/detail/range/multilinestring_segment_range.cuh @@ -0,0 +1,267 @@ +/* + * 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 + +namespace cuspatial { +namespace detail { + +/** + * @brief Computes the offsets to the starting segment per linestring + * + * The point indices and segment indices are correlated, but in a different index space. + * For example: + * ``` + * {0, 3} + * {0, 3, 3, 6} + * {A, B, C, X, Y, Z} + * ``` + * + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * points: A B C X Y Z + * pid: 0 1 2 3 4 5 + * ``` + * + * The original {0, 3, 3, 6} offsets are in the point index space. For example: + * The first and past the last point of the first linestring is at point index 0 and 3 (A, X). + * The first and past the last point of the second linestring is at point index 3 and 3 (empty), + * and so on. + * + * The transformed segment offsets {0, 2, 2, 4} are in the segment index space. For example: + * The first and past the last segment of the first linestring is at segment index 0 and 2 ((AB), + * (XY)). + * The first and past the last segment of the second linestring is at segment index 2 and 2 + * (empty), and so on. + * + * @tparam OffsetIterator Iterator type to the offset + */ +template +struct point_offset_to_segment_offset { + OffsetIterator part_offset_begin; + CountIterator non_empty_linestrings_count_begin; + + template + CUSPATIAL_HOST_DEVICE auto operator()(IndexType i) + { + return part_offset_begin[i] - non_empty_linestrings_count_begin[i]; + } +}; + +/// Deduction guide +template +point_offset_to_segment_offset(OffsetIterator, CountIterator) + -> point_offset_to_segment_offset; + +/** + * @brief Given a segment index, return the corresponding segment + * + * Given a segment index, first find its corresponding part index by performing a binary search in + * the segment offsets range. Then, skip the segment index by the number of non empty linestrings + * that precedes the current linestring to find point index to the first point of the segment. + * Dereference this point and the following point to construct the segment. + * + * @tparam OffsetIterator Iterator to the segment offsets + * @tparam CountIterator Iterator the the range of the prefix sum of non empty linestrings + * @tparam CoordinateIterator Iterator to the point range + */ +template +struct to_valid_segment_functor { + using element_t = iterator_vec_base_type; + + OffsetIterator segment_offset_begin; + OffsetIterator segment_offset_end; + CountIterator non_empty_partitions_begin; + CoordinateIterator point_begin; + + template + CUSPATIAL_HOST_DEVICE segment operator()(IndexType sid) + { + auto kit = + thrust::prev(thrust::upper_bound(thrust::seq, segment_offset_begin, segment_offset_end, sid)); + auto part_id = thrust::distance(segment_offset_begin, kit); + auto preceding_non_empty_linestrings = non_empty_partitions_begin[part_id]; + auto pid = sid + preceding_non_empty_linestrings; + + return segment{point_begin[pid], point_begin[pid + 1]}; + } +}; + +/// Deduction guide +template +to_valid_segment_functor(OffsetIterator, OffsetIterator, CountIterator, CoordinateIterator) + -> to_valid_segment_functor; + +/** + * @brief Non-owning object to `multilinestring_segment` class + * + * A `multilinestring_segment_range` provide views into the segments of a multilinestring. + * The segments of a multilinestring has a near 1-1 mapping to the points of the multilinestring, + * except that the last point of a linestring and the first point of the next linestring does not + * form a valid segment. For example, the below multilinestring (points are denoted a letters): + * + * ``` + * {0, 2} + * {0, 3, 6} + * {A, B, C, X, Y, Z} + * ``` + * + * contains 6 points, but only 4 segments. AB, BC, XY and YZ. + * If we assign an index to all four segments, and an index to all points: + * + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * points: A B C X Y Z + * pid: 0 1 2 3 4 5 + * ``` + * + * Notice that if we "skip" the segment index by a few steps, it can correctly find the + * corresponding point index of the starting point of the segment. For example: skipping sid==0 (AB) + * by 0 steps, finds the starting point of A (pid==0) skipping sid==2 (XY) by 1 step, finds the + * starting point of X (pid==3) + * + * Intuitively, the *steps to skip* equals the number of linestrings that precedes the linestring + * that the current segment is in. This is because every linestring adds an "invalid" segment to the + * preceding linestring. However, consider the following edge case that contains empty linestrings: + * + * ``` + * {0, 3} + * {0, 3, 3, 6} + * {A, B, C, X, Y, Z} + * ``` + * + * For segment XY, there are 2 linestrings that precedes its linestring ((0, 3) and (3, 3)). + * However, we cannot skip the sid of XY by 2 to get its starting point index. This is because the + * empty linestring in between does not introduce the "invalid" segment. Therefore, the correct + * steps to skip equals the number of *non-empty* linestrings that precedes the current linestring + * that the segment is in. + * + * Concretely, with the above example: + * ``` + * segments: AB BC XY YZ + * sid: 0 1 2 3 + * num_preceding_non_empty_linestrings: 0 0 1 1 + * skipped sid (pid): 0 0 3 4 + * starting point: A B X Y + * ``` + * + * @tparam ParentRange The multilinestring range to construct from + * @tparam IndexRange The range to the prefix sum of the non empty linestring counts + */ +template +class multilinestring_segment_range { + using index_t = typename IndexRange::value_type; + + public: + multilinestring_segment_range(ParentRange parent, + IndexRange non_empty_geometry_prefix_sum, + index_t num_segments) + : _parent(parent), + _non_empty_geometry_prefix_sum(non_empty_geometry_prefix_sum), + _num_segments(num_segments) + { + } + + /// Returns the number of segments in the multilinestring + CUSPATIAL_HOST_DEVICE index_t num_segments() { return _num_segments; } + + /// Returns starting iterator to the range of the starting segment index per + /// multilinestring or multipolygon + CUSPATIAL_HOST_DEVICE auto multigeometry_offset_begin() + { + return thrust::make_permutation_iterator(_per_linestring_offset_begin(), + _parent.geometry_offsets_begin()); + } + + /// Returns end iterator to the range of the starting segment index per multilinestring + /// or multipolygon + CUSPATIAL_HOST_DEVICE auto multigeometry_offset_end() + { + return multigeometry_offset_begin() + _parent.num_multilinestrings() + 1; + } + + /// Returns starting iterator to the range of the number of segments per multilinestring of + /// multipolygon + CUSPATIAL_HOST_DEVICE auto multigeometry_count_begin() + { + return make_element_count_iterator(multigeometry_offset_begin()); + } + + /// Returns end iterator to the range of the number of segments per multilinestring of + /// multipolygon + CUSPATIAL_HOST_DEVICE auto multigeometry_count_end() + { + return multigeometry_count_begin() + _parent.num_multilinestrings(); + } + + /// Returns the iterator to the first segment of the geometry range + /// See `to_valid_segment_functor` for implementation detail + CUSPATIAL_HOST_DEVICE auto begin() + { + return make_counting_transform_iterator( + 0, + to_valid_segment_functor{_per_linestring_offset_begin(), + _per_linestring_offset_end(), + _non_empty_geometry_prefix_sum.begin(), + _parent.point_begin()}); + } + + /// Returns the iterator to the past the last segment of the geometry range + CUSPATIAL_HOST_DEVICE auto end() { return begin() + _num_segments; } + + private: + ParentRange _parent; + IndexRange _non_empty_geometry_prefix_sum; + index_t _num_segments; + + /// Returns begin iterator to the index that points to the starting index for each linestring + /// See documentation of `to_segment_offset_iterator` for detail. + CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_begin() + { + return make_counting_transform_iterator( + 0, + point_offset_to_segment_offset{_parent.part_offset_begin(), + _non_empty_geometry_prefix_sum.begin()}); + } + + /// Returns end iterator to the index that points to the starting index for each linestring + CUSPATIAL_HOST_DEVICE auto _per_linestring_offset_end() + { + return _per_linestring_offset_begin() + _non_empty_geometry_prefix_sum.size(); + } +}; + +template +multilinestring_segment_range(ParentRange, IndexRange, typename IndexRange::value_type, bool) + -> multilinestring_segment_range; + +} // namespace detail + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh index 85781e55b..567105fd1 100644 --- a/cpp/include/cuspatial/detail/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipolygon_range.cuh @@ -17,7 +17,7 @@ #pragma once #include -#include +#include #include #include #include @@ -27,8 +27,11 @@ #include #include +#include + #include #include +#include #include #include #include @@ -154,16 +157,6 @@ multipolygon_range::n return thrust::distance(_point_begin, _point_end); } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::num_segments() -{ - return num_points() - num_rings(); -} - template :: auto multipolygon_point_offset_it = thrust::make_permutation_iterator( _ring_begin, thrust::make_permutation_iterator(_part_begin, _geometry_begin)); - auto point_offset_pair_it = thrust::make_zip_iterator(multipolygon_point_offset_it, - thrust::next(multipolygon_point_offset_it)); - - return thrust::make_transform_iterator(point_offset_pair_it, - detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(multipolygon_point_offset_it); } template :: auto multipolygon_ring_offset_it = thrust::make_permutation_iterator(_part_begin, _geometry_begin); - auto ring_offset_pair_it = thrust::make_zip_iterator(multipolygon_ring_offset_it, - thrust::next(multipolygon_ring_offset_it)); - - return thrust::make_transform_iterator(ring_offset_pair_it, - detail::offset_pair_to_count_functor{}); + return make_element_count_iterator(multipolygon_ring_offset_it); } template :: return multipolygon_ring_count_begin() + num_multipolygons(); } -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - multipolygon_segment_count_begin() -{ - auto multipolygon_point_ring_count_it = - thrust::make_zip_iterator(multipolygon_point_count_begin(), multipolygon_ring_count_begin()); - - return thrust::make_transform_iterator(multipolygon_point_ring_count_it, - detail::point_count_to_segment_count_functor{}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - multipolygon_segment_count_end() -{ - return multipolygon_segment_count_begin() + num_multipolygons(); -} - template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::segment_begin() -{ - return detail::make_counting_transform_iterator( - 0, - detail::to_valid_segment_functor{ - this->subtracted_ring_begin(), this->subtracted_ring_end(), _point_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::segment_end() -{ - return segment_begin() + num_segments(); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range:: - subtracted_ring_begin() -{ - return detail::make_counting_transform_iterator(0, detail::to_distance_iterator{_ring_begin}); -} - -template -CUSPATIAL_HOST_DEVICE auto -multipolygon_range::subtracted_ring_end() +auto multipolygon_range::_segments( + rmm::cuda_stream_view stream) { - return subtracted_ring_begin() + thrust::distance(_ring_begin, _ring_end); + auto multilinestring_range = this->as_multilinestring_range(); + return multilinestring_segment{multilinestring_range, stream}; } template 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1), \ - "Each linestring must have at least two vertices"); + "Polygon offsets must contain at least one (1) value"); /** * @brief Macro for validating the data array sizes for multilinestrings. @@ -46,7 +46,6 @@ * Raises an exception if any of the following are false: * - The number of multilinestring offsets is greater than zero. * - The number of linestring offsets is greater than zero. - * - There are at least two vertices per linestring offset. * * Multilinestrings follow [GeoArrow data layout][1]. Offsets arrays have one more element than the * number of items in the array. The last offset is always the sum of the previous offset and the @@ -69,8 +68,6 @@ * Raises an exception if any of the following are false: * - The number of polygon offsets is greater than zero. * - The number of ring offsets is greater than zero. - * - There is at least one ring offset per polygon offset. - * - There are at least four vertices per ring offset. * * Polygons follow [GeoArrow data layout][1]. Offsets arrays (polygons and rings) have one more * element than the number of items in the array. The last offset is always the sum of the previous @@ -78,22 +75,21 @@ * last ring offset plus the number of rings in the last polygon. See * [Arrow Variable-Size Binary layout](2). Note that an empty list still has one offset: {0}. * - * Rings are assumed to be closed (closed means the first and last vertices of - * each ring are equal). Therefore rings must have at least 4 vertices. + * The following are not explicitly checked in this macro, but is always assumed in cuspatial: + * + * 1. Polygon can contain zero or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are + * equal). Rings can also be empty. Therefore each ring must contain zero or four or more vertices. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout */ -#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( \ - num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ - "Polygon ring offsets must contain at least one (1) value"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets >= num_poly_offsets, \ - "Each polygon must have at least one (1) ring"); \ - CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_points >= 4 * (num_poly_ring_offsets - 1), \ - "Each ring must have at least four (4) vertices"); +#define CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( \ + num_poly_points, num_poly_offsets, num_poly_ring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_offsets > 0, \ + "Polygon offsets must contain at least one (1) value"); \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_poly_ring_offsets > 0, \ + "Polygon ring offsets must contain at least one (1) value"); /** * @brief Macro for validating the data array sizes for a multipolygon. @@ -102,8 +98,6 @@ * - The number of multipolygon offsets is greater than zero. * - The number of polygon offsets is greater than zero. * - The number of ring offsets is greater than zero. - * - There is at least one ring offset per polygon offset. - * - There are at least four vertices per ring offset. * * MultiPolygons follow [GeoArrow data layout][1]. Offsets arrays (polygons and rings) have one more * element than the number of items in the array. The last offset is always the sum of the previous @@ -111,8 +105,12 @@ * last ring offset plus the number of rings in the last polygon. See * [Arrow Variable-Size Binary layout](2). Note that an empty list still has one offset: {0}. * - * Rings are assumed to be closed (closed means the first and last vertices of - * each ring are equal). Therefore rings must have at least 4 vertices. + * The following are not explicitly checked in this macro, but is always assumed in cuspatial: + * + * 1. Polygon can contain zero or more rings. A polygon with zero rings is an empty polygon. + * 2. Rings are assumed to be closed (closed means the first and last vertices of each ring are + * equal). Rings can also be empty. Therefore each ring must contain zero or four or more + * coordinates. * * [1]: https://github.com/geoarrow/geoarrow/blob/main/format.md * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout diff --git a/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh b/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh index 872135f72..36304a7bc 100644 --- a/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh +++ b/cpp/include/cuspatial/geometry_collection/multilinestring_ref.cuh @@ -38,6 +38,9 @@ class multilinestring_ref { /// Return the number of linestrings in the multilinestring. CUSPATIAL_HOST_DEVICE auto size() const { return num_linestrings(); } + /// Return true if this multilinestring contains 0 linestrings. + CUSPATIAL_HOST_DEVICE bool is_empty() const { return num_linestrings() == 0; } + /// Return iterator to the first linestring. CUSPATIAL_HOST_DEVICE auto part_begin() const; /// Return iterator to one past the last linestring. diff --git a/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh b/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh index e4ff5010f..d9972d9b1 100644 --- a/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh +++ b/cpp/include/cuspatial/geometry_collection/multipolygon_ref.cuh @@ -41,6 +41,9 @@ class multipolygon_ref { /// Return the number of polygons in the multipolygon. CUSPATIAL_HOST_DEVICE auto size() const { return num_polygons(); } + /// Returns true if the multipolygon contains 0 geometries. + CUSPATIAL_HOST_DEVICE bool is_empty() const { return num_polygons() == 0; } + /// Return iterator to the first polygon. CUSPATIAL_HOST_DEVICE auto part_begin() const; /// Return iterator to one past the last polygon. diff --git a/cpp/include/cuspatial/iterator_factory.cuh b/cpp/include/cuspatial/iterator_factory.cuh index 1f026512c..49313298b 100644 --- a/cpp/include/cuspatial/iterator_factory.cuh +++ b/cpp/include/cuspatial/iterator_factory.cuh @@ -16,171 +16,14 @@ #pragma once -#include -#include -#include -#include +#include -#include -#include -#include #include #include #include #include -#include -#include - -#include namespace cuspatial { -namespace detail { - -/** - * @internal - * @brief Helper to create a `transform_iterator` that transforms sequential values. - */ -template -inline CUSPATIAL_HOST_DEVICE auto make_counting_transform_iterator(IndexType start, UnaryFunction f) -{ - return thrust::make_transform_iterator(thrust::make_counting_iterator(start), f); -} - -/** - * @internal - * @brief Helper to convert a tuple of elements into a `vec_2d` - */ -template > -struct tuple_to_vec_2d { - __device__ VectorType operator()(thrust::tuple const& pos) - { - return VectorType{thrust::get<0>(pos), thrust::get<1>(pos)}; - } -}; - -/** - * @internal - * @brief Helper to convert a `vec_2d` into a tuple of elements - */ -template > -struct vec_2d_to_tuple { - __device__ thrust::tuple operator()(VectorType const& xy) - { - return thrust::make_tuple(xy.x, xy.y); - } -}; - -/** - * @internal - * @brief Helper to convert a `box` into a tuple of elements - */ -template > -struct box_to_tuple { - __device__ thrust::tuple operator()(Box const& box) - { - return thrust::make_tuple(box.v1.x, box.v1.y, box.v2.x, box.v2.y); - } -}; - -/** - * @internal - * @brief Helper to convert a tuple of vec_2d into a `box` - */ -template > -struct vec_2d_tuple_to_box { - __device__ box operator()(thrust::tuple pair) - { - auto v1 = thrust::get<0>(pair); - auto v2 = thrust::get<1>(pair); - return {v1, v2}; - } -}; - -/** - * @internal - * @brief Generic to convert any iterator pointing to interleaved xy range into - * iterator of vec_2d. - * - * This generic version does not use of vectorized load. - * - * @pre `Iter` has operator[] defined. - * @pre std::iterator_traits::value_type is convertible to `T`. - */ -template -struct interleaved_to_vec_2d { - using element_t = typename std::iterator_traits::value_type; - using value_type = vec_2d; - Iter it; - constexpr interleaved_to_vec_2d(Iter it) : it{it} {} - - CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) - { - return vec_2d{it[2 * i], it[2 * i + 1]}; - } -}; - -/** - * @brief Specialization for thrust iterators conforming to `contiguous_iterator`. (including raw - * pointer) - * - * This iterator specific version uses vectorized load. - * - * @throw cuspatial::logic_error if `Iter` is not aligned to type `vec_2d` - * @pre `Iter` is a `contiguous_iterator` (including raw pointer). - */ -template -struct interleaved_to_vec_2d>> { - using element_t = typename std::iterator_traits::value_type; - using value_type = vec_2d; - - element_t const* ptr; - - constexpr interleaved_to_vec_2d(Iter it) : ptr{&thrust::raw_reference_cast(*it)} - { - CUSPATIAL_EXPECTS(!((intptr_t)ptr % alignof(vec_2d)), - "Misaligned interleaved data."); - } - - CUSPATIAL_HOST_DEVICE value_type operator()(std::size_t i) - { - auto const aligned = - static_cast(__builtin_assume_aligned(ptr + 2 * i, 2 * sizeof(element_t))); - return vec_2d{aligned[0], aligned[1]}; - } -}; - -/** - * @internal - * @brief Functor to transform an index to strided index. - */ -template -struct strided_functor { - auto __device__ operator()(std::size_t i) { return i * stride; } -}; - -/** - * @internal - * @brief Functor to transform an index into a geometry ID determined by a range of offsets. - */ -template -struct index_to_geometry_id { - GeometryIter geometry_begin; - GeometryIter geometry_end; - - index_to_geometry_id(GeometryIter begin, GeometryIter end) - : geometry_begin(begin), geometry_end(end) - { - } - - CUSPATIAL_HOST_DEVICE auto operator()(IndexT idx) - { - return thrust::distance(geometry_begin, - thrust::upper_bound(thrust::seq, geometry_begin, geometry_end, idx)); - } -}; - -} // namespace detail /** * @addtogroup type_factories @@ -424,6 +267,13 @@ auto make_geometry_id_iterator(GeometryIter geometry_offsets_begin, std::distance(geometry_offsets_begin, geometry_offsets_end))); } +template +CUSPATIAL_HOST_DEVICE auto make_element_count_iterator(OffsetIterator iter) +{ + auto zipped_it = thrust::make_zip_iterator(iter, thrust::next(iter)); + return thrust::make_transform_iterator(zipped_it, detail::offset_pair_to_count_functor{}); +} + /** * @} // end of doxygen group */ diff --git a/cpp/include/cuspatial/range/multilinestring_range.cuh b/cpp/include/cuspatial/range/multilinestring_range.cuh index b79b2ae77..b4bb6eaf5 100644 --- a/cpp/include/cuspatial/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/range/multilinestring_range.cuh @@ -23,6 +23,8 @@ #include #include +#include + #include namespace cuspatial { @@ -81,9 +83,6 @@ class multilinestring_range { /// Return the total number of points in the array. CUSPATIAL_HOST_DEVICE auto num_points(); - /// Return the total number of segments in the array. - CUSPATIAL_HOST_DEVICE auto num_segments(); - /// Return the iterator to the first multilinestring in the range. CUSPATIAL_HOST_DEVICE auto multilinestring_begin(); @@ -154,23 +153,16 @@ class multilinestring_range { /// Returns an iterator to the counts of segments per multilinestring CUSPATIAL_HOST_DEVICE auto multilinestring_point_count_end(); - /// Returns an iterator to the counts of segments per multilinestring - CUSPATIAL_HOST_DEVICE auto multilinestring_segment_count_begin(); - - /// Returns an iterator to the counts of points per multilinestring - CUSPATIAL_HOST_DEVICE auto multilinestring_segment_count_end(); - /// Returns an iterator to the counts of points per multilinestring CUSPATIAL_HOST_DEVICE auto multilinestring_linestring_count_begin(); /// Returns an iterator to the counts of points per multilinestring CUSPATIAL_HOST_DEVICE auto multilinestring_linestring_count_end(); - /// Returns an iterator to the start of the segment - CUSPATIAL_HOST_DEVICE auto segment_begin(); - - /// Returns an iterator to the end of the segment - CUSPATIAL_HOST_DEVICE auto segment_end(); + /// @internal + /// Returns the owning class that provides views into the segments of the multilinestring range + /// Can only be constructed on host + auto _segments(rmm::cuda_stream_view); /// Returns the `multilinestring_idx`th multilinestring in the range. template @@ -198,9 +190,6 @@ class multilinestring_range { VecIterator _point_begin; VecIterator _point_end; - CUSPATIAL_HOST_DEVICE auto segment_offset_begin(); - CUSPATIAL_HOST_DEVICE auto segment_offset_end(); - private: /// @internal /// Return the iterator to the part index where the point locates. diff --git a/cpp/include/cuspatial/range/multipolygon_range.cuh b/cpp/include/cuspatial/range/multipolygon_range.cuh index 1b8947345..99b2843b2 100644 --- a/cpp/include/cuspatial/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/range/multipolygon_range.cuh @@ -16,14 +16,16 @@ #pragma once -#include - #include #include #include #include #include +#include + +#include + namespace cuspatial { /** @@ -96,9 +98,6 @@ class multipolygon_range { /// Return the total number of points in the array. CUSPATIAL_HOST_DEVICE auto num_points(); - /// Return the total number of segments in the array. - CUSPATIAL_HOST_DEVICE auto num_segments(); - /// Return the iterator to the first multipolygon in the range. CUSPATIAL_HOST_DEVICE auto multipolygon_begin(); @@ -117,6 +116,12 @@ class multipolygon_range { /// Return the iterator to the one past the last point in the range. CUSPATIAL_HOST_DEVICE auto point_end(); + /// Return the iterator to the first geometry offset in the range. + CUSPATIAL_HOST_DEVICE auto geometry_offsets_begin() { return _part_begin; } + + /// Return the iterator to the one past the last geometry offset in the range. + CUSPATIAL_HOST_DEVICE auto geometry_offsets_end() { return _part_end; } + /// Return the iterator to the first part offset in the range. CUSPATIAL_HOST_DEVICE auto part_offset_begin() { return _part_begin; } @@ -175,16 +180,10 @@ class multipolygon_range { /// Returns the one past the iterator to the number of rings of the last multipolygon CUSPATIAL_HOST_DEVICE auto multipolygon_ring_count_end(); - /// Returns an iterator to the number of segments of the first multipolygon - CUSPATIAL_HOST_DEVICE auto multipolygon_segment_count_begin(); - /// Returns the one past the iterator to the number of segments of the last multipolygon - CUSPATIAL_HOST_DEVICE auto multipolygon_segment_count_end(); - - /// Returns an iterator to the start of the segment - CUSPATIAL_HOST_DEVICE auto segment_begin(); - - /// Returns an iterator to the end of the segment - CUSPATIAL_HOST_DEVICE auto segment_end(); + /// @internal + /// Returns the owning class that provides views into the segments of the multipolygon range + /// Can only be constructed on host. + auto _segments(rmm::cuda_stream_view); /// Range Casting @@ -205,10 +204,6 @@ class multipolygon_range { VecIterator _point_begin; VecIterator _point_end; - // TODO: find a better name - CUSPATIAL_HOST_DEVICE auto subtracted_ring_begin(); - CUSPATIAL_HOST_DEVICE auto subtracted_ring_end(); - private: template CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 ring_idx); diff --git a/cpp/include/cuspatial_test/vector_equality.hpp b/cpp/include/cuspatial_test/vector_equality.hpp index 5d72f6228..4abde035d 100644 --- a/cpp/include/cuspatial_test/vector_equality.hpp +++ b/cpp/include/cuspatial_test/vector_equality.hpp @@ -37,7 +37,7 @@ namespace cuspatial { namespace test { /** - * @brief Compare two floats are close within N ULPs + * @brief Compare two floats are close within N ULPs, nans are treated equal * * N is predefined by GoogleTest * https://google.github.io/googletest/reference/assertions.html#EXPECT_FLOAT_EQ @@ -46,22 +46,22 @@ template auto floating_eq_by_ulp(T val) { if constexpr (std::is_same_v) { - return ::testing::FloatEq(val); + return ::testing::NanSensitiveFloatEq(val); } else { - return ::testing::DoubleEq(val); + return ::testing::NanSensitiveDoubleEq(val); } } /** - * @brief Compare two floats are close within `abs_error` + * @brief Compare two floats are close within `abs_error`, nans are treated equal */ template auto floating_eq_by_abs_error(T val, T abs_error) { if constexpr (std::is_same_v) { - return ::testing::FloatNear(val, abs_error); + return ::testing::NanSensitiveFloatNear(val, abs_error); } else { - return ::testing::DoubleNear(val, abs_error); + return ::testing::NanSensitiveDoubleNear(val, abs_error); } } diff --git a/cpp/tests/distance/linestring_polygon_distance_test.cu b/cpp/tests/distance/linestring_polygon_distance_test.cu index bf74b731b..7cc3b6b55 100644 --- a/cpp/tests/distance/linestring_polygon_distance_test.cu +++ b/cpp/tests/distance/linestring_polygon_distance_test.cu @@ -413,3 +413,108 @@ TYPED_TEST(PairwiseLinestringPolygonDistanceTest, TwoPairsCrosses) P{5, 5}}, {std::sqrt(T{2}), 0.0}); } + +// Empty Geometries Tests + +/// Empty MultiLinestring vs Non-empty multipolygons +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, ThreePairEmptyMultiLinestring) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + + {0, 1, 2, 3}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-20, -20}, + P{-20, -21}, + P{-21, -21}, + P{-21, -20}, + P{-20, -20}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} + +/// Non-empty MultiLinestring vs Empty multipolygons +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, ThreePairEmptyMultiPolygon) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 2, 3}, + {0, 4, 7, 10}, + {P{0, 0}, + P{1, 1}, + P{2, 2}, + P{3, 3}, + P{20, 20}, + P{21, 21}, + P{22, 22}, + P{10, 10}, + P{11, 11}, + P{12, 12}}, + + {0, 1, 1, 2}, + {0, 1, 2}, + {0, 4, 9}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} + +/// FIXME: Empty MultiLinestring vs Empty multipolygons +/// This example fails at distance util, where point-polyogn intersection kernel doesn't handle +/// empty multipoint/multipolygons. +TYPED_TEST(PairwiseLinestringPolygonDistanceTest, + DISABLED_ThreePairEmptyMultiLineStringEmptyMultiPolygon) +{ + using T = TypeParam; + using P = vec_2d; + + CUSPATIAL_RUN_TEST(this->run_single, + + {0, 1, 1, 3}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + + {0, 1, 1, 2}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {P{-1, -1}, + P{-2, -2}, + P{-2, -1}, + P{-1, -1}, + + P{-10, -10}, + P{-10, -11}, + P{-11, -11}, + P{-11, -10}, + P{-10, -10}}, + {std::sqrt(T{2}), std::numeric_limits::quiet_NaN(), 20 * std::sqrt(T{2})}); +} diff --git a/cpp/tests/range/multilinestring_range_test.cu b/cpp/tests/range/multilinestring_range_test.cu index ec374ef8f..4d891ac06 100644 --- a/cpp/tests/range/multilinestring_range_test.cu +++ b/cpp/tests/range/multilinestring_range_test.cu @@ -13,6 +13,7 @@ * See the License for the specific language governing permissions and * limitations under the License. */ +#include #include #include @@ -38,11 +39,13 @@ struct MultilinestringRangeTest : public BaseFixture { { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); - auto rng = multilinestring_array.range(); - - rmm::device_uvector> got(rng.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), rng.segment_begin(), rng.segment_end(), got.begin()); + rmm::device_uvector> got(segments_range.num_segments(), stream()); + thrust::copy( + rmm::exec_policy(stream()), segments_range.begin(), segments_range.end(), got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); CUSPATIAL_EXPECT_VEC2D_PAIRS_EQUIVALENT(d_expected, got); @@ -76,15 +79,39 @@ struct MultilinestringRangeTest : public BaseFixture { { auto multilinestring_array = make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); auto d_expected = thrust::device_vector(expected.begin(), expected.end()); - auto rng = multilinestring_array.range(); - rmm::device_uvector got(rng.num_multilinestrings(), stream()); thrust::copy(rmm::exec_policy(stream()), - rng.multilinestring_segment_count_begin(), - rng.multilinestring_segment_count_end(), + segments_range.multigeometry_count_begin(), + segments_range.multigeometry_count_end(), + got.begin()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); + } + + void run_multilinestring_segment_offset_test(std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list> coordinates, + std::initializer_list expected) + { + auto multilinestring_array = + make_multilinestring_array(geometry_offset, part_offset, coordinates); + auto rng = multilinestring_array.range(); + auto segments = rng._segments(stream()); + auto segments_range = segments.view(); + + auto d_expected = thrust::device_vector(expected.begin(), expected.end()); + + rmm::device_uvector got(rng.num_multilinestrings() + 1, stream()); + + thrust::copy(rmm::exec_policy(stream()), + segments_range.multigeometry_offset_begin(), + segments_range.multigeometry_offset_end(), got.begin()); CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); @@ -190,22 +217,26 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) CUSPATIAL_RUN_TEST(this->run_segment_test_single, {0, 1, 2, 3}, {0, 6, 11, 14}, - {P{0, 0}, - P{1, 1}, - P{2, 2}, - P{3, 3}, - P{4, 4}, - P{5, 5}, - - P{10, 10}, - P{11, 11}, - P{12, 12}, - P{13, 13}, - P{14, 14}, - - P{20, 20}, - P{21, 21}, - P{22, 22}}, + { + + P{0, 0}, + P{1, 1}, + P{2, 2}, + P{3, 3}, + P{4, 4}, + P{5, 5}, + + P{10, 10}, + P{11, 11}, + P{12, 12}, + P{13, 13}, + P{14, 14}, + + P{20, 20}, + P{21, 21}, + P{22, 22} + + }, {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, @@ -220,8 +251,7 @@ TYPED_TEST(MultilinestringRangeTest, SegmentIteratorManyPairTest) S{P{21, 21}, P{22, 22}}}); } -/// FIXME: Currently, segment iterator doesn't handle empty linestrings. -TYPED_TEST(MultilinestringRangeTest, DISABLED_SegmentIteratorWithEmptyLineTest) +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyLineTest) { using T = TypeParam; using P = vec_2d; @@ -235,6 +265,38 @@ TYPED_TEST(MultilinestringRangeTest, DISABLED_SegmentIteratorWithEmptyLineTest) {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, S{P{10, 10}, P{11, 11}}, S{P{11, 11}, P{12, 12}}}); } +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST( + this->run_segment_test_single, + {0, 1, 1, 2}, + {0, 3, 6}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, + {S{P{0, 0}, P{1, 1}}, S{P{1, 1}, P{2, 2}}, S{P{10, 10}, P{11, 11}}, S{P{11, 11}, P{12, 12}}}); +} + +TYPED_TEST(MultilinestringRangeTest, SegmentIteratorWithEmptyMultiLineStringTest2) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_segment_test_single, + + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {S{P{0, 0}, P{1, 1}}, + S{P{1, 1}, P{2, 2}}, + S{P{2, 2}, P{3, 3}}, + S{P{10, 10}, P{11, 11}}, + S{P{11, 11}, P{12, 12}}}); +} + TYPED_TEST(MultilinestringRangeTest, PerMultilinestringCountTest) { using T = TypeParam; @@ -271,8 +333,7 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest) this->run_multilinestring_segment_count_test, {0, 1}, {0, 3}, {P{0, 0}, P{1, 1}, P{2, 2}}, {2}); } -// FIXME: contains empty linestring -TYPED_TEST(MultilinestringRangeTest, DISABLED_MultilinestringSegmentCountTest2) +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest2) { using T = TypeParam; using P = vec_2d; @@ -311,8 +372,8 @@ TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest4) {2, 2}); } -// FIXME: contains empty linestring -TYPED_TEST(MultilinestringRangeTest, DISABLED_MultilinestringSegmentCountTest5) +// contains empty linestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest5) { using T = TypeParam; using P = vec_2d; @@ -325,6 +386,48 @@ TYPED_TEST(MultilinestringRangeTest, DISABLED_MultilinestringSegmentCountTest5) {2, 0, 2}); } +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest6) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + {0, 1, 1, 2}, + {0, 3, 6}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{10, 10}, P{11, 11}, P{12, 12}}, + {2, 0, 2}); +} + +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentCountTest7) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_count_test, + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {3, 0, 2}); +} + +// contains empty multilinestring +TYPED_TEST(MultilinestringRangeTest, MultilinestringSegmentOffsetTest) +{ + using T = TypeParam; + using P = vec_2d; + using S = segment; + + CUSPATIAL_RUN_TEST(this->run_multilinestring_segment_offset_test, + {0, 1, 1, 2}, + {0, 4, 7}, + {P{0, 0}, P{1, 1}, P{2, 2}, P{3, 3}, P{10, 10}, P{11, 11}, P{12, 12}}, + {0, 3, 3, 5}); +} + TYPED_TEST(MultilinestringRangeTest, MultilinestringLinestringCountTest) { using T = TypeParam; diff --git a/cpp/tests/range/multipolygon_range_test.cu b/cpp/tests/range/multipolygon_range_test.cu index 90ffa7e18..fc0afffc7 100644 --- a/cpp/tests/range/multipolygon_range_test.cu +++ b/cpp/tests/range/multipolygon_range_test.cu @@ -33,19 +33,23 @@ using namespace cuspatial::test; template struct MultipolygonRangeTest : public BaseFixture { - void run_multipolygon_segment_iterator_single(std::initializer_list geometry_offset, - std::initializer_list part_offset, - std::initializer_list ring_offset, - std::initializer_list> coordinates, - std::initializer_list> expected) + void run_multipolygon_segment_method_iterator_single( + std::initializer_list geometry_offset, + std::initializer_list part_offset, + std::initializer_list ring_offset, + std::initializer_list> coordinates, + std::initializer_list> expected) { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); + auto rng = multipolygon_array.range(); + auto segments = rng._segments(stream()); + auto segment_range = segments.view(); - auto got = rmm::device_uvector>(rng.num_segments(), stream()); + auto got = rmm::device_uvector>(segment_range.num_segments(), stream()); - thrust::copy(rmm::exec_policy(stream()), rng.segment_begin(), rng.segment_end(), got.begin()); + thrust::copy( + rmm::exec_policy(stream()), segment_range.begin(), segment_range.end(), got.begin()); auto d_expected = thrust::device_vector>(expected.begin(), expected.end()); @@ -76,7 +80,7 @@ struct MultipolygonRangeTest : public BaseFixture { CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, d_expected); } - void run_multipolygon_segment_count_single( + void run_multipolygon_segment_method_count_single( std::initializer_list geometry_offset, std::initializer_list part_offset, std::initializer_list ring_offset, @@ -85,13 +89,15 @@ struct MultipolygonRangeTest : public BaseFixture { { auto multipolygon_array = make_multipolygon_array(geometry_offset, part_offset, ring_offset, coordinates); - auto rng = multipolygon_array.range(); + auto rng = multipolygon_array.range(); + auto segments = rng._segments(stream()); + auto segment_range = segments.view(); auto got = rmm::device_uvector(rng.num_multipolygons(), stream()); thrust::copy(rmm::exec_policy(stream()), - rng.multipolygon_segment_count_begin(), - rng.multipolygon_segment_count_end(), + segment_range.multigeometry_count_begin(), + segment_range.multigeometry_count_end(), got.begin()); auto d_expected = thrust::device_vector(expected_segment_counts.begin(), @@ -158,17 +164,17 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators) using T = TypeParam; using P = vec_2d; using S = segment; - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1}, {0, 1}, {0, 4}, {{0, 0}, {1, 0}, {1, 1}, {0, 0}}, - {S{P{0, 0}, P{1, 0}}, S{P{1, 0}, P{1, 1}}, S{P{1, 1}, P{0, 0}}}); + {S{{0, 0}, P{1, 0}}, S{P{1, 0}, P{1, 1}}, S{P{1, 1}, P{0, 0}}}); } TYPED_TEST(MultipolygonRangeTest, SegmentIterators2) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1}, {0, 2}, {0, 4, 8}, @@ -183,7 +189,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators2) TYPED_TEST(MultipolygonRangeTest, SegmentIterators3) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 2}, {0, 1, 2}, {0, 4, 8}, @@ -198,7 +204,7 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators3) TYPED_TEST(MultipolygonRangeTest, SegmentIterators4) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_iterator_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, {0, 1, 2}, {0, 1, 2}, {0, 4, 8}, @@ -211,6 +217,87 @@ TYPED_TEST(MultipolygonRangeTest, SegmentIterators4) {{11, 11}, {10, 10}}}); } +TYPED_TEST(MultipolygonRangeTest, SegmentIterators5) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 2, 3}, + {0, 1, 2, 3}, + {0, 4, 9, 14}, + {{-1, -1}, + {-2, -2}, + {-2, -1}, + {-1, -1}, + + {-20, -20}, + {-20, -21}, + {-21, -21}, + {-21, -20}, + {-20, -20}, + + {-10, -10}, + {-10, -11}, + {-11, -11}, + {-11, -10}, + {-10, -10}}, + + {{{-1, -1}, {-2, -2}}, + {{-2, -2}, {-2, -1}}, + {{-2, -1}, {-1, -1}}, + {{-20, -20}, {-20, -21}}, + {{-20, -21}, {-21, -21}}, + {{-21, -21}, {-21, -20}}, + {{-21, -20}, {-20, -20}}, + {{-10, -10}, {-10, -11}}, + {{-10, -11}, {-11, -11}}, + {{-11, -11}, {-11, -10}}, + {{-11, -10}, {-10, -10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators5EmptyRing) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 2}, + {0, 1, 3}, + {0, 4, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators6EmptyPolygon) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 3}, + {0, 1, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + +TYPED_TEST(MultipolygonRangeTest, SegmentIterators7EmptyMultiPolygon) +{ + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_iterator_single, + {0, 1, 1, 2}, + {0, 1, 2}, + {0, 4, 8}, + {{0, 0}, {1, 0}, {1, 1}, {0, 0}, {10, 10}, {11, 10}, {11, 11}, {10, 10}}, + {{{0, 0}, {1, 0}}, + {{1, 0}, {1, 1}}, + {{1, 1}, {0, 0}}, + {{10, 10}, {11, 10}}, + {{11, 10}, {11, 11}}, + {{11, 11}, {10, 10}}}); +} + TYPED_TEST(MultipolygonRangeTest, MultipolygonCountIterator) { CUSPATIAL_RUN_TEST(this->run_multipolygon_point_count_iterator_single, @@ -280,7 +367,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonCountIterator4) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 1}, {0, 1}, {0, 4}, @@ -291,7 +378,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount2) { CUSPATIAL_RUN_TEST( - this->run_multipolygon_segment_count_single, + this->run_multipolygon_segment_method_count_single, {0, 1}, {0, 2}, {0, 4, 8}, @@ -301,7 +388,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount2) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount3) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2}, {0, 2, 3}, {0, 4, 8, 12}, @@ -322,7 +409,7 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount3) TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount4) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2, 3}, {0, 2, 3, 4}, {0, 4, 8, 12, 16}, @@ -345,16 +432,19 @@ TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount4) {9, 3}); } -// FIXME: multipolygon doesn't constructor doesn't allow empty rings, should it? -TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyRing) +// FIXME: multipolygon constructor doesn't allow empty rings, should it? +TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount_ContainsEmptyRing) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 2, 3}, {0, 2, 3, 4}, - {0, 4, 4, 8, 12}, + {0, 7, 7, 11, 18}, {{0, 0}, {1, 0}, {1, 1}, + {0.5, 1.5}, + {0, 1.0}, + {0.5, 0.5}, {0, 0}, {0.2, 0.2}, {0.2, 0.3}, @@ -363,16 +453,19 @@ TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmp {0, 0}, {1, 0}, {1, 1}, + {0.5, 1.5}, + {0, 1.0}, + {0.5, 0.5}, {0, 0}}, - {6, 3}); + {9, 6}); } -// FIXME: multipolygon doesn't constructor doesn't allow empty rings, should it? -TYPED_TEST(MultipolygonRangeTest, DISABLED_MultipolygonSegmentCount_ConatainsEmptyPart) +// FIXME: multipolygon constructor doesn't allow empty rings, should it? +TYPED_TEST(MultipolygonRangeTest, MultipolygonSegmentCount_ContainsEmptyPart) { - CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_count_single, + CUSPATIAL_RUN_TEST(this->run_multipolygon_segment_method_count_single, {0, 3, 4}, - {0, 2, 2, 3, 4}, + {0, 1, 1, 2, 3}, {0, 4, 8, 12}, {{0, 0}, {1, 0},