From 9d91cd3063acc49a3d86e8bb49640d5958c6b615 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Tue, 20 Jun 2023 13:13:41 -0400 Subject: [PATCH 1/3] Updates Build Instructions to Adjust for Devcontainer Instructions and Remove Stale Infomation (#1179) closes #1174 Minor improvement: - updates build development environment to point to devcontainers, not rapids-compose. - Move point-linestring python test to correct location. Authors: - Michael Wang (https://github.com/isVoid) Approvers: - Ben Jarmak (https://github.com/jarmak-nv) - Mark Harris (https://github.com/harrism) URL: https://github.com/rapidsai/cuspatial/pull/1179 --- docs/source/developer_guide/build.md | 60 ++++++++++++------ .../development_environment.md | 36 +++++++---- .../test_point_linestring_distance.py | 0 test_fixtures/shapefiles/empty_poly.cpg | 1 - test_fixtures/shapefiles/empty_poly.dbf | Bin 66 -> 0 bytes test_fixtures/shapefiles/empty_poly.sbn | Bin 108 -> 0 bytes test_fixtures/shapefiles/empty_poly.sbx | Bin 108 -> 0 bytes test_fixtures/shapefiles/empty_poly.shp | Bin 100 -> 0 bytes test_fixtures/shapefiles/empty_poly.shx | Bin 100 -> 0 bytes test_fixtures/shapefiles/one_poly.cpg | 1 - test_fixtures/shapefiles/one_poly.dbf | Bin 86 -> 0 bytes test_fixtures/shapefiles/one_poly.sbn | Bin 132 -> 0 bytes test_fixtures/shapefiles/one_poly.sbx | Bin 116 -> 0 bytes test_fixtures/shapefiles/one_poly.shp | Bin 236 -> 0 bytes test_fixtures/shapefiles/one_poly.shx | Bin 108 -> 0 bytes test_fixtures/shapefiles/two_polys.dbf | Bin 90 -> 0 bytes test_fixtures/shapefiles/two_polys.sbn | Bin 164 -> 0 bytes test_fixtures/shapefiles/two_polys.sbx | Bin 124 -> 0 bytes test_fixtures/shapefiles/two_polys.shp | Bin 372 -> 0 bytes test_fixtures/shapefiles/two_polys.shx | Bin 116 -> 0 bytes 20 files changed, 63 insertions(+), 35 deletions(-) rename python/cuspatial/cuspatial/tests/{ => spatial/distance}/test_point_linestring_distance.py (100%) delete mode 100644 test_fixtures/shapefiles/empty_poly.cpg delete mode 100644 test_fixtures/shapefiles/empty_poly.dbf delete mode 100644 test_fixtures/shapefiles/empty_poly.sbn delete mode 100644 test_fixtures/shapefiles/empty_poly.sbx delete mode 100644 test_fixtures/shapefiles/empty_poly.shp delete mode 100644 test_fixtures/shapefiles/empty_poly.shx delete mode 100644 test_fixtures/shapefiles/one_poly.cpg delete mode 100644 test_fixtures/shapefiles/one_poly.dbf delete mode 100644 test_fixtures/shapefiles/one_poly.sbn delete mode 100644 test_fixtures/shapefiles/one_poly.sbx delete mode 100644 test_fixtures/shapefiles/one_poly.shp delete mode 100644 test_fixtures/shapefiles/one_poly.shx delete mode 100644 test_fixtures/shapefiles/two_polys.dbf delete mode 100644 test_fixtures/shapefiles/two_polys.sbn delete mode 100644 test_fixtures/shapefiles/two_polys.sbx delete mode 100644 test_fixtures/shapefiles/two_polys.shp delete mode 100644 test_fixtures/shapefiles/two_polys.shx diff --git a/docs/source/developer_guide/build.md b/docs/source/developer_guide/build.md index ea6da6bc3..c7d56d41d 100644 --- a/docs/source/developer_guide/build.md +++ b/docs/source/developer_guide/build.md @@ -18,30 +18,50 @@ git clone https://github.com/rapidsai/cuspatial.git $CUSPATIAL_HOME 2. clone the cuSpatial repo ```shell -conda env update --file conda/environments/all_cuda-118_arch-x86_64.yaml +conda env create -n cuspatial --file conda/environments/all_cuda-118_arch-x86_64.yaml ``` -## Build and install cuSpatial +## Build cuSpatial -1. Compile and install - ```shell - cd $CUSPATIAL_HOME && \ - chmod +x ./build.sh && \ - ./build.sh - ``` +### From the cuSpatial Dev Container: -2. Run C++/Python test code +Execute `build-cuspatial-cpp to build `libcuspatial`. The following options may be added. + - `-DBUILD_TESTS=ON`: build `libcuspatial` unit tests. + - `-DBUILD_BENCHMARKS=ON`: build `libcuspatial` benchmarks. + - `-DCMAKE_BUILD_TYPE=Debug`: Create a Debug build of `libcuspatial` (default is Release). +In addition, `build-cuspatial-python` to build cuspatial cython components. - Some tests using inline data can be run directly, e.g.: +### From Bare Metal: - ```shell - $CUSPATIAL_HOME/cpp/build/gtests/LEGACY_HAUSDORFF_TEST - $CUSPATIAL_HOME/cpp/build/gtests/POINT_IN_POLYGON_TEST - python python/cuspatial/cuspatial/tests/legacy/test_hausdorff_distance.py - python python/cuspatial/cuspatial/tests/test_pip.py - ``` +Compile libcuspatial (C++), cuspatial (cython) and C++ tests: +```shell +cd $CUSPATIAL_HOME && \ +chmod +x ./build.sh && \ +./build.sh libcuspatial cuspatial tests +``` +Additionally, the following options are also commonly used: +- `benchmarks`: build libcuspatial benchmarks +- `clean`: remove all existing build artifacts and configuration +Execute `./build.sh -h` for full list of available options. + +## Validate Installation with C++ and Python Tests + +```{note} +To manage difference between branches and build types, the build directories are located at +`$CUSPATIAL_HOME/cpp/build/[release|debug]` depending on build type, and `$CUSPATIAL_HOME/cpp/build/latest`. +is a symbolic link to the most recent build directory. On bare metal builds, remove the extra `latest` level in +the path below. +``` + +- C++ tests are located within the `$CUSPATIAL_HOME/cpp/build/latest/gtests` directory. +- Python tests are located within the `$CUSPATIAL_HOME/python/cuspatial/cuspatial/tests` directory. - Some other tests involve I/O from data files under `$CUSPATIAL_HOME/test_fixtures`. - For example, `$CUSPATIAL_HOME/cpp/build/gtests/SHAPEFILE_READER_TEST` requires three - pre-generated polygon shapefiles that contain 0, 1 and 2 polygons, respectively. They are available at - `$CUSPATIAL_HOME/test_fixtures/shapefiles`
+Execute C++ tests: +```shell +ninja -C $CUSPATIAL_HOME/cpp/build/latest test +``` + +Execute Python tests: +```shell +pytest $CUSPATIAL_HOME/python/cuspatial/cuspatial/tests/ +``` diff --git a/docs/source/developer_guide/development_environment.md b/docs/source/developer_guide/development_environment.md index abac8ea8a..0a0e47992 100644 --- a/docs/source/developer_guide/development_environment.md +++ b/docs/source/developer_guide/development_environment.md @@ -1,15 +1,25 @@ # Creating a Development Environment -cuSpatial follows the RAPIDS release schedule, so developers are encouraged to develop -using the latest development branch of RAPIDS libraries that cuspatial depends on. Other -cuspatial dependencies can be found in `conda/environments/`. - -Maintaining a local development environment can be an arduous task, especially after each -RAPIDS release. Most cuspatial developers today use -[rapids-compose](https://github.com/trxcllnt/rapids-compose) to setup their development environment. -It contains helpful scripts to build a RAPIDS development container image with the required -dependencies and RAPIDS libraries automatically fetched and correctly versioned. It also provides -script commands for simple building and testing of all RAPIDS libraries, including cuSpatial. -`rapids-compose` is the recommended way to set up your environment to develop for cuspatial. - -For developers who would like to build from conda or from source, see the [build page](https://docs.rapids.ai/api/cuspatial/stable/developer_guide/build.html). +cuSpatial recommends using [Dev Containers](https://containers.dev/) to setup the development environment. +To setup Dev Containers for cuspatial, please refer to [documentation](https://github.com/rapidsai/cuspatial/tree/main/.devcontainer). + +## From Bare Metal + +RAPIDS keeps a single source of truth for library dependencies in `dependencies.yaml`. This file divides +the dependencies into several dimensions: building, testing, documentations, notebooks etc. As a developer, +you generally want to generate an environment recipe that includes everything that the library *may* use. + +To do so, install the rapids-dependency-file-generator via pip: +```shell +pip install rapids-dependency-file-generator +``` + +And run under the repo root: +```shell +rapids-dependency-file-generator --clean +``` + +The environment recipe is generated within the `conda/environments` directory. To continue the next step of building, +see the [build page](https://docs.rapids.ai/api/cuspatial/stable/developer_guide/build.html). + +For more information about how RAPIDS manages dependencies, see [README of rapids-dependency-file-generator repo](https://github.com/rapidsai/dependency-file-generator). diff --git a/python/cuspatial/cuspatial/tests/test_point_linestring_distance.py b/python/cuspatial/cuspatial/tests/spatial/distance/test_point_linestring_distance.py similarity index 100% rename from python/cuspatial/cuspatial/tests/test_point_linestring_distance.py rename to python/cuspatial/cuspatial/tests/spatial/distance/test_point_linestring_distance.py diff --git a/test_fixtures/shapefiles/empty_poly.cpg b/test_fixtures/shapefiles/empty_poly.cpg deleted file mode 100644 index 3ad133c04..000000000 --- a/test_fixtures/shapefiles/empty_poly.cpg +++ /dev/null @@ -1 +0,0 @@ -UTF-8 \ No newline at end of file diff --git a/test_fixtures/shapefiles/empty_poly.dbf b/test_fixtures/shapefiles/empty_poly.dbf deleted file mode 100644 index 857052728e1b88ef889c11929fd6749c3d9dc7f6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 66 jcmZQB=VoC50!IcB5QPEUJYC`qA);;|N|+l}39l3YYZwB} diff --git a/test_fixtures/shapefiles/empty_poly.sbn b/test_fixtures/shapefiles/empty_poly.sbn deleted file mode 100644 index 79a5619d86913d9d7eaf9d660a8be93fb8b9b547..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 108 jcmZQzQ0Myp|6c(ECNKjD)&BrXFyf*ywP6)u1c?Fw`F;tt diff --git a/test_fixtures/shapefiles/empty_poly.sbx b/test_fixtures/shapefiles/empty_poly.sbx deleted file mode 100644 index f700f0577969bd4799b832197e88c5914758eb9e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 108 jcmZQzQ0Myp|6c(ECNKjD)&BrXFyf*ywP6)80*L|u`f~~9 diff --git a/test_fixtures/shapefiles/empty_poly.shp b/test_fixtures/shapefiles/empty_poly.shp deleted file mode 100644 index 25194a3bb47bef31d2f713f52cfc698e5f38740a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 100 hcmZQzQ0HR64vbzfGcd4XmjjA^*bk9{(Kr<{0084X1hN1C diff --git a/test_fixtures/shapefiles/empty_poly.shx b/test_fixtures/shapefiles/empty_poly.shx deleted file mode 100644 index 25194a3bb47bef31d2f713f52cfc698e5f38740a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 100 hcmZQzQ0HR64vbzfGcd4XmjjA^*bk9{(Kr<{0084X1hN1C diff --git a/test_fixtures/shapefiles/one_poly.cpg b/test_fixtures/shapefiles/one_poly.cpg deleted file mode 100644 index 3ad133c04..000000000 --- a/test_fixtures/shapefiles/one_poly.cpg +++ /dev/null @@ -1 +0,0 @@ -UTF-8 \ No newline at end of file diff --git a/test_fixtures/shapefiles/one_poly.dbf b/test_fixtures/shapefiles/one_poly.dbf deleted file mode 100644 index 89544e2b6fe51d79a60de28bd79645915f7badf0..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 86 wcmZQB=VoPOU|?`$5CM{yz|GSo-Vh?}2BL(yQPuD&C@2`{86ZHawt;~Z0LknIjQ{`u diff --git a/test_fixtures/shapefiles/one_poly.sbn b/test_fixtures/shapefiles/one_poly.sbn deleted file mode 100644 index 652a33e277aab888fcfec4a4ad284f5e844d037c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 132 zcmZQzQ0Myp|6c(ECU61@F&-Y=m%nO!+21dR{;SG5KIL_DFDas1*8A~ diff --git a/test_fixtures/shapefiles/two_polys.sbn b/test_fixtures/shapefiles/two_polys.sbn deleted file mode 100644 index 448266c9872379b266c57093010fb3c81b920b82..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 164 zcmZQzQ0Myp|6c(ECI|uwF`)}Ps6a$uG`c(vJ|j?^2Z%xDG5!Yvs5)j4A1cnmzyKBp OsR4-t`3(&~nh^l8-3!xd05$Zsu!R Date: Thu, 22 Jun 2023 12:50:25 -0700 Subject: [PATCH 2/3] Fix a small typo in pairwise_linestring_distance (#1199) This PR fixes a small typo in pairwise_linestring_distance Authors: - Michael Wang (https://github.com/isVoid) Approvers: - Mark Harris (https://github.com/harrism) - H. Thomson Comer (https://github.com/thomcom) URL: https://github.com/rapidsai/cuspatial/pull/1199 --- cpp/include/cuspatial/distance.cuh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/include/cuspatial/distance.cuh b/cpp/include/cuspatial/distance.cuh index 6f7da67d4..6c1ebb821 100644 --- a/cpp/include/cuspatial/distance.cuh +++ b/cpp/include/cuspatial/distance.cuh @@ -235,9 +235,9 @@ OutputIt pairwise_point_polygon_distance(MultiPointRange multipoints, * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator * "LegacyRandomAccessIterator" */ -template +template OutputIt pairwise_linestring_distance(MultiLinestringRange1 multilinestrings1, - MultiLinstringRange2 multilinestrings2, + MultiLinestringRange2 multilinestrings2, OutputIt distances_first, rmm::cuda_stream_view stream = rmm::cuda_stream_default); From d5499d2082e69b715abf4b44de2809d8d7412303 Mon Sep 17 00:00:00 2001 From: Michael Wang Date: Thu, 22 Jun 2023 15:22:09 -0700 Subject: [PATCH 3/3] Simplify point-in-polygon header only APIs (#1192) closes #707 This PR simplifies `point_in_polygon` and `pairwise_point_in_polygon` API using `multipoint_range` and `multipolygon_range`. While these range methods supports a pair of multi geometry semantics, in this PR I'm not really targeting to create an API that works with multipoint-in-multipolygon semantics (such as "any point in any polygon", or "all points in any polygon", that sort of semantics). Therefore, this PR introduces `contains_only_single_geometry` method in `multipoint_range` and `multipolygon_range` that provides compile-time check for API that only want to work single-type geometry ranges. There are caveats to this introduction, see *caveats* section below. This refactor results in a net decrease in LOC, and the compatibility layer for `is_point_in_polygon` is removed. The API dries up to a point where there's no raw kernel anymore. In addition, previous `pairwise_point_in_polygon` API stores the row wise result in an int32_t per pair. Since it's only storing booleans, this PR uses `uint8_t` instead. This saves memory and increases throughputs. ### Caveats `contains_only_single_geometry` method is an over-constrained method comparing to what it's really testing. A multipoint range can have a materialized column of integer sequences and still represent a single geometry column. The reason behind this refactor is that the original point-in-polygon test explicitly requires (by the number of arguments) that only single geometry columns accepted by the API. This means developers know at compile time that the geometry used in the API are constructed this way. `contains_only_single_geometry` is a `constexpr` method and verifies developers construction at compile time as well. This maintains the developer expectation while allowing `multi*_range` to be retrofit into the modern APIs. Authors: - Michael Wang (https://github.com/isVoid) Approvers: - Mark Harris (https://github.com/harrism) URL: https://github.com/rapidsai/cuspatial/pull/1192 --- .../point_in_polygon/point_in_polygon.cu | 136 ++---- .../detail/algorithm/is_point_in_polygon.cuh | 24 - .../cuspatial/detail/point_in_polygon.cuh | 267 +++-------- .../detail/range/multipoint_range.cuh | 11 +- cpp/include/cuspatial/point_in_polygon.cuh | 145 ++---- .../cuspatial/range/multipoint_range.cuh | 9 +- .../cuspatial/range/multipolygon_range.cuh | 62 +++ .../cuspatial_test/vector_factories.cuh | 24 +- cpp/src/point_in_polygon/point_in_polygon.cu | 44 +- .../quadtree_point_in_polygon_test_large.cu | 25 +- .../pairwise_point_in_polygon_test.cpp | 2 +- .../pairwise_point_in_polygon_test.cu | 333 ++++++++------ .../point_in_polygon/point_in_polygon_test.cu | 420 +++++++----------- 13 files changed, 612 insertions(+), 890 deletions(-) diff --git a/cpp/benchmarks/point_in_polygon/point_in_polygon.cu b/cpp/benchmarks/point_in_polygon/point_in_polygon.cu index 13033f8c9..7f41dd371 100644 --- a/cpp/benchmarks/point_in_polygon/point_in_polygon.cu +++ b/cpp/benchmarks/point_in_polygon/point_in_polygon.cu @@ -17,6 +17,8 @@ #include #include +#include + #include #include @@ -30,94 +32,17 @@ using namespace cuspatial; -constexpr double PI = 3.141592653589793; auto constexpr radius = 10.0; -auto constexpr num_polygons = 31; -auto constexpr num_rings_per_polygon = 1; // only 1 ring for now - -/** - * @brief Generate a random point within a window of [minXY, maxXY] - */ -template -vec_2d random_point(vec_2d minXY, vec_2d maxXY) -{ - auto x = minXY.x + (maxXY.x - minXY.x) * rand() / static_cast(RAND_MAX); - auto y = minXY.y + (maxXY.y - minXY.y) * rand() / static_cast(RAND_MAX); - return vec_2d{x, y}; -} - -/** - * @brief Helper to generate 31 simple polygons used for benchmarks. - * - * The polygons are generated by setting a centroid and a radius. The vertices of the - * polygons are generated by rotating a circle around the centroid. The centroid of - * the polygon is randomly sampled from window [minXY, maxXY]. - * - * @tparam T The floating point type for the coordinates - * @param num_sides Number of sides of the polygon - * @param radius The radius of the circle from which the vertices are sampled - * @param minXY The minimum xy coordinates of the window from which the centroid is sampled - * @param maxXY The maximum xy coordinates of the window from which the centroid is sampled - * @return 32 polygons in structure of arrays: - * [polygon offset, poly ring offset, point coordinates] - * - */ -template -std::tuple, rmm::device_vector, rmm::device_vector>> -generate_polygon(int32_t num_sides, T radius, vec_2d minXY, vec_2d maxXY) -{ - std::vector polygon_offsets(num_polygons); - std::vector ring_offsets(num_polygons * num_rings_per_polygon); - std::vector> polygon_points(31 * (num_sides + 1)); - - std::iota(polygon_offsets.begin(), polygon_offsets.end(), 0); - std::iota(ring_offsets.begin(), ring_offsets.end(), 0); - std::transform( - ring_offsets.begin(), ring_offsets.end(), ring_offsets.begin(), [num_sides](int32_t i) { - return i * (num_sides + 1); - }); - - for (int32_t i = 0; i < num_polygons; i++) { - auto it = thrust::make_counting_iterator(0); - auto begin = i * num_sides + polygon_points.begin(); - auto center = random_point(minXY, maxXY); - std::transform(it, it + num_sides + 1, begin, [num_sides, radius, center](int32_t j) { - return center + - radius * - vec_2d{ - static_cast(std::cos(2 * PI * (j % num_sides) / static_cast(num_sides))), - static_cast(std::sin(2 * PI * (j % num_sides) / static_cast(num_sides)))}; - }); - } - - // Implicitly convert to device_vector - return std::make_tuple(polygon_offsets, ring_offsets, polygon_points); -} - -/** - * @brief Randomly generate `num_test_points` points within window `minXY` and `maxXY` - * - * @tparam T The floating point type for the coordinates - */ -template -rmm::device_vector> generate_points(int32_t num_test_points, - vec_2d minXY, - vec_2d maxXY) -{ - std::vector> points(num_test_points); - std::generate( - points.begin(), points.end(), [minXY, maxXY]() { return random_point(minXY, maxXY); }); - // Implicitly convert to device_vector - return points; -} +auto constexpr num_polygons = 31ul; +auto constexpr num_rings_per_polygon = 1ul; // only 1 ring for now template void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list) { // TODO: to be replaced by nvbench fixture once it's ready cuspatial::rmm_pool_raii rmm_pool; + rmm::cuda_stream_view stream(rmm::cuda_stream_default); - std::srand(0); // For reproducibility auto const minXY = vec_2d{-radius * 2, -radius * 2}; auto const maxXY = vec_2d{radius * 2, radius * 2}; @@ -128,14 +53,33 @@ void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list) auto const num_polygon_points = num_rings * (num_sides_per_ring + 1); // +1 for the overlapping start and end point of the ring - auto test_points = generate_points(num_test_points, minXY, maxXY); - auto [polygon_offsets, ring_offsets, polygon_points] = - generate_polygon(num_sides_per_ring, radius, minXY, maxXY); - rmm::device_vector result(num_test_points); + auto point_gen_param = test::multipoint_generator_parameter{ + static_cast(num_test_points), 1, minXY, maxXY}; + auto poly_gen_param = + test::multipolygon_generator_parameter{static_cast(num_polygons), + 1, + 0, + static_cast(num_sides_per_ring), + vec_2d{0, 0}, + radius}; + auto test_points = test::generate_multipoint_array(point_gen_param, stream); + auto test_polygons = test::generate_multipolygon_array(poly_gen_param, stream); + + auto [_, points] = test_points.release(); + auto [__, part_offset_array, ring_offset_array, poly_coords] = test_polygons.release(); + + auto points_range = make_multipoint_range( + num_test_points, thrust::make_counting_iterator(0), points.size(), points.begin()); + auto polys_range = make_multipolygon_range(num_polygons, + thrust::make_counting_iterator(0), + part_offset_array.size() - 1, + part_offset_array.begin(), + ring_offset_array.size() - 1, + ring_offset_array.begin(), + poly_coords.size(), + poly_coords.begin()); - auto polygon_offsets_begin = polygon_offsets.begin(); - auto ring_offsets_begin = ring_offsets.begin(); - auto polygon_points_begin = polygon_points.begin(); + rmm::device_vector result(num_test_points); state.add_element_count(num_polygon_points, "NumPolygonPoints"); state.add_global_memory_reads(num_test_points * 2, "TotalMemoryReads"); @@ -145,22 +89,8 @@ void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list) state.add_global_memory_writes(num_test_points, "TotalMemoryWrites"); state.exec(nvbench::exec_tag::sync, - [&test_points, - polygon_offsets_begin, - ring_offsets_begin, - &num_rings, - polygon_points_begin, - &num_polygon_points, - &result](nvbench::launch& launch) { - point_in_polygon(test_points.begin(), - test_points.end(), - polygon_offsets_begin, - polygon_offsets_begin + num_polygons, - ring_offsets_begin, - ring_offsets_begin + num_rings, - polygon_points_begin, - polygon_points_begin + num_polygon_points, - result.begin()); + [points_range, polys_range, &result, stream](nvbench::launch& launch) { + point_in_polygon(points_range, polys_range, result.begin(), stream); }); } diff --git a/cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh b/cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh index 51258e05f..54ad6b499 100644 --- a/cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh +++ b/cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh @@ -99,29 +99,5 @@ __device__ inline bool is_point_in_polygon(vec_2d const& test_point, PolygonR return point_is_within; } -/** - * @brief Compatibility layer with non-OOP style input - */ -template ::difference_type, - class Cart2dItDiffType = typename std::iterator_traits::difference_type> -__device__ inline bool is_point_in_polygon(Cart2d const& test_point, - OffsetType poly_begin, - OffsetType poly_end, - OffsetIterator ring_offsets_first, - OffsetItDiffType const& num_rings, - Cart2dIt poly_points_first, - Cart2dItDiffType const& num_poly_points) -{ - auto polygon = polygon_ref{thrust::next(ring_offsets_first, poly_begin), - thrust::next(ring_offsets_first, poly_end + 1), - poly_points_first, - thrust::next(poly_points_first, num_poly_points)}; - return is_point_in_polygon(test_point, polygon); -} - } // namespace detail } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/point_in_polygon.cuh b/cpp/include/cuspatial/detail/point_in_polygon.cuh index ecf6dd0ba..72c21872a 100644 --- a/cpp/include/cuspatial/detail/point_in_polygon.cuh +++ b/cpp/include/cuspatial/detail/point_in_polygon.cuh @@ -24,240 +24,109 @@ #include +#include + +#include #include +#include #include #include namespace cuspatial { -namespace detail { - -// Get the begin and end offsets of a polygon -template > -__device__ auto poly_begin_end(OffsetIteratorA poly_offsets_first, - OffsetItADiffType const num_polys, - OffsetItADiffType const num_rings, - OffsetType const poly_idx) -{ - auto poly_idx_next = poly_idx + 1; - OffsetType poly_begin = poly_offsets_first[poly_idx]; - OffsetType poly_end = (poly_idx_next < num_polys) ? poly_offsets_first[poly_idx_next] : num_rings; - return std::make_pair(poly_begin, poly_end); -} - -template ::difference_type, - class Cart2dItBDiffType = typename std::iterator_traits::difference_type, - class OffsetItADiffType = typename std::iterator_traits::difference_type, - class OffsetItBDiffType = typename std::iterator_traits::difference_type> -__global__ void point_in_polygon_kernel(Cart2dItA test_points_first, - Cart2dItADiffType const num_test_points, - OffsetIteratorA poly_offsets_first, - OffsetItADiffType const num_polys, - OffsetIteratorB ring_offsets_first, - OffsetItBDiffType const num_rings, - Cart2dItB poly_points_first, - Cart2dItBDiffType const num_poly_points, - OutputIt result) -{ - using Cart2d = iterator_value_type; - using OffsetType = iterator_value_type; - - auto idx = blockIdx.x * blockDim.x + threadIdx.x; - - if (idx >= num_test_points) { return; } - - int32_t hit_mask = 0; - - Cart2d const test_point = test_points_first[idx]; - - // for each polygon - for (auto poly_idx = 0; poly_idx < num_polys; poly_idx++) { - auto [poly_begin, poly_end] = - poly_begin_end(poly_offsets_first, num_polys, num_rings, poly_idx); - - bool const point_is_within = is_point_in_polygon(test_point, - poly_begin, - poly_end, - ring_offsets_first, - num_rings, - poly_points_first, - num_poly_points); - - hit_mask |= point_is_within << poly_idx; - } - result[idx] = hit_mask; -} - -template ::difference_type, - class Cart2dItBDiffType = typename std::iterator_traits::difference_type, - class OffsetItADiffType = typename std::iterator_traits::difference_type, - class OffsetItBDiffType = typename std::iterator_traits::difference_type> -__global__ void pairwise_point_in_polygon_kernel(Cart2dItA test_points_first, - Cart2dItADiffType const num_test_points, - OffsetIteratorA poly_offsets_first, - OffsetItADiffType const num_polys, - OffsetIteratorB ring_offsets_first, - OffsetItBDiffType const num_rings, - Cart2dItB poly_points_first, - Cart2dItBDiffType const num_poly_points, - OutputIt result) -{ - using Cart2d = iterator_value_type; - using OffsetType = iterator_value_type; - for (auto idx = threadIdx.x + blockIdx.x * blockDim.x; idx < num_test_points; - idx += gridDim.x * blockDim.x) { - Cart2d const test_point = test_points_first[idx]; - // for the matching polygon - auto [poly_begin, poly_end] = poly_begin_end(poly_offsets_first, num_polys, num_rings, idx); - - bool const point_is_within = is_point_in_polygon(test_point, - poly_begin, - poly_end, - ring_offsets_first, - num_rings, - poly_points_first, - num_poly_points); - result[idx] = point_is_within; +/** + * @brief Computes point-in-polygon result of a single point to up to 32 polygons + * Result is stored in an `int32_t` integer. + */ +template +struct pip_functor { + PointRange multipoints; + PolygonRange multipolygons; + + int32_t __device__ operator()(std::size_t i) + { + using T = typename PointRange::element_t; + vec_2d point = multipoints[i][0]; + int32_t hit_mask = 0; + for (auto poly_idx = 0; poly_idx < multipolygons.size(); ++poly_idx) + hit_mask |= (is_point_in_polygon(point, multipolygons[poly_idx][0]) << poly_idx); + return hit_mask; } -} +}; -} // namespace detail +template +pip_functor(PointRange, PolygonRange) -> pip_functor; -template -OutputIt point_in_polygon(Cart2dItA test_points_first, - Cart2dItA test_points_last, - OffsetIteratorA polygon_offsets_first, - OffsetIteratorA polygon_offsets_last, - OffsetIteratorB poly_ring_offsets_first, - OffsetIteratorB poly_ring_offsets_last, - Cart2dItB polygon_points_first, - Cart2dItB polygon_points_last, +template +OutputIt point_in_polygon(PointRange points, + PolygonRange polygons, OutputIt output, rmm::cuda_stream_view stream) { - using T = iterator_vec_base_type; + using T = typename PointRange::element_t; - static_assert(is_same_floating_point>(), - "Underlying type of Cart2dItA and Cart2dItB must be the same floating point type"); - static_assert( - is_same, iterator_value_type, iterator_value_type>(), - "Inputs must be cuspatial::vec_2d"); - - static_assert(cuspatial::is_integral, - iterator_value_type>(), - "OffsetIterators must point to integral type."); + static_assert(is_same_floating_point(), + "points and polygons must have the same coordinate type."); static_assert(std::is_same_v, int32_t>, "OutputIt must point to 32 bit integer type."); - auto const num_test_points = std::distance(test_points_first, test_points_last); - - if (num_test_points > 0) { - auto const num_polys = std::distance(polygon_offsets_first, polygon_offsets_last) - 1; - auto const num_rings = std::distance(poly_ring_offsets_first, poly_ring_offsets_last) - 1; - auto const num_poly_points = std::distance(polygon_points_first, polygon_points_last); + CUSPATIAL_EXPECTS(points.num_multipoints() == points.num_points(), + "Point in polygon API only support single point - single polygon tests. " + "Multipoint input is not accepted."); - CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( - num_poly_points, - std::distance(polygon_offsets_first, polygon_offsets_last), - std::distance(poly_ring_offsets_first, poly_ring_offsets_last)); + CUSPATIAL_EXPECTS(polygons.num_multipolygons() == polygons.num_polygons(), + "Point in polygon API only support single point - single polygon tests. " + "MultiPolygon input is not accepted."); - CUSPATIAL_EXPECTS(num_polys <= std::numeric_limits::digits, - "Number of polygons cannot exceed 31"); + CUSPATIAL_EXPECTS(polygons.size() <= std::numeric_limits::digits, + "Number of polygons cannot exceed 31"); - auto [threads_per_block, num_blocks] = grid_1d(num_test_points); + if (points.size() == 0) return output; - detail::point_in_polygon_kernel<<>>( - test_points_first, - num_test_points, - polygon_offsets_first, - num_polys, - poly_ring_offsets_first, - num_rings, - polygon_points_first, - num_poly_points, - output); - CUSPATIAL_CHECK_CUDA(stream.value()); - } + thrust::tabulate( + rmm::exec_policy(stream), output, output + points.size(), pip_functor{points, polygons}); - return output + num_test_points; + return output + points.size(); } -template -OutputIt pairwise_point_in_polygon(Cart2dItA test_points_first, - Cart2dItA test_points_last, - OffsetIteratorA polygon_offsets_first, - OffsetIteratorA polygon_offsets_last, - OffsetIteratorB poly_ring_offsets_first, - OffsetIteratorB poly_ring_offsets_last, - Cart2dItB polygon_points_first, - Cart2dItB polygon_points_last, +template +OutputIt pairwise_point_in_polygon(PointRange points, + PolygonRange polygons, OutputIt output, rmm::cuda_stream_view stream) { - using T = iterator_vec_base_type; + using T = typename PointRange::element_t; - static_assert(is_same_floating_point>(), - "Underlying type of Cart2dItA and Cart2dItB must be the same floating point type"); - static_assert( - is_same, iterator_value_type, iterator_value_type>(), - "Inputs must be cuspatial::vec_2d"); + static_assert(is_same_floating_point(), + "points and polygons must have the same coordinate type."); - static_assert(cuspatial::is_integral, - iterator_value_type>(), - "OffsetIterators must point to integral type."); - - static_assert(std::is_same_v, int32_t>, - "OutputIt must point to 32 bit integer type."); + static_assert(std::is_same_v, uint8_t>, + "OutputIt must be iterator to a uint8_t range."); - auto const num_test_points = std::distance(test_points_first, test_points_last); - auto const num_polys = std::distance(polygon_offsets_first, polygon_offsets_last) - 1; - auto const num_rings = std::distance(poly_ring_offsets_first, poly_ring_offsets_last) - 1; - auto const num_poly_points = std::distance(polygon_points_first, polygon_points_last); + CUSPATIAL_EXPECTS(points.num_multipoints() == points.num_points(), + "Point in polygon API only supports single point - single polygon tests. " + "Multipoint input is not accepted."); - CUSPATIAL_EXPECTS_VALID_POLYGON_SIZES( - num_poly_points, - std::distance(polygon_offsets_first, polygon_offsets_last), - std::distance(poly_ring_offsets_first, poly_ring_offsets_last)); + CUSPATIAL_EXPECTS(polygons.num_multipolygons() == polygons.num_polygons(), + "Point in polygon API only supports single point - single polygon tests. " + "MultiPolygon input is not accepted."); - CUSPATIAL_EXPECTS(num_test_points == num_polys, - "Must pass in an equal number of points and polygons"); + CUSPATIAL_EXPECTS(points.size() == polygons.size(), + "Must pass in an equal number of (multi)points and (multi)polygons"); - auto [threads_per_block, num_blocks] = grid_1d(num_test_points); - detail::pairwise_point_in_polygon_kernel<<>>( - test_points_first, - num_test_points, - polygon_offsets_first, - num_polys, - poly_ring_offsets_first, - num_rings, - polygon_points_first, - num_poly_points, - output); - CUSPATIAL_CHECK_CUDA(stream.value()); + if (points.size() == 0) return output; - return output + num_test_points; + return thrust::transform(rmm::exec_policy(stream), + points.begin(), + points.end(), + polygons.begin(), + output, + [] __device__(auto multipoint, auto multipolygon) { + return is_point_in_polygon(static_cast>(multipoint[0]), + multipolygon[0]); + }); } } // namespace cuspatial diff --git a/cpp/include/cuspatial/detail/range/multipoint_range.cuh b/cpp/include/cuspatial/detail/range/multipoint_range.cuh index 7e9c18cb3..04c3abc34 100644 --- a/cpp/include/cuspatial/detail/range/multipoint_range.cuh +++ b/cpp/include/cuspatial/detail/range/multipoint_range.cuh @@ -38,6 +38,7 @@ struct to_multipoint_functor { GeometryIterator _offset_iter; VecIterator _points_begin; + CUSPATIAL_HOST_DEVICE to_multipoint_functor(GeometryIterator offset_iter, VecIterator points_begin) : _offset_iter(offset_iter), _points_begin(points_begin) { @@ -66,6 +67,9 @@ CUSPATIAL_HOST_DEVICE multipoint_range::multipoin { static_assert(is_vec_2d>, "Coordinate range must be constructed with iterators to vec_2d."); + + static_assert(std::is_integral_v>, + "Offset range must be constructed with iterators to integers."); } template @@ -81,14 +85,14 @@ CUSPATIAL_HOST_DEVICE auto multipoint_range::num_ } template -auto multipoint_range::multipoint_begin() +CUSPATIAL_HOST_DEVICE auto multipoint_range::multipoint_begin() { return cuspatial::detail::make_counting_transform_iterator( 0, detail::to_multipoint_functor(_geometry_begin, _points_begin)); } template -auto multipoint_range::multipoint_end() +CUSPATIAL_HOST_DEVICE auto multipoint_range::multipoint_end() { return multipoint_begin() + size(); } @@ -122,8 +126,7 @@ template CUSPATIAL_HOST_DEVICE auto multipoint_range::operator[]( IndexType idx) { - return multipoint_ref{_points_begin + _geometry_begin[idx], - _points_begin + _geometry_begin[idx + 1]}; + return *(thrust::next(begin(), idx)); } template diff --git a/cpp/include/cuspatial/point_in_polygon.cuh b/cpp/include/cuspatial/point_in_polygon.cuh index 6ae058acd..814d40f2c 100644 --- a/cpp/include/cuspatial/point_in_polygon.cuh +++ b/cpp/include/cuspatial/point_in_polygon.cuh @@ -35,42 +35,27 @@ namespace cuspatial { * represents a hit or miss for each of the input polygons in least-significant-bit order. i.e. * `output[3] & 0b0010` indicates a hit or miss for the 3rd point against the 2nd polygon. * - * - * @tparam Cart2dItA iterator type for point array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam Cart2dItB iterator type for point array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIteratorA iterator type for offset array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIteratorB iterator type for offset array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OutputIt iterator type for output array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI], be device-accessible, mutable and - * iterate on `int32_t` type. - * - * @param test_points_first begin of range of test points - * @param test_points_last end of range of test points - * @param polygon_offsets_first begin of range of indices to the first ring in each polygon - * @param polygon_offsets_last end of range of indices to the first ring in each polygon - * @param ring_offsets_first begin of range of indices to the first point in each ring - * @param ring_offsets_last end of range of indices to the first point in each ring - * @param polygon_points_first begin of range of polygon points - * @param polygon_points_last end of range of polygon points + * Note that the input must be a single geometry column, that is a (multi*)geometry_range + * initialized with counting iterator as the geometry offsets iterator. + * + * @tparam PointRange an instance of template type `multipoint_range`, where + * `GeometryIterator` must be a counting iterator + * @tparam PolygonRange an instance of template type `multipolygon_range`, where + * `GeometryIterator` must be a counting iterator + * @tparam OutputIt iterator type for output array. Must meet the requirements of + * [LegacyRandomAccessIterator][LinkLRAI], be device-accessible, mutable and iterate on `int32_t` + * type. + * + * @param points Range of points, one per computed point-in-polygon pair, + * @param polygons Range of polygons, one per comptued point-in-polygon pair * @param output begin iterator to the output buffer * @param stream The CUDA stream to use for kernel launches. * @return iterator to one past the last element in the output buffer * - * @note Limit 31 polygons per call. Polygons may contain multiple rings. * @note Direction of rings does not matter. - * @note This algorithm supports the ESRI shapefile format, but assumes all polygons are "clean" (as - * defined by the format), and does _not_ verify whether the input adheres to the shapefile format. - * @note The points of the rings can be either explicitly closed (the first and last vertex - * overlaps), or implicitly closed (not overlaps). Either input format is supported. + * @note The points of the rings must be explicitly closed. * @note Overlapping rings negate each other. This behavior is not limited to a single negation, * allowing for "islands" within the same polygon. - * @note `poly_ring_offsets` must contain only the rings that make up the polygons indexed by - * `poly_offsets`. If there are rings in `poly_ring_offsets` that are not part of the polygons in - * `poly_offsets`, results are likely to be incorrect and behavior is undefined. * * ``` * poly w/two rings poly w/four rings @@ -85,78 +70,44 @@ namespace cuspatial { * +-----------+ +------------------------+ * ``` * - * @pre All point iterators must have the same `vec_2d` value type, with the same underlying - * floating-point coordinate type (e.g. `cuspatial::vec_2d`). - * @pre All offset iterators must have the same integral value type. * @pre Output iterator must be mutable and iterate on int32_t type. * - * @throw cuspatial::logic_error if the number of polygons or rings exceeds 31. - * @throw cuspatial::logic_error polygon has less than 1 ring. - * @throw cuspatial::logic_error polygon has less than 4 vertices. - * * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator * "LegacyRandomAccessIterator" */ -template -OutputIt point_in_polygon(Cart2dItA test_points_first, - Cart2dItA test_points_last, - OffsetIteratorA polygon_offsets_first, - OffsetIteratorA polygon_offsets_last, - OffsetIteratorB poly_ring_offsets_first, - OffsetIteratorB poly_ring_offsets_last, - Cart2dItB polygon_points_first, - Cart2dItB polygon_points_last, +template +OutputIt point_in_polygon(PointRange points, + PolygonRange polygons, OutputIt output, rmm::cuda_stream_view stream = rmm::cuda_stream_default); /** - * @brief Given (point, polygon) pairs, tests whether the point of each pair is inside the polygon - * of the pair. - * - * Tests whether each point is inside a corresponding polygon. Points on the edges of the - * polygon are not considered to be inside. - * Polygons are a collection of one or more rings. Rings are a collection of three or more vertices. - * - * Each input point will map to one `int32_t` element in the output. - * - * - * @tparam Cart2dItA iterator type for point array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam Cart2dItB iterator type for point array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIteratorA iterator type for offset array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OffsetIteratorB iterator type for offset array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI] and be device-accessible. - * @tparam OutputIt iterator type for output array. Must meet - * the requirements of [LegacyRandomAccessIterator][LinkLRAI], be device-accessible, mutable and - * iterate on `int32_t` type. - * - * @param test_points_first begin of range of test points - * @param test_points_last end of range of test points - * @param polygon_offsets_first begin of range of indices to the first ring in each polygon - * @param polygon_offsets_last end of range of indices to the first ring in each polygon - * @param ring_offsets_first begin of range of indices to the first point in each ring - * @param ring_offsets_last end of range of indices to the first point in each ring - * @param polygon_points_first begin of range of polygon points - * @param polygon_points_last end of range of polygon points + * @brief Given (point, polygon) pairs, tests whether the point in the pair is in the polygon in the + * pair. + * + * Note that the input must be a single geometry column, that is a (multi*)geometry_range + * initialized with counting iterator as the geometry offsets iterator. + * + * Each input point will map to one `uint8_t` element in the output. + * + * @tparam PointRange an instance of template type `multipoint_range`, where + * `GeometryIterator` must be a counting iterator + * @tparam PolygonRange an instance of template type `multipolygon_range`, where + * `GeometryIterator` must be a counting iterator + * @tparam OutputIt iterator type for output array. Must meet the requirements of + * [LegacyRandomAccessIterator][LinkLRAI], be device-accessible, mutable and iterate on `int32_t` + * type. + * + * @param points Range of points, one per computed point-in-polygon pair, + * @param polygons Range of polygons, one per comptued point-in-polygon pair * @param output begin iterator to the output buffer * @param stream The CUDA stream to use for kernel launches. * @return iterator to one past the last element in the output buffer * * @note Direction of rings does not matter. - * @note This algorithm supports the ESRI shapefile format, but assumes all polygons are "clean" (as - * defined by the format), and does _not_ verify whether the input adheres to the shapefile format. * @note The points of the rings must be explicitly closed. * @note Overlapping rings negate each other. This behavior is not limited to a single negation, * allowing for "islands" within the same polygon. - * @note `poly_ring_offsets` must contain only the rings that make up the polygons indexed by - * `poly_offsets`. If there are rings in `poly_ring_offsets` that are not part of the polygons in - * `poly_offsets`, results are likely to be incorrect and behavior is undefined. * * ``` * poly w/two rings poly w/four rings @@ -171,31 +122,15 @@ OutputIt point_in_polygon(Cart2dItA test_points_first, * +-----------+ +------------------------+ * ``` * - * @pre All point iterators must have the same `vec_2d` value type, with the same underlying - * floating-point coordinate type (e.g. `cuspatial::vec_2d`). - * @pre All offset iterators must have the same integral value type. - * @pre Output iterator must be mutable and iterate on int32_t type. - * - * @throw cuspatial::logic_error polygon has less than 1 ring. - * @throw cuspatial::logic_error polygon has less than 4 vertices. + * @pre Output iterator must be mutable and iterate on uint8_t type. * * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator * "LegacyRandomAccessIterator" */ -template -OutputIt pairwise_point_in_polygon(Cart2dItA test_points_first, - Cart2dItA test_points_last, - OffsetIteratorA polygon_offsets_first, - OffsetIteratorA polygon_offsets_last, - OffsetIteratorB poly_ring_offsets_first, - OffsetIteratorB poly_ring_offsets_last, - Cart2dItB polygon_points_first, - Cart2dItB polygon_points_last, - OutputIt output, +template +OutputIt pairwise_point_in_polygon(PointRange points, + PolygonRange polygons, + OutputIt results, rmm::cuda_stream_view stream = rmm::cuda_stream_default); /** diff --git a/cpp/include/cuspatial/range/multipoint_range.cuh b/cpp/include/cuspatial/range/multipoint_range.cuh index b41edbd7e..92f7bff76 100644 --- a/cpp/include/cuspatial/range/multipoint_range.cuh +++ b/cpp/include/cuspatial/range/multipoint_range.cuh @@ -65,6 +65,7 @@ class multipoint_range { GeometryIterator geometry_end, VecIterator points_begin, VecIterator points_end); + /** * @brief Returns the number of multipoints in the array. */ @@ -83,22 +84,22 @@ class multipoint_range { /** * @brief Returns the iterator to the first multipoint in the multipoint array. */ - auto multipoint_begin(); + CUSPATIAL_HOST_DEVICE auto multipoint_begin(); /** * @brief Returns the iterator past the last multipoint in the multipoint array. */ - auto multipoint_end(); + CUSPATIAL_HOST_DEVICE auto multipoint_end(); /** * @brief Returns the iterator to the start of the multipoint array. */ - auto begin() { return multipoint_begin(); } + CUSPATIAL_HOST_DEVICE auto begin() { return multipoint_begin(); } /** * @brief Returns the iterator past the last multipoint in the multipoint array. */ - auto end() { return multipoint_end(); } + CUSPATIAL_HOST_DEVICE auto end() { return multipoint_end(); } /** * @brief Returns the iterator to the start of the underlying point array. diff --git a/cpp/include/cuspatial/range/multipolygon_range.cuh b/cpp/include/cuspatial/range/multipolygon_range.cuh index 99b2843b2..2af70e066 100644 --- a/cpp/include/cuspatial/range/multipolygon_range.cuh +++ b/cpp/include/cuspatial/range/multipolygon_range.cuh @@ -209,6 +209,68 @@ class multipolygon_range { CUSPATIAL_HOST_DEVICE bool is_valid_segment_id(IndexType1 segment_idx, IndexType2 ring_idx); }; +/** + * @brief Create a multipoylgon_range object of from size and start iterators + * + * @tparam GeometryIteratorDiffType Integer type of the size of the geometry offset array + * @tparam PartIteratorDiffType Integer type of the size of the part offset array + * @tparam RingIteratorDiffType Integer type of the size of the ring offset array + * @tparam VecIteratorDiffType Integer type of the size of the point array + * @tparam GeometryIterator iterator type for offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam PartIterator iterator type for offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam RingIterator iterator type for offset array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * @tparam VecIterator iterator type for the point array. Must meet + * the requirements of [LegacyRandomAccessIterator][LinkLRAI]. + * + * @note Iterators should be device-accessible if the view is intended to be + * used on device. + * + * @param num_multipolygons Number of multipolygons in the array + * @param geometry_begin Iterator to the start of the geometry offset array + * @param num_polygons Number of polygons in the array + * @param part_begin Iterator to the start of the part offset array + * @param num_rings Number of rings in the array + * @param ring_begin Iterator to the start of the ring offset array + * @param num_points Number of underlying points in the multipoint array + * @param point_begin Iterator to the start of the points array + * @return range to multipolygon array + * + * [LinkLRAI]: https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator + * "LegacyRandomAccessIterator" + */ +template +multipolygon_range +make_multipolygon_range(GeometryIteratorDiffType num_multipolygons, + GeometryIterator geometry_begin, + PartIteratorDiffType num_polygons, + PartIterator part_begin, + RingIteratorDiffType num_rings, + RingIterator ring_begin, + VecIteratorDiffType num_points, + VecIterator point_begin) +{ + return multipolygon_range{ + geometry_begin, + thrust::next(geometry_begin, num_multipolygons + 1), + part_begin, + thrust::next(part_begin, num_polygons + 1), + ring_begin, + thrust::next(ring_begin, num_rings + 1), + point_begin, + thrust::next(point_begin, num_points), + }; +} + /** * @brief Create a range object of multipolygon from cuspatial::geometry_column_view. * Specialization for polygons column. diff --git a/cpp/include/cuspatial_test/vector_factories.cuh b/cpp/include/cuspatial_test/vector_factories.cuh index f2572c17e..91d2c04d7 100644 --- a/cpp/include/cuspatial_test/vector_factories.cuh +++ b/cpp/include/cuspatial_test/vector_factories.cuh @@ -93,22 +93,22 @@ class multipolygon_array { multipolygon_array(thrust::device_vector geometry_offsets_array, thrust::device_vector part_offsets_array, thrust::device_vector ring_offsets_array, - thrust::device_vector coordinate_offsets_array) + thrust::device_vector coordinates_array) : _geometry_offsets_array(geometry_offsets_array), _part_offsets_array(part_offsets_array), _ring_offsets_array(ring_offsets_array), - _coordinate_offsets_array(coordinate_offsets_array) + _coordinates_array(coordinates_array) { } multipolygon_array(rmm::device_uvector&& geometry_offsets_array, rmm::device_uvector&& part_offsets_array, rmm::device_uvector&& ring_offsets_array, - rmm::device_uvector&& coordinate_offsets_array) + rmm::device_uvector&& coordinates_array) : _geometry_offsets_array(std::move(geometry_offsets_array)), _part_offsets_array(std::move(part_offsets_array)), _ring_offsets_array(std::move(ring_offsets_array)), - _coordinate_offsets_array(std::move(coordinate_offsets_array)) + _coordinates_array(std::move(coordinates_array)) { } @@ -124,8 +124,8 @@ class multipolygon_array { _part_offsets_array.end(), _ring_offsets_array.begin(), _ring_offsets_array.end(), - _coordinate_offsets_array.begin(), - _coordinate_offsets_array.end()); + _coordinates_array.begin(), + _coordinates_array.end()); } /** @@ -136,11 +136,19 @@ class multipolygon_array { auto geometry_offsets = cuspatial::test::to_host(_geometry_offsets_array); auto part_offsets = cuspatial::test::to_host(_part_offsets_array); auto ring_offsets = cuspatial::test::to_host(_ring_offsets_array); - auto coordinate_offsets = cuspatial::test::to_host(_coordinate_offsets_array); + auto coordinate_offsets = cuspatial::test::to_host(_coordinates_array); return std::tuple{geometry_offsets, part_offsets, ring_offsets, coordinate_offsets}; } + auto release() + { + return std::tuple{std::move(_geometry_offsets_array), + std::move(_part_offsets_array), + std::move(_ring_offsets_array), + std::move(_coordinates_array)}; + } + /** * @brief Output stream operator for `multipolygon_array` for human-readable formatting */ @@ -160,7 +168,7 @@ class multipolygon_array { GeometryArray _geometry_offsets_array; PartArray _part_offsets_array; RingArray _ring_offsets_array; - CoordinateArray _coordinate_offsets_array; + CoordinateArray _coordinates_array; }; template #include #include +#include +#include #include #include @@ -57,7 +59,7 @@ struct point_in_polygon_functor { rmm::mr::device_memory_resource* mr) { auto size = test_points_x.size(); - auto tid = cudf::type_to_id(); + auto tid = pairwise ? cudf::type_to_id() : cudf::type_to_id(); auto type = cudf::data_type{tid}; auto results = cudf::make_fixed_width_column(type, size, cudf::mask_state::UNALLOCATED, stream, mr); @@ -70,31 +72,27 @@ struct point_in_polygon_functor { auto ring_offsets_begin = poly_ring_offsets.begin(); auto polygon_points_begin = cuspatial::make_vec_2d_iterator(poly_points_x.begin(), poly_points_y.begin()); - auto results_begin = results->mutable_view().begin(); - if (pairwise) { - cuspatial::pairwise_point_in_polygon(points_begin, - points_begin + test_points_x.size(), - polygon_offsets_begin, - polygon_offsets_begin + poly_offsets.size(), - ring_offsets_begin, - ring_offsets_begin + poly_ring_offsets.size(), - polygon_points_begin, - polygon_points_begin + poly_points_x.size(), - results_begin, - stream); + auto multipoints_range = + make_multipoint_range(size, thrust::make_counting_iterator(0), size, points_begin); + + auto polygon_size = poly_offsets.size() - 1; + auto multipolygon_range = make_multipolygon_range(polygon_size, + thrust::make_counting_iterator(0), + polygon_size, + polygon_offsets_begin, + poly_ring_offsets.size() - 1, + ring_offsets_begin, + poly_points_x.size(), + polygon_points_begin); + if (pairwise) { + auto results_begin = results->mutable_view().begin(); + cuspatial::pairwise_point_in_polygon( + multipoints_range, multipolygon_range, results_begin, stream); } else { - cuspatial::point_in_polygon(points_begin, - points_begin + test_points_x.size(), - polygon_offsets_begin, - polygon_offsets_begin + poly_offsets.size(), - ring_offsets_begin, - ring_offsets_begin + poly_ring_offsets.size(), - polygon_points_begin, - polygon_points_begin + poly_points_x.size(), - results_begin, - stream); + auto results_begin = results->mutable_view().begin(); + cuspatial::point_in_polygon(multipoints_range, multipolygon_range, results_begin, stream); } return results; diff --git a/cpp/tests/join/quadtree_point_in_polygon_test_large.cu b/cpp/tests/join/quadtree_point_in_polygon_test_large.cu index 8da2a2aed..1be866afc 100644 --- a/cpp/tests/join/quadtree_point_in_polygon_test_large.cu +++ b/cpp/tests/join/quadtree_point_in_polygon_test_large.cu @@ -170,16 +170,21 @@ TYPED_TEST(PIPRefineTestLarge, TestLarge) { // verify rmm::device_uvector hits(points.size(), this->stream()); - auto hits_end = cuspatial::point_in_polygon(points.begin(), - points.end(), - multipolygons.part_offset_begin(), - multipolygons.part_offset_end(), - multipolygons.ring_offset_begin(), - multipolygons.ring_offset_end(), - multipolygons.point_begin(), - multipolygons.point_end(), - hits.begin(), - this->stream()); + + auto points_range = make_multipoint_range( + points.size(), thrust::make_counting_iterator(0), points.size(), points.begin()); + + auto polygons_range = make_multipolygon_range(multipolygons.size(), + thrust::make_counting_iterator(0), + multipolygons.size(), + multipolygons.part_offset_begin(), + multipolygons.num_rings(), + multipolygons.ring_offset_begin(), + multipolygons.num_points(), + multipolygons.point_begin()); + + auto hits_end = + cuspatial::point_in_polygon(points_range, polygons_range, hits.begin(), this->stream()); auto hits_host = cuspatial::test::to_host(hits); diff --git a/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cpp b/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cpp index 477d8ce90..21f939214 100644 --- a/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cpp +++ b/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cpp @@ -50,7 +50,7 @@ TYPED_TEST(PairwisePointInPolygonTest, Empty) auto poly_point_xs = wrapper({}); auto poly_point_ys = wrapper({}); - auto expected = wrapper({}); + auto expected = wrapper({}); auto actual = cuspatial::pairwise_point_in_polygon( test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys); diff --git a/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cu b/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cu index c05c24f63..0dc955204 100644 --- a/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cu +++ b/cpp/tests/point_in_polygon/pairwise_point_in_polygon_test.cu @@ -14,11 +14,14 @@ * limitations under the License. */ +#include +#include +#include + #include #include #include #include -#include #include @@ -26,50 +29,87 @@ #include #include -#include +#include using namespace cuspatial; using namespace cuspatial::test; template -struct PairwisePointInPolygonTest : public ::testing::Test {}; +struct PairwisePointInPolygonTest : public BaseFixture { + void run_test(std::initializer_list> points, + std::initializer_list polygon_offsets, + std::initializer_list ring_offsets, + std::initializer_list> polygon_points, + std::initializer_list expected) + { + auto d_points = make_device_vector>(points); + auto d_polygon_offsets = make_device_vector(polygon_offsets); + auto d_ring_offsets = make_device_vector(ring_offsets); + auto d_polygon_points = make_device_vector>(polygon_points); + + auto mpoints = make_multipoint_range( + d_points.size(), thrust::make_counting_iterator(0), d_points.size(), d_points.begin()); + auto mpolys = make_multipolygon_range(polygon_offsets.size() - 1, + thrust::make_counting_iterator(0), + d_polygon_offsets.size() - 1, + d_polygon_offsets.begin(), + d_ring_offsets.size() - 1, + d_ring_offsets.begin(), + d_polygon_points.size(), + d_polygon_points.begin()); + + auto d_expected = make_device_vector(expected); + + auto got = rmm::device_uvector(points.size(), stream()); + + auto ret = pairwise_point_in_polygon(mpoints, mpolys, got.begin(), stream()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); + EXPECT_EQ(ret, got.end()); + } +}; // float and double are logically the same but would require separate tests due to precision. -using TestTypes = ::testing::Types; -TYPED_TEST_CASE(PairwisePointInPolygonTest, TestTypes); +TYPED_TEST_CASE(PairwisePointInPolygonTest, FloatingPointTypes); TYPED_TEST(PairwisePointInPolygonTest, OnePolygonOneRing) { - using T = TypeParam; - auto point_list = std::vector>{{-2.0, 0.0}, - {2.0, 0.0}, - {0.0, -2.0}, - {0.0, 2.0}, - {-0.5, 0.0}, - {0.5, 0.0}, - {0.0, -0.5}, - {0.0, 0.5}}; + using T = TypeParam; + auto point_list = std::vector>{{-2.0, 0.0}, + {2.0, 0.0}, + {0.0, -2.0}, + {0.0, 2.0}, + {-0.5, 0.0}, + {0.5, 0.0}, + {0.0, -0.5}, + {0.0, 0.5}}; + auto poly_offsets = make_device_vector({0, 1}); auto poly_ring_offsets = make_device_vector({0, 5}); auto poly_point = make_device_vector>( {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}, {-1.0, -1.0}}); - auto got = rmm::device_vector(1); + auto polygon_range = make_multipolygon_range(poly_offsets.size() - 1, + thrust::make_counting_iterator(0), + poly_offsets.size() - 1, + poly_offsets.begin(), + poly_ring_offsets.size() - 1, + poly_ring_offsets.begin(), + poly_point.size(), + poly_point.begin()); + + auto got = rmm::device_vector(1); auto expected = cuspatial::test::make_host_vector({false, false, false, false, true, true, true, true}); for (size_t i = 0; i < point_list.size(); ++i) { - auto point = make_device_vector>({{point_list[i][0], point_list[i][1]}}); - auto ret = pairwise_point_in_polygon(point.begin(), - point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - EXPECT_EQ(got, std::vector({expected[i]})); + auto p = point_list[i]; + auto d_point = make_device_vector>({{p.x, p.y}}); + auto point_range = make_multipoint_range( + d_point.size(), thrust::make_counting_iterator(0), d_point.size(), d_point.begin()); + + auto ret = pairwise_point_in_polygon(point_range, polygon_range, got.begin(), this->stream()); + EXPECT_EQ(got, std::vector({expected[i]})); EXPECT_EQ(ret, got.end()); } } @@ -77,14 +117,14 @@ TYPED_TEST(PairwisePointInPolygonTest, OnePolygonOneRing) TYPED_TEST(PairwisePointInPolygonTest, TwoPolygonsOneRingEach) { using T = TypeParam; - auto point_list = std::vector>{{-2.0, 0.0}, - {2.0, 0.0}, - {0.0, -2.0}, - {0.0, 2.0}, - {-0.5, 0.0}, - {0.5, 0.0}, - {0.0, -0.5}, - {0.0, 0.5}}; + auto point_list = std::vector>{{-2.0, 0.0}, + {2.0, 0.0}, + {0.0, -2.0}, + {0.0, 2.0}, + {-0.5, 0.0}, + {0.5, 0.0}, + {0.0, -0.5}, + {0.0, 0.5}}; auto poly_offsets = make_device_vector({0, 1, 2}); auto poly_ring_offsets = make_device_vector({0, 5, 10}); @@ -99,23 +139,27 @@ TYPED_TEST(PairwisePointInPolygonTest, TwoPolygonsOneRingEach) {-1.0, 0.0}, {0.0, 1.0}}); - auto got = rmm::device_vector(2); - auto expected = std::vector({false, false, false, false, true, true, true, true}); + auto polygon_range = make_multipolygon_range(poly_offsets.size() - 1, + thrust::make_counting_iterator(0), + poly_offsets.size() - 1, + poly_offsets.begin(), + poly_ring_offsets.size() - 1, + poly_ring_offsets.begin(), + poly_point.size(), + poly_point.begin()); + + auto got = rmm::device_vector(2); + auto expected = std::vector({false, false, false, false, true, true, true, true}); for (size_t i = 0; i < point_list.size() / 2; i = i + 2) { auto points = make_device_vector>( - {{point_list[i][0], point_list[i][1]}, {point_list[i + 1][0], point_list[i + 1][1]}}); - auto ret = pairwise_point_in_polygon(points.begin(), - points.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, std::vector({expected[i], expected[i + 1]})); + {{point_list[i].x, point_list[i].y}, {point_list[i + 1].x, point_list[i + 1].y}}); + auto points_range = make_multipoint_range( + points.size(), thrust::make_counting_iterator(0), points.size(), points.begin()); + + auto ret = pairwise_point_in_polygon(points_range, polygon_range, got.begin(), this->stream()); + + EXPECT_EQ(got, std::vector({expected[i], expected[i + 1]})); EXPECT_EQ(ret, got.end()); } } @@ -125,7 +169,9 @@ TYPED_TEST(PairwisePointInPolygonTest, OnePolygonTwoRings) using T = TypeParam; auto point_list = std::vector>{{0.0, 0.0}, {-0.4, 0.0}, {-0.6, 0.0}, {0.0, 0.4}, {0.0, -0.6}}; - auto poly_offsets = make_device_vector({0, 1}); + + auto poly_offsets = make_device_vector({0, 2}); + auto num_polys = poly_offsets.size() - 1; auto poly_ring_offsets = make_device_vector({0, 5, 10}); auto poly_point = make_device_vector>({{-1.0, -1.0}, {1.0, -1.0}, @@ -138,90 +184,64 @@ TYPED_TEST(PairwisePointInPolygonTest, OnePolygonTwoRings) {0.5, -0.5}, {-0.5, -0.5}}); - auto got = rmm::device_vector(1); - auto expected = std::vector{0b0, 0b0, 0b1, 0b0, 0b1}; + auto polygon_range = make_multipolygon_range(num_polys, + thrust::make_counting_iterator(0), + num_polys, + poly_offsets.begin(), + poly_ring_offsets.size() - 1, + poly_ring_offsets.begin(), + poly_point.size(), + poly_point.begin()); + + auto expected = std::vector{0b0, 0b0, 0b1, 0b0, 0b1}; for (size_t i = 0; i < point_list.size(); ++i) { - auto point = make_device_vector>({{point_list[i][0], point_list[i][1]}}); - auto ret = pairwise_point_in_polygon(point.begin(), - point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, std::vector{expected[i]}); + auto got = rmm::device_vector(1); + + auto point = make_device_vector>({{point_list[i][0], point_list[i][1]}}); + auto points_range = make_multipoint_range( + point.size(), thrust::make_counting_iterator(0), point.size(), point.begin()); + + auto ret = pairwise_point_in_polygon(points_range, polygon_range, got.begin(), this->stream()); + + EXPECT_EQ(got, std::vector{expected[i]}); EXPECT_EQ(ret, got.end()); } } TYPED_TEST(PairwisePointInPolygonTest, EdgesOfSquare) { - auto test_point = - make_device_vector>({{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}); - auto poly_offsets = make_device_vector({0, 1, 2, 3, 4}); - auto poly_ring_offsets = make_device_vector({0, 5, 10, 15, 20}); - // 0: rect on min x side // 1: rect on max x side // 2: rect on min y side // 3: rect on max y side - auto poly_point = make_device_vector>( + CUSPATIAL_RUN_TEST( + this->run_test, + {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}, + {0, 1, 2, 3, 4}, + {0, 5, 10, 15, 20}, {{-1.0, -1.0}, {0.0, -1.0}, {0.0, 1.0}, {-1.0, 1.0}, {-1.0, -1.0}, {0.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {0.0, 1.0}, {0.0, -1.0}, {-1.0, -1.0}, {-1.0, 0.0}, {1.0, 0.0}, {1.0, -1.0}, - {-1.0, 1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {-1.0, 0.0}}); - - auto expected = std::vector{0b0, 0b0, 0b0, 0b0}; - auto got = rmm::device_vector(test_point.size()); - - auto ret = pairwise_point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + {-1.0, 1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {-1.0, 0.0}}, + {0b0, 0b0, 0b0, 0b0}); } TYPED_TEST(PairwisePointInPolygonTest, CornersOfSquare) { - auto test_point = - make_device_vector>({{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}); - auto poly_offsets = make_device_vector({0, 1, 2, 3, 4}); - auto poly_ring_offsets = make_device_vector({0, 5, 10, 15, 20}); - // 0: min x min y corner // 1: min x max y corner // 2: max x min y corner // 3: max x max y corner - auto poly_point = make_device_vector>( + + CUSPATIAL_RUN_TEST( + this->run_test, + {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}, + {0, 1, 2, 3, 4}, + {0, 5, 10, 15, 20}, {{-1.0, -1.0}, {-1.0, 0.0}, {0.0, 0.0}, {0.0, -1.0}, {-1.0, -1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {0.0, 1.0}, {-1.0, 0.0}, {-1.0, 0.0}, {0.0, -1.0}, {0.0, 0.0}, {1.0, 0.0}, {1.0, -1.0}, - {0.0, -1.0}, {0.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {0.0, 0.0}}); - - auto expected = std::vector{0b0, 0b0, 0b0, 0b0}; - auto got = rmm::device_vector(test_point.size()); - - auto ret = pairwise_point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + {0.0, -1.0}, {0.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {0.0, 0.0}}, + {0b0, 0b0, 0b0, 0b0}); } struct OffsetIteratorFunctor { @@ -271,6 +291,10 @@ TYPED_TEST(PairwisePointInPolygonTest, 32PolygonSupport) {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}}); + + auto points_range = make_multipoint_range( + test_point.size(), thrust::make_counting_iterator(0), test_point.size(), test_point.begin()); + auto offsets_iter = thrust::make_counting_iterator(0); auto poly_ring_offsets_iter = thrust::make_transform_iterator(offsets_iter, OffsetIteratorFunctor{}); @@ -280,21 +304,22 @@ TYPED_TEST(PairwisePointInPolygonTest, 32PolygonSupport) thrust::make_transform_iterator(offsets_iter, PolyPointIteratorFunctorB{}); auto poly_point_iter = make_vec_2d_iterator(poly_point_xs_iter, poly_point_ys_iter); - auto expected = std::vector({1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, - 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}); - auto got = rmm::device_vector(test_point.size()); - - auto ret = pairwise_point_in_polygon(test_point.begin(), - test_point.end(), - offsets_iter, - offsets_iter + num_polys + 1, - poly_ring_offsets_iter, - poly_ring_offsets_iter + num_polys + 1, - poly_point_iter, - poly_point_iter + num_poly_points, - got.begin()); - - EXPECT_EQ(got, expected); + auto polygons_range = make_multipolygon_range(num_polys, + thrust::make_counting_iterator(0), + num_polys, + offsets_iter, + num_polys, + poly_ring_offsets_iter, + num_poly_points, + poly_point_iter); + + auto expected = std::vector({1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, + 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}); + auto got = rmm::device_vector(test_point.size()); + + auto ret = pairwise_point_in_polygon(points_range, polygons_range, got.begin(), this->stream()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(got, expected); EXPECT_EQ(ret, got.end()); } @@ -304,21 +329,27 @@ TEST_F(PairwisePointInPolygonErrorTest, InsufficientPoints) { using T = double; - auto test_point = make_device_vector>({{0.0, 0.0}, {0.0, 0.0}}); + auto test_point = make_device_vector>({{0.0, 0.0}, {0.0, 0.0}}); + auto points_range = make_multipoint_range( + test_point.size(), thrust::make_counting_iterator(0), test_point.size(), test_point.begin()); + auto poly_offsets = make_device_vector({0, 1}); + auto num_polys = poly_offsets.size() - 1; auto poly_ring_offsets = make_device_vector({0, 3}); auto poly_point = make_device_vector>({{0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}}); - auto got = rmm::device_vector(test_point.size()); - - EXPECT_THROW(pairwise_point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()), + + auto polygons_range = make_multipolygon_range(num_polys, + thrust::make_counting_iterator(0), + num_polys, + poly_offsets.begin(), + num_polys, + poly_ring_offsets.begin(), + poly_point.size(), + poly_point.begin()); + + auto got = rmm::device_vector(test_point.size()); + + EXPECT_THROW(pairwise_point_in_polygon(points_range, polygons_range, got.begin(), this->stream()), cuspatial::logic_error); } @@ -326,21 +357,27 @@ TEST_F(PairwisePointInPolygonErrorTest, InsufficientPolyOffsets) { using T = double; - auto test_point = make_device_vector>({{0.0, 0.0}, {0.0, 0.0}}); + auto test_point = make_device_vector>({{0.0, 0.0}, {0.0, 0.0}}); + auto points_range = make_multipoint_range( + test_point.size(), thrust::make_counting_iterator(0), test_point.size(), test_point.begin()); + auto poly_offsets = make_device_vector({0}); + auto num_polys = poly_offsets.size() - 1; auto poly_ring_offsets = make_device_vector({0, 4}); auto poly_point = make_device_vector>({{0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}, {0.0, 1.0}}); - auto got = rmm::device_vector(test_point.size()); - - EXPECT_THROW(pairwise_point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()), + + auto polygons_range = make_multipolygon_range(num_polys, + thrust::make_counting_iterator(0), + num_polys, + poly_offsets.begin(), + num_polys, + poly_ring_offsets.begin(), + poly_point.size(), + poly_point.begin()); + + auto got = rmm::device_vector(test_point.size()); + + EXPECT_THROW(pairwise_point_in_polygon(points_range, polygons_range, got.begin(), this->stream()), cuspatial::logic_error); } diff --git a/cpp/tests/point_in_polygon/point_in_polygon_test.cu b/cpp/tests/point_in_polygon/point_in_polygon_test.cu index 00e4229ed..c4705cad9 100644 --- a/cpp/tests/point_in_polygon/point_in_polygon_test.cu +++ b/cpp/tests/point_in_polygon/point_in_polygon_test.cu @@ -14,6 +14,7 @@ * limitations under the License. */ +#include #include #include @@ -32,55 +33,61 @@ #include using namespace cuspatial; +using namespace cuspatial::test; template -struct PointInPolygonTest : public ::testing::Test { +struct PointInPolygonTest : public BaseFixture { public: - rmm::device_vector> make_device_points(std::initializer_list> pts) + void run_test(std::initializer_list> points, + std::initializer_list polygon_offsets, + std::initializer_list ring_offsets, + std::initializer_list> polygon_points, + std::initializer_list expected) { - return rmm::device_vector>(pts.begin(), pts.end()); - } - - rmm::device_vector make_device_offsets(std::initializer_list pts) - { - return rmm::device_vector(pts.begin(), pts.end()); + auto d_points = make_device_vector>(points); + auto d_polygon_offsets = make_device_vector(polygon_offsets); + auto d_ring_offsets = make_device_vector(ring_offsets); + auto d_polygon_points = make_device_vector>(polygon_points); + + auto mpoints = make_multipoint_range( + d_points.size(), thrust::make_counting_iterator(0), d_points.size(), d_points.begin()); + auto mpolys = make_multipolygon_range(polygon_offsets.size() - 1, + thrust::make_counting_iterator(0), + d_polygon_offsets.size() - 1, + d_polygon_offsets.begin(), + d_ring_offsets.size() - 1, + d_ring_offsets.begin(), + d_polygon_points.size(), + d_polygon_points.begin()); + + auto d_expected = make_device_vector(expected); + + auto got = rmm::device_uvector(points.size(), stream()); + + auto ret = point_in_polygon(mpoints, mpolys, got.begin(), stream()); + + CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(d_expected, got); + EXPECT_EQ(ret, got.end()); } }; -// float and double are logically the same but would require separate tests due to precision. -using TestTypes = ::testing::Types; -TYPED_TEST_CASE(PointInPolygonTest, TestTypes); +TYPED_TEST_CASE(PointInPolygonTest, FloatingPointTypes); TYPED_TEST(PointInPolygonTest, OnePolygonOneRing) { - auto test_point = this->make_device_points({{-2.0, 0.0}, - {2.0, 0.0}, - {0.0, -2.0}, - {0.0, 2.0}, - {-0.5, 0.0}, - {0.5, 0.0}, - {0.0, -0.5}, - {0.0, 0.5}}); - auto poly_offsets = this->make_device_offsets({0, 1}); - auto poly_ring_offsets = this->make_device_offsets({0, 5}); - auto poly_point = - this->make_device_points({{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}, {-1.0, -1.0}}); - - auto got = rmm::device_vector(test_point.size()); - auto expected = std::vector{false, false, false, false, true, true, true, true}; - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + CUSPATIAL_RUN_TEST(this->run_test, + {{-2.0, 0.0}, + {2.0, 0.0}, + {0.0, -2.0}, + {0.0, 2.0}, + {-0.5, 0.0}, + {0.5, 0.0}, + {0.0, -0.5}, + {0.0, 0.5}}, + {0, 1}, + {0, 5}, + {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}, {-1.0, -1.0}}, + {false, false, false, false, true, true, true, true}); } // cuspatial expects closed rings, however algorithms may work OK with unclosed rings @@ -90,172 +97,106 @@ TYPED_TEST(PointInPolygonTest, OnePolygonOneRing) // uses a polygon ring with 4 vertices so it doesn't fail polygon validation. TYPED_TEST(PointInPolygonTest, OnePolygonOneRingUnclosed) { - auto test_point = this->make_device_points({{-2.0, 0.0}, - {2.0, 0.0}, - {0.0, -2.0}, - {0.0, 2.0}, - {-0.5, 0.0}, - {0.5, 0.0}, - {0.0, -0.5}, - {0.0, 0.5}}); - auto poly_offsets = this->make_device_offsets({0, 1}); - auto poly_ring_offsets = this->make_device_offsets({0, 4}); - auto poly_point = this->make_device_points({{-1.0, -1.0}, {1.0, -1.0}, {1.0, 0.0}, {1.0, 1.0}}); - - auto got = rmm::device_vector(test_point.size()); - auto expected = std::vector{false, false, false, false, false, true, true, false}; - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + CUSPATIAL_RUN_TEST(this->run_test, + {{-2.0, 0.0}, + {2.0, 0.0}, + {0.0, -2.0}, + {0.0, 2.0}, + {-0.5, 0.0}, + {0.5, 0.0}, + {0.0, -0.5}, + {0.0, 0.5}}, + {0, 1}, + {0, 4}, + {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 0.0}, {1.0, 1.0}}, + {false, false, false, false, false, true, true, false}); } TYPED_TEST(PointInPolygonTest, TwoPolygonsOneRingEach) { - auto test_point = this->make_device_points({{-2.0, 0.0}, - {2.0, 0.0}, - {0.0, -2.0}, - {0.0, 2.0}, - {-0.5, 0.0}, - {0.5, 0.0}, - {0.0, -0.5}, - {0.0, 0.5}}); - - auto poly_offsets = this->make_device_offsets({0, 1, 2}); - auto poly_ring_offsets = this->make_device_offsets({0, 5, 10}); - auto poly_point = this->make_device_points({{-1.0, -1.0}, - {-1.0, 1.0}, - {1.0, 1.0}, - {1.0, -1.0}, - {-1.0, -1.0}, - {0.0, 1.0}, - {1.0, 0.0}, - {0.0, -1.0}, - {-1.0, 0.0}, - {0.0, 1.0}}); - - auto got = rmm::device_vector(test_point.size()); - auto expected = std::vector({0b00, 0b00, 0b00, 0b00, 0b11, 0b11, 0b11, 0b11}); - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + CUSPATIAL_RUN_TEST(this->run_test, + + {{-2.0, 0.0}, + {2.0, 0.0}, + {0.0, -2.0}, + {0.0, 2.0}, + {-0.5, 0.0}, + {0.5, 0.0}, + {0.0, -0.5}, + {0.0, 0.5}}, + + {0, 1, 2}, + {0, 5, 10}, + {{-1.0, -1.0}, + {-1.0, 1.0}, + {1.0, 1.0}, + {1.0, -1.0}, + {-1.0, -1.0}, + {0.0, 1.0}, + {1.0, 0.0}, + {0.0, -1.0}, + {-1.0, 0.0}, + {0.0, 1.0}}, + + {0b00, 0b00, 0b00, 0b00, 0b11, 0b11, 0b11, 0b11}); } TYPED_TEST(PointInPolygonTest, OnePolygonTwoRings) { - auto test_point = - this->make_device_points({{0.0, 0.0}, {-0.4, 0.0}, {-0.6, 0.0}, {0.0, 0.4}, {0.0, -0.6}}); - auto poly_offsets = this->make_device_offsets({0, 2}); - auto poly_ring_offsets = this->make_device_offsets({0, 5, 10}); - auto poly_point = this->make_device_points({{-1.0, -1.0}, - {1.0, -1.0}, - {1.0, 1.0}, - {-1.0, 1.0}, - {-1.0, -1.0}, - {-0.5, -0.5}, - {-0.5, 0.5}, - {0.5, 0.5}, - {0.5, -0.5}, - {-0.5, -0.5}}); - - auto got = rmm::device_vector(test_point.size()); - auto expected = std::vector{0b0, 0b0, 0b1, 0b0, 0b1}; - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + CUSPATIAL_RUN_TEST(this->run_test, + {{0.0, 0.0}, {-0.4, 0.0}, {-0.6, 0.0}, {0.0, 0.4}, {0.0, -0.6}}, + {0, 2}, + {0, 5, 10}, + {{-1.0, -1.0}, + {1.0, -1.0}, + {1.0, 1.0}, + {-1.0, 1.0}, + {-1.0, -1.0}, + {-0.5, -0.5}, + {-0.5, 0.5}, + {0.5, 0.5}, + {0.5, -0.5}, + {-0.5, -0.5}}, + + {0b0, 0b0, 0b1, 0b0, 0b1}); } TYPED_TEST(PointInPolygonTest, EdgesOfSquare) { - auto test_point = this->make_device_points({{0.0, 0.0}}); - auto poly_offsets = this->make_device_offsets({0, 1, 2, 3, 4}); - auto poly_ring_offsets = this->make_device_offsets({0, 5, 10, 15, 20}); - - // 0: rect on min x side - // 1: rect on max x side - // 2: rect on min y side - // 3: rect on max y side - auto poly_point = this->make_device_points( + CUSPATIAL_RUN_TEST( + this->run_test, + {{0.0, 0.0}}, + {0, 1, 2, 3, 4}, + {0, 5, 10, 15, 20}, + + // 0: rect on min x side + // 1: rect on max x side + // 2: rect on min y side + // 3: rect on max y side {{-1.0, -1.0}, {0.0, -1.0}, {0.0, 1.0}, {-1.0, 1.0}, {-1.0, -1.0}, {0.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {0.0, 1.0}, {0.0, -1.0}, {-1.0, -1.0}, {-1.0, 0.0}, {1.0, 0.0}, {1.0, -1.0}, - {-1.0, 1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {-1.0, 0.0}}); - - auto expected = std::vector{0b0000}; - auto got = rmm::device_vector(test_point.size()); - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); + {-1.0, 1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {-1.0, 0.0}}, - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + {0b0000}); } TYPED_TEST(PointInPolygonTest, CornersOfSquare) { - auto test_point = this->make_device_points({{0.0, 0.0}}); - auto poly_offsets = this->make_device_offsets({0, 1, 2, 3, 4}); - auto poly_ring_offsets = this->make_device_offsets({0, 5, 10, 15, 20}); - - // 0: min x min y corner - // 1: min x max y corner - // 2: max x min y corner - // 3: max x max y corner - auto poly_point = this->make_device_points( + CUSPATIAL_RUN_TEST( + this->run_test, + {{0.0, 0.0}}, + {0, 1, 2, 3, 4}, + {0, 5, 10, 15, 20}, + + // 0: min x min y corner + // 1: min x max y corner + // 2: max x min y corner + // 3: max x max y corner {{-1.0, -1.0}, {-1.0, 0.0}, {0.0, 0.0}, {0.0, -1.0}, {-1.0, -1.0}, {-1.0, 0.0}, {-1.0, 1.0}, {0.0, 1.0}, {-1.0, 0.0}, {-1.0, 0.0}, {0.0, -1.0}, {0.0, 0.0}, {1.0, 0.0}, {1.0, -1.0}, - {0.0, -1.0}, {0.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {0.0, 0.0}}); - - auto expected = std::vector{0b0000}; - auto got = rmm::device_vector(test_point.size()); + {0.0, -1.0}, {0.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}, {1.0, 0.0}, {0.0, 0.0}}, - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(got, expected); - EXPECT_EQ(ret, got.end()); + {0b0000}); } struct OffsetIteratorFunctor { @@ -299,7 +240,7 @@ TYPED_TEST(PointInPolygonTest, 31PolygonSupport) auto constexpr num_polys = 31; auto constexpr num_poly_points = num_polys * 5; - auto test_point = this->make_device_points({{0.0, 0.0}, {2.0, 0.0}}); + auto test_point = make_device_vector>({{0.0, 0.0}, {2.0, 0.0}}); auto offsets_iter = thrust::make_counting_iterator(0); auto poly_ring_offsets_iter = thrust::make_transform_iterator(offsets_iter, OffsetIteratorFunctor{}); @@ -309,103 +250,60 @@ TYPED_TEST(PointInPolygonTest, 31PolygonSupport) thrust::make_transform_iterator(offsets_iter, PolyPointIteratorFunctorB{}); auto poly_point_iter = make_vec_2d_iterator(poly_point_xs_iter, poly_point_ys_iter); + auto points_range = make_multipoint_range( + test_point.size(), thrust::make_counting_iterator(0), test_point.size(), test_point.begin()); + + auto polygons_range = make_multipolygon_range(num_polys, + thrust::make_counting_iterator(0), + num_polys, + offsets_iter, + num_polys, + poly_ring_offsets_iter, + num_poly_points, + poly_point_iter); + auto expected = std::vector({0b1111111111111111111111111111111, 0b0000000000000000000000000000000}); auto got = rmm::device_vector(test_point.size()); - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - offsets_iter, - offsets_iter + num_polys + 1, - poly_ring_offsets_iter, - poly_ring_offsets_iter + num_polys + 1, - poly_point_iter, - poly_point_iter + num_poly_points, - got.begin()); + auto ret = point_in_polygon(points_range, polygons_range, got.begin()); EXPECT_EQ(got, expected); EXPECT_EQ(ret, got.end()); } -struct PointInPolygonErrorTest : public PointInPolygonTest {}; - TYPED_TEST(PointInPolygonTest, SelfClosingLoopLeftEdgeMissing) { - using T = TypeParam; - auto test_point = this->make_device_points({{-2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}}); - auto poly_offsets = this->make_device_offsets({0, 1}); - auto poly_ring_offsets = this->make_device_offsets({0, 4}); - // "left" edge missing - auto poly_point = this->make_device_points({{-1, 1}, {1, 1}, {1, -1}, {-1, -1}}); - auto expected = std::vector{0b0, 0b1, 0b0}; - auto got = rmm::device_vector(test_point.size()); - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(expected, got); - EXPECT_EQ(got.end(), ret); + CUSPATIAL_RUN_TEST(this->run_test, + + {{-2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}}, + {0, 1}, + {0, 4}, + // "left" edge missing + {{-1, 1}, {1, 1}, {1, -1}, {-1, -1}}, + {0b0, 0b1, 0b0}); } TYPED_TEST(PointInPolygonTest, SelfClosingLoopRightEdgeMissing) { - using T = TypeParam; - auto test_point = this->make_device_points({{-2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}}); - auto poly_offsets = this->make_device_offsets({0, 1}); - auto poly_ring_offsets = this->make_device_offsets({0, 4}); - // "right" edge missing - auto poly_point = this->make_device_points({{1, -1}, {-1, -1}, {-1, 1}, {1, 1}}); - auto expected = std::vector{0b0, 0b1, 0b0}; - auto got = rmm::device_vector(test_point.size()); - - auto ret = point_in_polygon(test_point.begin(), - test_point.end(), - poly_offsets.begin(), - poly_offsets.end(), - poly_ring_offsets.begin(), - poly_ring_offsets.end(), - poly_point.begin(), - poly_point.end(), - got.begin()); - - EXPECT_EQ(expected, got); - EXPECT_EQ(got.end(), ret); + using T = TypeParam; + CUSPATIAL_RUN_TEST(this->run_test, + {{-2.0, 0.0}, {0.0, 0.0}, {2.0, 0.0}}, + {0, 1}, + {0, 4}, + // "right" edge missing + {{1, -1}, {-1, -1}, {-1, 1}, {1, 1}}, + {0b0, 0b1, 0b0}); } TYPED_TEST(PointInPolygonTest, ContainsButCollinearWithBoundary) { using T = TypeParam; - - auto point = cuspatial::test::make_multipoint_array({{{0.5, 0.5}}}); - auto polygon = cuspatial::test::make_multipolygon_array( - {0, 1}, + CUSPATIAL_RUN_TEST( + this->run_test, + {{0.5, 0.5}}, {0, 1}, {0, 9}, - {{0, 0}, {0, 1}, {1, 1}, {1, 0.5}, {1.5, 0.5}, {1.5, 1}, {2, 1}, {2, 0}, {0, 0}}); - - auto point_range = point.range(); - auto polygon_range = polygon.range(); - - auto res = rmm::device_uvector(1, rmm::cuda_stream_default); - - cuspatial::point_in_polygon(point_range.point_begin(), - point_range.point_end(), - polygon_range.part_offset_begin(), - polygon_range.part_offset_end(), - polygon_range.ring_offset_begin(), - polygon_range.ring_offset_end(), - polygon_range.point_begin(), - polygon_range.point_end(), - res.begin()); - - auto expect = cuspatial::test::make_device_vector({0b1}); - - CUSPATIAL_EXPECT_VECTORS_EQUIVALENT(res, expect); + {{0, 0}, {0, 1}, {1, 1}, {1, 0.5}, {1.5, 0.5}, {1.5, 1}, {2, 1}, {2, 0}, {0, 0}}, + {0b1}); }