uber · nrabinowitz · Apr 23, 2022 · Apr 21, 2022 · Apr 22, 2022 · Apr 23, 2022
diff --git a/src/apps/testapps/testPolygon.c b/src/apps/testapps/testPolygon.c
@@ -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},

diff --git 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, &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);
+    }
 }
diff --git a/src/h3lib/include/polygonAlgos.h b/src/h3lib/include/polygonAlgos.h
@@ -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;
+        }
 if (!pointInsidePolygon(geoPolygon, bboxes, &hexCenter)) { 
 pointInsideLinkedGeoLoop(polygons[i]->first, bboxes[i], 
 if (!pointInsidePolygon(geoPolygon, bboxes, &hexCenter)) { 
 pointInsideLinkedGeoLoop(polygons[i]->first, bboxes[i], 
+
         // 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) {