Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Header-only Refactor of point_in_polygon #587

Merged
merged 31 commits into from
Jul 27, 2022
Merged

Header-only Refactor of point_in_polygon #587

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
6cb6180
Initial conversion
isVoid Jul 6, 2022
7133e34
add first header only test
isVoid Jul 6, 2022
4efbe67
add second test
isVoid Jul 6, 2022
a6209e1
Add more tests
isVoid Jul 8, 2022
c8afc71
Separate the two offset iterator types
isVoid Jul 8, 2022
f7ab8c2
style
isVoid Jul 8, 2022
aef327d
Remove some type tests
isVoid Jul 8, 2022
acef20e
Replace iterator_trait with helpers
isVoid Jul 8, 2022
b3d5b7a
Add folded is_integral
isVoid Jul 8, 2022
e1d0e30
move iterator type helpers to traits.hpp
isVoid Jul 8, 2022
0432721
Remove overlapping tests in cudf API
isVoid Jul 8, 2022
27fb73c
documentation
isVoid Jul 8, 2022
794f6ad
Adopt explict type for references to iterator values
isVoid Jul 13, 2022
d4fd3dc
Factor out the *single* point in polygon logic
isVoid Jul 13, 2022
56fdbd8
Refactor single point test kernel interface
isVoid Jul 14, 2022
fe9c220
remove experiemental prefix in doxygen_groups.h
isVoid Jul 14, 2022
f192862
Add todo
isVoid Jul 14, 2022
ae44ef3
fix typo
isVoid Jul 14, 2022
0c5298c
Merge branch 'branch-22.08' of https://github.com/rapidsai/cuspatial …
isVoid Jul 20, 2022
61366d0
Introduce point in polygon benchmark
isVoid Jul 20, 2022
399e80d
Implement Eric's style of "wrap around" technique
isVoid Jul 21, 2022
041d722
Further remove the division in the pip test
isVoid Jul 21, 2022
1db2574
Reintroduce the type tests for cudf API pip test.
isVoid Jul 21, 2022
c5eaa2d
style
isVoid Jul 21, 2022
f015f28
Update docs
isVoid Jul 21, 2022
047551d
More docstrings
isVoid Jul 22, 2022
e3a2134
Add comment on potential impact on divergence
isVoid Jul 27, 2022
faaf9a3
docs update
isVoid Jul 27, 2022
8a7b34b
Merge branch 'branch-22.08' of https://github.com/rapidsai/cuspatial …
isVoid Jul 27, 2022
841b69d
Update cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh
isVoid Jul 27, 2022
901a4d8
style
isVoid Jul 27, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 21 additions & 20 deletions cpp/include/cuspatial/experimental/detail/point_in_polygon.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,29 +62,30 @@ __device__ inline bool is_point_in_polygon(Cart2d const& test_point,
bool point_is_within = false;
// for each ring
for (auto ring_idx = poly_begin; ring_idx < poly_end; ring_idx++) {
auto ring_idx_next = ring_idx + 1;
auto ring_begin = ring_offsets_first[ring_idx];
auto ring_end =
int32_t ring_idx_next = ring_idx + 1;
int32_t ring_begin = ring_offsets_first[ring_idx];
int32_t ring_end =
(ring_idx_next < num_rings) ? ring_offsets_first[ring_idx_next] : num_poly_points;
auto ring_len = ring_end - ring_begin;

Cart2d b = poly_points_first[ring_end - 1];
bool y0_flag = b.y > test_point.y;
bool y1_flag;
// for each line segment, including the segment between the last and first vertex
for (auto point_idx = 0; point_idx < ring_len; point_idx++) {
Cart2d const a = poly_points_first[ring_begin + ((point_idx + 0) % ring_len)];
Cart2d const b = poly_points_first[ring_begin + ((point_idx + 1) % ring_len)];

bool y_between_ay_by =
a.y <= test_point.y && test_point.y < b.y; // is y in range [ay, by) when ay < by?
bool y_between_by_ay =
b.y <= test_point.y && test_point.y < a.y; // is y in range [by, ay) when by < ay?
bool y_in_bounds = y_between_ay_by || y_between_by_ay; // is y in range [by, ay]?
T run = b.x - a.x;
T rise = b.y - a.y;
T rise_to_point = test_point.y - a.y;

if (y_in_bounds && test_point.x < (run / rise) * rise_to_point + a.x) {
point_is_within = not point_is_within;
for (auto point_idx = ring_begin; point_idx < ring_end; point_idx++) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to know how much of the speedup came just from eliminating the division vs. all the other changes. Because there is potential for more thread divergence in the new code.

Copy link
Contributor Author

@isVoid isVoid Jul 27, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will include the comparison data of division in the PR description. I will also include a comment stating the assumption we have here is that the cycles it costs on division is greater than that of warp syncs. It's apparently dependent on many factors such as whether user compiles the code with fast-math, or how fast the hardware synchronizes divergence.

That said, I don't think it's appropriate to mention benchmark results in the comments as it's hardware dependent and will get stale very quickly.

Cart2d const a = poly_points_first[point_idx];
y1_flag = a.y > test_point.y;
if (y1_flag != y0_flag) {
T run = b.x - a.x;
T rise = b.y - a.y;
T rise_to_point = test_point.y - a.y;

if (rise > 0 && (test_point.x - a.x) * rise < run * rise_to_point)
point_is_within = not point_is_within;
else if (rise < 0 && (test_point.x - a.x) * rise > run * rise_to_point)
point_is_within = not point_is_within;
}
b = a;
y0_flag = y1_flag;
}
}

Expand Down Expand Up @@ -115,7 +116,7 @@ __global__ void point_in_polygon_kernel(Cart2dItA test_points_first,

auto idx = blockIdx.x * blockDim.x + threadIdx.x;

if (idx > num_test_points) { return; }
if (idx >= num_test_points) { return; }

int32_t hit_mask = 0;

Expand Down
134 changes: 134 additions & 0 deletions cpp/tests/spatial/point_in_polygon_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,3 +58,137 @@ TYPED_TEST(PointInPolygonTest, Empty)

expect_columns_equal(expected, actual->view(), verbosity);
}

template <typename T>
struct PointInPolygonUnsupportedTypesTest : public BaseFixture {
};

using UnsupportedTestTypes = RemoveIf<ContainedIn<TestTypes>, NumericTypes>;
TYPED_TEST_CASE(PointInPolygonUnsupportedTypesTest, UnsupportedTestTypes);

TYPED_TEST(PointInPolygonUnsupportedTypesTest, UnsupportedPointType)
{
using T = TypeParam;

auto test_point_xs = wrapper<T>({0.0, 0.0});
auto test_point_ys = wrapper<T>({0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T>({0.0, 1.0, 0.0, -1.0});
auto poly_point_ys = wrapper<T>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

template <typename T>
struct PointInPolygonUnsupportedChronoTypesTest : public BaseFixture {
};

TYPED_TEST_CASE(PointInPolygonUnsupportedChronoTypesTest, ChronoTypes);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't remember why we have unsupported types tests specifically for Chrono types but not for other non-numerical types (e.g. nested types). Is it necessary to test these at all? Every type we test increases compile time, so I wonder if we can just test one unsupported type as a proxy for all the others?

Don't necessarily need to decide this in the present PR!


TYPED_TEST(PointInPolygonUnsupportedChronoTypesTest, UnsupportedPointChronoType)
{
using T = TypeParam;
using R = typename T::rep;

auto test_point_xs = wrapper<T, R>({R{0}, R{0}});
auto test_point_ys = wrapper<T, R>({R{0}});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T, R>({R{0}, R{1}, R{0}, R{-1}});
auto poly_point_ys = wrapper<T, R>({R{1}, R{0}, R{-1}, R{0}});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

struct PointInPolygonErrorTest : public BaseFixture {
};

TEST_F(PointInPolygonErrorTest, MismatchTestPointXYLength)
{
using T = double;

auto test_point_xs = wrapper<T>({0.0, 0.0});
auto test_point_ys = wrapper<T>({0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T>({0.0, 1.0, 0.0, -1.0});
auto poly_point_ys = wrapper<T>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

TEST_F(PointInPolygonErrorTest, MismatchTestPointType)
{
using T = double;

auto test_point_xs = wrapper<T>({0.0, 0.0});
auto test_point_ys = wrapper<float>({0.0, 0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T>({0.0, 1.0, 0.0});
auto poly_point_ys = wrapper<T>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

TEST_F(PointInPolygonErrorTest, MismatchPolyPointXYLength)
{
using T = double;

auto test_point_xs = wrapper<T>({0.0, 0.0});
auto test_point_ys = wrapper<T>({0.0, 0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T>({0.0, 1.0, 0.0});
auto poly_point_ys = wrapper<T>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

TEST_F(PointInPolygonErrorTest, MismatchPolyPointType)
{
using T = double;

auto test_point_xs = wrapper<T>({0.0, 0.0});
auto test_point_ys = wrapper<T>({0.0, 0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<T>({0.0, 1.0, 0.0});
auto poly_point_ys = wrapper<float>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}

TEST_F(PointInPolygonErrorTest, MismatchPointTypes)
{
auto test_point_xs = wrapper<float>({0.0, 0.0});
auto test_point_ys = wrapper<float>({0.0, 0.0});
auto poly_offsets = wrapper<cudf::size_type>({0});
auto poly_ring_offsets = wrapper<cudf::size_type>({0});
auto poly_point_xs = wrapper<double>({0.0, 1.0, 0.0, -1.0});
auto poly_point_ys = wrapper<double>({1.0, 0.0, -1.0, 0.0});

EXPECT_THROW(
cuspatial::point_in_polygon(
test_point_xs, test_point_ys, poly_offsets, poly_ring_offsets, poly_point_xs, poly_point_ys),
cuspatial::logic_error);
}