diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 7a5b8acae..b13d53cfc 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -136,6 +136,7 @@ add_library(cuspatial src/spatial/point_linestring_distance.cu src/spatial/point_polygon_distance.cu src/spatial/linestring_polygon_distance.cu + src/spatial/polygon_distance.cu src/spatial/point_linestring_nearest_points.cu src/spatial/sinusoidal_projection.cu src/trajectory/derive_trajectories.cu diff --git a/cpp/include/cuspatial/assert.cuh b/cpp/include/cuspatial/assert.cuh new file mode 100644 index 000000000..c5e13d6f4 --- /dev/null +++ b/cpp/include/cuspatial/assert.cuh @@ -0,0 +1,35 @@ +/* + * 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 + +/** + * @brief `assert`-like macro for device code + * + * This is effectively the same as the standard `assert` macro, except it + * relies on the `__PRETTY_FUNCTION__` macro which is specific to GCC and Clang + * to produce better assert messages. + */ +#if !defined(NDEBUG) && defined(__CUDA_ARCH__) && (defined(__clang__) || defined(__GNUC__)) +#define __ASSERT_STR_HELPER(x) #x +#define cuspatial_assert(e) \ + ((e) ? static_cast(0) \ + : __assert_fail(__ASSERT_STR_HELPER(e), __FILE__, __LINE__, __PRETTY_FUNCTION__)) +#else +#define cuspatial_assert(e) (static_cast(0)) +#endif diff --git a/cpp/include/cuspatial/detail/utility/validation.hpp b/cpp/include/cuspatial/detail/utility/validation.hpp index 7fe5742ea..1a4bd2001 100644 --- a/cpp/include/cuspatial/detail/utility/validation.hpp +++ b/cpp/include/cuspatial/detail/utility/validation.hpp @@ -35,10 +35,10 @@ * [2]: https://arrow.apache.org/docs/format/Columnar.html#variable-size-binary-layout */ #define CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets) \ - CUSPATIAL_EXPECTS(num_linestring_offsets > 0, \ - "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_linestring_points >= 2 * (num_linestring_offsets - 1), \ - "Each linestring must have at least two vertices"); + CUSPATIAL_HOST_DEVICE_EXPECTS(num_linestring_offsets > 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"); /** * @brief Macro for validating the data array sizes for multilinestrings. @@ -57,10 +57,10 @@ * [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_MULTILINESTRING_SIZES( \ - num_linestring_points, num_multilinestring_offsets, num_linestring_offsets) \ - CUSPATIAL_EXPECTS(num_multilinestring_offsets > 0, \ - "Multilinestring offsets must contain at least one (1) value"); \ +#define CUSPATIAL_EXPECTS_VALID_MULTILINESTRING_SIZES( \ + num_linestring_points, num_multilinestring_offsets, num_linestring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_multilinestring_offsets > 0, \ + "Multilinestring offsets must contain at least one (1) value"); \ CUSPATIAL_EXPECTS_VALID_LINESTRING_SIZES(num_linestring_points, num_linestring_offsets); /** @@ -84,15 +84,16 @@ * [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_EXPECTS(num_poly_offsets > 0, "Polygon offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_poly_ring_offsets > 0, \ - "Polygon ring offsets must contain at least one (1) value"); \ - CUSPATIAL_EXPECTS(num_poly_ring_offsets >= num_poly_offsets, \ - "Each polygon must have at least one (1) ring"); \ - CUSPATIAL_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"); \ + 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"); /** * @brief Macro for validating the data array sizes for a multipolygon. @@ -116,8 +117,8 @@ * [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_MULTIPOLYGON_SIZES( \ - num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets) \ - CUSPATIAL_EXPECTS(num_multipoly_offsets > 0, \ - "Multipolygon offsets must contain at least one (1) value"); \ +#define CUSPATIAL_EXPECTS_VALID_MULTIPOLYGON_SIZES( \ + num_poly_points, num_multipoly_offsets, num_poly_offsets, num_poly_ring_offsets) \ + CUSPATIAL_HOST_DEVICE_EXPECTS(num_multipoly_offsets > 0, \ + "Multipolygon offsets must contain at least one (1) value"); \ CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES(num_poly_points, num_poly_offsets, num_poly_ring_offsets); diff --git a/cpp/include/cuspatial/distance/polygon_distance.hpp b/cpp/include/cuspatial/distance/polygon_distance.hpp new file mode 100644 index 000000000..4b1a291c9 --- /dev/null +++ b/cpp/include/cuspatial/distance/polygon_distance.hpp @@ -0,0 +1,44 @@ +/* + * 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 + +namespace cuspatial { + +/** + * @brief Compute pairwise (multi)polygon-to-(multi)polygon Cartesian distance + * + * Computes the cartesian distance between each pair of the multipolygons. + * + * @param lhs Geometry column of the multipolygons to compute distance from + * @param rhs Geometry column of the multipolygons to compute distance to + * @param mr Device memory resource used to allocate the returned column. + * + * @return Column of distances between each pair of input geometries, same type as input coordinate + * types. + */ +std::unique_ptr pairwise_polygon_distance( + geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::mr::device_memory_resource* mr = rmm::mr::get_current_device_resource()); + +} // namespace cuspatial diff --git a/cpp/include/cuspatial/error.hpp b/cpp/include/cuspatial/error.hpp index b1190c8ac..d2973349f 100644 --- a/cpp/include/cuspatial/error.hpp +++ b/cpp/include/cuspatial/error.hpp @@ -1,5 +1,5 @@ /* - * Copyright (c) 2020-2022, NVIDIA CORPORATION. + * Copyright (c) 2020-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. @@ -16,6 +16,8 @@ #pragma once +#include + #include #include @@ -76,6 +78,31 @@ struct cuda_error : public std::runtime_error { : throw cuspatial::logic_error("cuSpatial failure at: " __FILE__ \ ":" CUSPATIAL_STRINGIFY(__LINE__) ": " reason) +/**---------------------------------------------------------------------------* + * @brief Macro for checking (pre-)conditions that throws an exception when + * a condition is violated. + * + * Example usage: + * + * @code + * CUSPATIAL_HOST_DEVICE_EXPECTS(lhs->dtype == rhs->dtype, "Column type mismatch"); + * @endcode + * + * @param[in] cond Expression that evaluates to true or false + * @param[in] reason String literal description of the reason that cond is + * expected to be true + * + * (if on host) + * @throw cuspatial::logic_error if the condition evaluates to false. + * (if on device) + * program terminates and assertion error message is printed to stderr. + *---------------------------------------------------------------------------**/ +#ifndef __CUDA_ARCH__ +#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) CUSPATIAL_EXPECTS(cond, reason) +#else +#define CUSPATIAL_HOST_DEVICE_EXPECTS(cond, reason) cuspatial_assert(cond&& reason) +#endif + /**---------------------------------------------------------------------------* * @brief Indicates that an erroneous code path has been taken. * diff --git a/cpp/include/cuspatial/range/multilinestring_range.cuh b/cpp/include/cuspatial/range/multilinestring_range.cuh index 04c280e04..18a30a2bf 100644 --- a/cpp/include/cuspatial/range/multilinestring_range.cuh +++ b/cpp/include/cuspatial/range/multilinestring_range.cuh @@ -59,12 +59,12 @@ class multilinestring_range { using point_t = iterator_value_type; using element_t = iterator_vec_base_type; - multilinestring_range(GeometryIterator geometry_begin, - GeometryIterator geometry_end, - PartIterator part_begin, - PartIterator part_end, - VecIterator points_begin, - VecIterator points_end); + CUSPATIAL_HOST_DEVICE multilinestring_range(GeometryIterator geometry_begin, + GeometryIterator geometry_end, + PartIterator part_begin, + PartIterator part_end, + VecIterator points_begin, + VecIterator points_end); /// Return the number of multilinestrings in the array. CUSPATIAL_HOST_DEVICE auto size() { return num_multilinestrings(); } diff --git a/cpp/src/spatial/polygon_distance.cu b/cpp/src/spatial/polygon_distance.cu new file mode 100644 index 000000000..8d3192e13 --- /dev/null +++ b/cpp/src/spatial/polygon_distance.cu @@ -0,0 +1,113 @@ +/* + * 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 "../utility/iterator.hpp" +#include "../utility/multi_geometry_dispatch.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace cuspatial { + +namespace detail { + +namespace { + +template +struct pairwise_polygon_distance_impl { + using SizeType = cudf::device_span::size_type; + + template )> + std::unique_ptr operator()(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + auto lhs_range = make_multipolygon_range(lhs); + auto rhs_range = make_multipolygon_range(rhs); + + auto output = cudf::make_numeric_column( + lhs.coordinate_type(), lhs.size(), cudf::mask_state::UNALLOCATED, stream, mr); + + cuspatial::pairwise_polygon_distance( + lhs_range, rhs_range, output->mutable_view().begin(), stream); + return output; + } + + template ), typename... Args> + std::unique_ptr operator()(Args&&...) + + { + CUSPATIAL_FAIL("polygon distance API only supports floating point coordinates."); + } +}; + +} // namespace + +template +struct pairwise_polygon_distance { + std::unique_ptr operator()(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::cuda_stream_view stream, + rmm::mr::device_memory_resource* mr) + { + return cudf::type_dispatcher( + lhs.coordinate_type(), + pairwise_polygon_distance_impl{}, + lhs, + rhs, + stream, + mr); + } +}; + +} // namespace detail + +std::unique_ptr pairwise_polygon_distance(geometry_column_view const& lhs, + geometry_column_view const& rhs, + rmm::mr::device_memory_resource* mr) +{ + CUSPATIAL_EXPECTS(lhs.geometry_type() == geometry_type_id::POLYGON && + rhs.geometry_type() == geometry_type_id::POLYGON, + "Unexpected input geometry types."); + + CUSPATIAL_EXPECTS(lhs.coordinate_type() == rhs.coordinate_type(), + "Input geometries must have the same coordinate data types."); + + return multi_geometry_double_dispatch( + lhs.collection_type(), rhs.collection_type(), lhs, rhs, rmm::cuda_stream_default, mr); +} + +} // namespace cuspatial diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 3ec61c6ce..4f9533537 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -89,12 +89,15 @@ ConfigureTest(POINT_LINESTRING_DISTANCE_TEST ConfigureTest(LINESTRING_DISTANCE_TEST spatial/distance/linestring_distance_test.cpp) -ConfigureTest(LINESTRING_POLYGON_DISTANCE_TEST - spatial/linestring_polygon_distance_test.cpp) - ConfigureTest(POINT_POLYGON_DISTANCE_TEST spatial/distance/point_polygon_distance_test.cpp) +ConfigureTest(LINESTRING_POLYGON_DISTANCE_TEST + spatial/distance/linestring_polygon_distance_test.cpp) + +ConfigureTest(POLYGON_DISTANCE_TEST + spatial/distance/polygon_distance_test.cpp) + # equality ConfigureTest(PAIRWISE_MULTIPOINT_EQUALS_COUNT_TEST spatial/equality/pairwise_multipoint_equals_count_test.cpp) diff --git a/cpp/tests/spatial/linestring_polygon_distance_test.cpp b/cpp/tests/spatial/distance/linestring_polygon_distance_test.cpp similarity index 100% rename from cpp/tests/spatial/linestring_polygon_distance_test.cpp rename to cpp/tests/spatial/distance/linestring_polygon_distance_test.cpp diff --git a/cpp/tests/spatial/distance/polygon_distance_test.cpp b/cpp/tests/spatial/distance/polygon_distance_test.cpp new file mode 100644 index 000000000..1bea37f43 --- /dev/null +++ b/cpp/tests/spatial/distance/polygon_distance_test.cpp @@ -0,0 +1,149 @@ +/* + * 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 + +#include +#include + +#include + +#include + +#include +#include + +using namespace cuspatial; +using namespace cuspatial::test; + +using namespace cudf; +using namespace cudf::test; + +template +struct PairwisePolygonDistanceTestBase : ::testing::Test { + void run_single(geometry_column_view lhs, + geometry_column_view rhs, + std::initializer_list expected) + { + auto got = pairwise_polygon_distance(lhs, rhs); + CUDF_TEST_EXPECT_COLUMNS_EQUIVALENT(*got, fixed_width_column_wrapper(expected)); + } + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } +}; + +template +struct PairwisePolygonDistanceTestEmpty : PairwisePolygonDistanceTestBase { + void SetUp() + { + [[maybe_unused]] collection_type_id _; + std::tie(_, empty_polygon_column) = make_polygon_column({0}, {0}, {}, this->stream()); + std::tie(_, empty_multipolygon_column) = + make_polygon_column({0}, {0}, {0}, {}, this->stream()); + } + + geometry_column_view empty_polygon() + { + return geometry_column_view( + empty_polygon_column->view(), collection_type_id::SINGLE, geometry_type_id::POLYGON); + } + + geometry_column_view empty_multipolygon() + { + return geometry_column_view( + empty_multipolygon_column->view(), collection_type_id::MULTI, geometry_type_id::POLYGON); + } + + std::unique_ptr empty_polygon_column; + std::unique_ptr empty_multipolygon_column; +}; + +using TestTypes = ::testing::Types; +TYPED_TEST_CASE(PairwisePolygonDistanceTestEmpty, TestTypes); + +struct PairwisePolygonDistanceTestUntyped : testing::Test { + rmm::cuda_stream_view stream() { return cudf::get_default_stream(); } +}; +TYPED_TEST(PairwisePolygonDistanceTestEmpty, SingleToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_polygon(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, SingleToMultiEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_polygon(), this->empty_multipolygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, MultiToSingleEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_multipolygon(), this->empty_polygon(), {}); +}; + +TYPED_TEST(PairwisePolygonDistanceTestEmpty, MultiToMultiEmpty) +{ + CUSPATIAL_RUN_TEST(this->run_single, this->empty_multipolygon(), this->empty_multipolygon(), {}); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, SizeMismatch) +{ + auto [ptype, polygons1] = make_polygon_column( + {0, 1, 2}, {0, 4, 8}, {0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0}, this->stream()); + + auto [polytype, polygons2] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto polygons1_view = + geometry_column_view(polygons1->view(), ptype, geometry_type_id::LINESTRING); + auto polygons2_view = + geometry_column_view(polygons1->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(polygons1_view, polygons2_view), cuspatial::logic_error); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, TypeMismatch) +{ + auto [ptype, polygons1] = make_polygon_column( + {0, 1, 2}, {0, 4, 8}, {0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0}, this->stream()); + + auto [polytype, polygons2] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto polygons1_view = + geometry_column_view(polygons1->view(), ptype, geometry_type_id::LINESTRING); + auto polygons2_view = + geometry_column_view(polygons2->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(polygons1_view, polygons2_view), cuspatial::logic_error); +}; + +TEST_F(PairwisePolygonDistanceTestUntyped, WrongGeometryType) +{ + auto [ptype, points] = make_point_column({0, 1}, {0.0, 0.0}, this->stream()); + + auto [polytype, polygons] = + make_polygon_column({0, 1}, {0, 1}, {0, 4}, {1, 1, 1, 2, 2, 2, 1, 1}, this->stream()); + + auto points_view = geometry_column_view(points->view(), ptype, geometry_type_id::POINT); + auto polygons_view = geometry_column_view(polygons->view(), polytype, geometry_type_id::POLYGON); + + EXPECT_THROW(pairwise_polygon_distance(points_view, polygons_view), cuspatial::logic_error); +};