8000 FIX: Fix polyfill bug when vertex latitude exactly matches cell center by nrabinowitz · Pull Request #603 · uber/h3 · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

FIX: Fix polyfill bug when vertex latitude exactly matches cell center #603

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

Merged
merged 3 commits into from
Apr 23, 2022
Merged

FIX: Fix polyfill bug when vertex latitude exactly matches cell center #603

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8b21565
Fix polyfill bug when vertex latitude exactly matches cell center
Apr 21, 2022
6d5d92b
Add an explanatory note about north pole edge case
Apr 22, 2022
ea14a6b
Add additional point-in-poly tests
Apr 23, 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
74 changes: 71 additions & 3 deletions src/apps/testapps/testPolygon.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,16 +48,84 @@ SUITE(polygon) {
BBox bbox;
bboxFromGeoLoop(&geoloop, &bbox);

// For exact points on the polygon, we bias west and south
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[0]),
"contains exact");
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[4]),
"contains exact 4");
"does not contain exact vertex 0");
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[3]),
"contains exact vertex 3");

t_assert(pointInsideGeoLoop(&geoloop, &bbox, &inside),
"contains point inside");
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &somewhere),
"contains somewhere else");
}

TEST(pointInsideGeoLoopCornerCases) {
LatLng verts[] = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
GeoLoop geoloop = {.numVerts = 4, .verts = verts};

BBox bbox;
bboxFromGeoLoop(&geoloop, &bbox);

LatLng point = {0, 0};

// Test corners. For exact points on the polygon, we bias west and
// north, so only the southeast corner is contained.
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
"does not contain sw corner");
point.lat = 1;
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
"does not contain nw corner ");
point.lng = 1;
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
"does not contain ne corner ");
point.lat = 0;
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &point),
"contains se corner ");
}

TEST(pointInsideGeoLoopEdgeCases) {
LatLng verts[] = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
GeoLoop geoloop = {.numVerts = 4, .verts = verts};

BBox bbox;
bboxFromGeoLoop(&geoloop, &bbox);

LatLng point;

// Test edges. Only points on south and east edges are contained.
point.lat = 0.5;
point.lng = 0;
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
"does not contain point on west edge");
point.lat = 1;
point.lng = 0.5;
t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
"does not contain point on north edge");
point.lat = 0.5;
point.lng = 1;
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &point),
"contains point on east edge");
point.lat = 0;
point.lng = 0.5;
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &point),
"contains point on south edge");
}

TEST(pointInsideGeoLoopExtraEdgeCase) {
// This is a carefully crafted shape + point to hit an otherwise
// missed branch in coverage
LatLng verts[] = {{0, 0}, {1, 0.5}, {0, 1}};
GeoLoop geoloop = {.numVerts = 4, .verts = verts};

BBox bbox;
bboxFromGeoLoop(&geoloop, &bbox);

LatLng point = {0.5, 0.5};
t_assert(pointInsideGeoLoop(&geoloop, &bbox, &point),
"contains inside point matching longitude of a vertex");
}

TEST(pointInsideGeoLoopTransmeridian) {
LatLng verts[] = {{0.01, -M_PI + 0.01},
{0.01, M_PI - 0.01},
Expand Down
33 changes: 33 additions & 0 deletions src/apps/testapps/testPolygonToCellsReported.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,4 +170,37 @@ SUITE(polygonToCells_reported) {
t_assert(actualNumIndexes == 4353, "got expected polygonToCells size");
free(hexagons);
}

// https://github.com/uber/h3/issues/595
TEST(h3_136) {
H3Index center = 0x85283473fffffff;
LatLng centerLatLng;
H3_EXPORT(cellToLatLng)(center, &centerLatLng);

// This polygon should include the center cell. The issue here arises
// when one of the polygon vertexes is to the east of the index center,
// with exactly the same latitude
LatLng testVerts[] = {{.lat = centerLatLng.lat, -2.121207808248113},
{0.6565301558937859, -2.1281107217935986},
{0.6515463604919347, -2.1345342663428695},
{0.6466583305904194, -2.1276313527973842}};

GeoLoop testGeoLoop = {.numVerts = 4, .verts = testVerts};
GeoPolygon testPolygon;
testPolygon.geoloop = testGeoLoop;
testPolygon.numHoles = 0;

int res = 5;
int64_t numHexagons;
t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0,
&numHexagons));
H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));

t_assertSuccess(
H3_EXPORT(polygonToCells)(&testPolygon, res, 0, hexagons));
int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);

t_assert(actualNumIndexes == 8, "got expected polygonToCells size");
free(hexagons);
}
}
14 changes: 14 additions & 0 deletions src/h3lib/include/polygonAlgos.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,20 @@ bool GENERIC_LOOP_ALGO(pointInside)(const TYPE *loop, const BBox *bbox,
b = tmp;
}

// If the latitude matches exactly, we'll hit an edge case where
// the ray passes through the vertex twice on successive segment
// checks. To avoid this, adjust the latiude northward if needed.
//
// NOTE: This currently means that a point at the north pole cannot
// be contained in any polygon. This is acceptable in current usage,
// because the point we test in this function at present is always
// a cell center or vertex, and no cell has a center or vertex on the
// north pole. If we need to expand this algo to more generic uses we
// might need to handle this edge case.
if (lat == a.lat || lat == b.lat) {
lat += DBL_EPSILON;
}
Comment on lines +105 to +107
Copy link
Collaborator

Choose a reason for hiding this comment

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

What happens if the latitude is maximally north already? I assume there won't be any issue here?

Copy link
Collaborator Author
@nrabinowitz nrabinowitz Apr 22, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

You may be correct - if we're checking the latlng point [90, 0] (i.e. the north pole), this algorithm will I think fail to contain it in any polygon. But because we're using Cartesian polygons to begin with, there's no pole-enclosing polygon anyway. This is all moot for the polyfill usage, however, because the north pole cell at res 15 has a center latitude of 89.99999581738042, so it's not possible for this issue to arise. If we offered point-in-poly as an external function, this might be an issue, but probably one we could address in documentation.

I suppose we could address this issue with

if (lat > M_PI_2) {
  lat = M_PI_2 - DBL_EPSILON
}

but I don't think this is necessary for our current usage, and not worth the extra branch.

Copy link
Collaborator

Choose a reason for hiding this comment

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

So the lat and lng here is provided by the caller, so this is theoretically a problem, but we just use the centerpoint of the polygon in reality:

h3/src/h3lib/lib/algos.c

Line 1004 in ba89a19

if (!pointInsidePolygon(geoPolygon, bboxes, &hexCenter)) {

There might be a problem with pointInsideLinkedGeoLoop, though, since it just uses the first point which could be the northernmost one:

h3/src/h3lib/lib/linkedGeo.c

Line 191 in adde6da

pointInsideLinkedGeoLoop(polygons[i]->first, bboxes[i],

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The pointInsideLinkedGeoLoop function is looking at whether one loop is inside another. So for this issue to occur, you'd need a fairly pathological shape, where not only was one vertex of the outer loop on the north pole, but one vertex of the inner loop was also on the north pole, eg.

image

But again, the loops involved here are derived from H3 cell boundaries, and there is no cell with a center or a vertex on the north pole.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Perhaps add a comment explaining the assumption that only cell centers or vertexes are expected to be compared here against the user input?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Added a long comment noting the issue.


// If we're totally above or below the latitude ranges, the test
// ray cannot intersect the line segment, so let's move on
if (lat < a.lat || lat > b.lat) {
Expand Down
0