Skip to content

Commit

Permalink
Merge branch 'branch-23.08' into bug/improve-polygon-linestring-contains
Browse files Browse the repository at this point in the history
  • Loading branch information
@thomcom
thomcom committed Jun 23, 2023
2 parents aa62649 + d5499d2 commit 8098da6
Show file tree
Hide file tree
Showing 34 changed files with 677 additions and 927 deletions.
136 changes: 33 additions & 103 deletions cpp/benchmarks/point_in_polygon/point_in_polygon.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@
#include <benchmarks/fixture/rmm_pool_raii.hpp>
#include <nvbench/nvbench.cuh>

#include <cuspatial_test/geometry_generator.cuh>

#include <cuspatial/geometry/vec_2d.hpp>
#include <cuspatial/point_in_polygon.cuh>

Expand All @@ -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 <typename T>
vec_2d<T> random_point(vec_2d<T> minXY, vec_2d<T> maxXY)
{
auto x = minXY.x + (maxXY.x - minXY.x) * rand() / static_cast<T>(RAND_MAX);
auto y = minXY.y + (maxXY.y - minXY.y) * rand() / static_cast<T>(RAND_MAX);
return vec_2d<T>{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 <typename T>
std::tuple<rmm::device_vector<int32_t>, rmm::device_vector<int32_t>, rmm::device_vector<vec_2d<T>>>
generate_polygon(int32_t num_sides, T radius, vec_2d<T> minXY, vec_2d<T> maxXY)
{
std::vector<int32_t> polygon_offsets(num_polygons);
std::vector<int32_t> ring_offsets(num_polygons * num_rings_per_polygon);
std::vector<vec_2d<T>> 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<T>{
static_cast<T>(std::cos(2 * PI * (j % num_sides) / static_cast<T>(num_sides))),
static_cast<T>(std::sin(2 * PI * (j % num_sides) / static_cast<T>(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 <typename T>
rmm::device_vector<vec_2d<T>> generate_points(int32_t num_test_points,
vec_2d<T> minXY,
vec_2d<T> maxXY)
{
std::vector<vec_2d<T>> 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 <typename T>
void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list<T>)
{
// 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<T>{-radius * 2, -radius * 2};
auto const maxXY = vec_2d<T>{radius * 2, radius * 2};

Expand All @@ -128,14 +53,33 @@ void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list<T>)
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<T>(num_test_points, minXY, maxXY);
auto [polygon_offsets, ring_offsets, polygon_points] =
generate_polygon<T>(num_sides_per_ring, radius, minXY, maxXY);
rmm::device_vector<int32_t> result(num_test_points);
auto point_gen_param = test::multipoint_generator_parameter<T>{
static_cast<std::size_t>(num_test_points), 1, minXY, maxXY};
auto poly_gen_param =
test::multipolygon_generator_parameter<T>{static_cast<std::size_t>(num_polygons),
1,
0,
static_cast<std::size_t>(num_sides_per_ring),
vec_2d<T>{0, 0},
radius};
auto test_points = test::generate_multipoint_array<T>(point_gen_param, stream);
auto test_polygons = test::generate_multipolygon_array<T>(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<int32_t> result(num_test_points);

state.add_element_count(num_polygon_points, "NumPolygonPoints");
state.add_global_memory_reads<T>(num_test_points * 2, "TotalMemoryReads");
Expand All @@ -145,22 +89,8 @@ void point_in_polygon_benchmark(nvbench::state& state, nvbench::type_list<T>)
state.add_global_memory_writes<int32_t>(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);
});
}

Expand Down
24 changes: 0 additions & 24 deletions cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,29 +99,5 @@ __device__ inline bool is_point_in_polygon(vec_2d<T> const& test_point, PolygonR
return point_is_within;
}

/**
* @brief Compatibility layer with non-OOP style input
*/
template <class Cart2d,
class OffsetType,
class OffsetIterator,
class Cart2dIt,
class OffsetItDiffType = typename std::iterator_traits<OffsetIterator>::difference_type,
class Cart2dItDiffType = typename std::iterator_traits<Cart2dIt>::difference_type>
__device__ inline bool is_point_in_polygon(Cart2d const& test_point,
OffsetType poly_begin,
OffsetType poly_end,
OffsetIterator ring_offsets_first,
OffsetItDiffType const& num_rings,
Cart2dIt poly_points_first,
Cart2dItDiffType const& num_poly_points)
{
auto polygon = polygon_ref{thrust::next(ring_offsets_first, poly_begin),
thrust::next(ring_offsets_first, poly_end + 1),
poly_points_first,
thrust::next(poly_points_first, num_poly_points)};
return is_point_in_polygon(test_point, polygon);
}

} // namespace detail
} // namespace cuspatial
Loading

0 comments on commit 8098da6

Please sign in to comment.