From 8b21565a8146fa6d2b4c7051cdad14e9616a7c1b Mon Sep 17 00:00:00 2001 From: Nick Rabinowitz Date: Thu, 21 Apr 2022 15:38:03 -0700 Subject: [PATCH 1/3] Fix polyfill bug when vertex latitude exactly matches cell center --- src/apps/testapps/testPolygon.c | 18 ++++++++-- .../testapps/testPolygonToCellsReported.c | 33 +++++++++++++++++++ src/h3lib/include/polygonAlgos.h | 7 ++++ 3 files changed, 55 insertions(+), 3 deletions(-) diff --git a/src/apps/testapps/testPolygon.c b/src/apps/testapps/testPolygon.c index 9f3015d90..4a638da3e 100644 --- a/src/apps/testapps/testPolygon.c +++ b/src/apps/testapps/testPolygon.c @@ -48,10 +48,22 @@ SUITE(polygon) { BBox bbox; bboxFromGeoLoop(&geoloop, &bbox); + // For exact points on the polygon, we bias west and south, + // so only the southwest corner is considered in the polygon t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[0]), - "contains exact"); - t_assert(pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[4]), - "contains exact 4"); + "!contains exact 0"); + t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[1]), + "!contains exact 1"); + t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[2]), + "!contains exact 2"); + t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[4]), + "!contains exact 4"); + t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[5]), + "!contains exact 5"); + + t_assert(pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[3]), + "contains exact 3 (southwest corner)"); + t_assert(pointInsideGeoLoop(&geoloop, &bbox, &inside), "contains point inside"); t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &somewhere), diff --git a/src/apps/testapps/testPolygonToCellsReported.c b/src/apps/testapps/testPolygonToCellsReported.c index 644122f86..8e4d9a9ed 100644 --- a/src/apps/testapps/testPolygonToCellsReported.c +++ b/src/apps/testapps/testPolygonToCellsReported.c @@ -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, ¢erLatLng); + + // 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); + } } diff --git a/src/h3lib/include/polygonAlgos.h b/src/h3lib/include/polygonAlgos.h index b24c99c3d..d1e555d9c 100644 --- a/src/h3lib/include/polygonAlgos.h +++ b/src/h3lib/include/polygonAlgos.h @@ -92,6 +92,13 @@ 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 + if (lat == a.lat || lat == b.lat) { + lat += DBL_EPSILON; + } + // 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) { From 6d5d92b74f3104edb236b1cd56c0f0549313bb76 Mon Sep 17 00:00:00 2001 From: Nick Rabinowitz Date: Fri, 22 Apr 2022 16:58:26 -0700 Subject: [PATCH 2/3] Add an explanatory note about north pole edge case --- src/h3lib/include/polygonAlgos.h | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/h3lib/include/polygonAlgos.h b/src/h3lib/include/polygonAlgos.h index d1e555d9c..e6ce98403 100644 --- a/src/h3lib/include/polygonAlgos.h +++ b/src/h3lib/include/polygonAlgos.h @@ -94,7 +94,14 @@ bool GENERIC_LOOP_ALGO(pointInside)(const TYPE *loop, const BBox *bbox, // 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 + // 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; } From ea14a6b263f2a9f2841b2e96d071d6e572bddc0c Mon Sep 17 00:00:00 2001 From: Nick Rabinowitz Date: Fri, 22 Apr 2022 18:12:28 -0700 Subject: [PATCH 3/3] Add additional point-in-poly tests --- src/apps/testapps/testPolygon.c | 82 +++++++++++++++++++++++++++------ 1 file changed, 69 insertions(+), 13 deletions(-) diff --git a/src/apps/testapps/testPolygon.c b/src/apps/testapps/testPolygon.c index 4a638da3e..df9bc39d7 100644 --- a/src/apps/testapps/testPolygon.c +++ b/src/apps/testapps/testPolygon.c @@ -48,21 +48,11 @@ SUITE(polygon) { BBox bbox; bboxFromGeoLoop(&geoloop, &bbox); - // For exact points on the polygon, we bias west and south, - // so only the southwest corner is considered in the polygon + // For exact points on the polygon, we bias west and south t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[0]), - "!contains exact 0"); - t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[1]), - "!contains exact 1"); - t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[2]), - "!contains exact 2"); - t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[4]), - "!contains exact 4"); - t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[5]), - "!contains exact 5"); - + "does not contain exact vertex 0"); t_assert(pointInsideGeoLoop(&geoloop, &bbox, &sfVerts[3]), - "contains exact 3 (southwest corner)"); + "contains exact vertex 3"); t_assert(pointInsideGeoLoop(&geoloop, &bbox, &inside), "contains point inside"); @@ -70,6 +60,72 @@ SUITE(polygon) { "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},