Fix case where parent bounding box contains hole (#806)

This issue came up when comparing old/new output for countries: When the parent bounding box is fully contained by the outer polygon but fully contains a hole, we were still including the parent in the output. The fix is to check whether the bounding box contains a vertex of the hole, and fail fast if it does.

src/apps/testapps/testPolygonInternal.c

@@ -854,4 +854,47 @@ SUITE(polygonInternal) {
                  "not inside when within hole");
         free(bboxes);
     }
+
+    TEST(cellBoundaryInsidePolygon_notInsideContains) {
+        LatLng verts[] = {{0.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}};
+        GeoLoop geoloop = {.numVerts = 4, .verts = verts};
+
+        GeoPolygon polygon = {.geoloop = geoloop, .numHoles = 0};
+        BBox *bboxes = calloc(sizeof(BBox), 1);
+        bboxesFromGeoPolygon(&polygon, bboxes);
+
+        CellBoundary boundary = {.numVerts = 4,
+                                 .verts = {{0, 0}, {0, 1}, {1, 1}, {1, 0}}};
+        BBox boundaryBBox = {0, 1, 0, 1};
+
+        t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary,
+                                            &boundaryBBox),
+                 "not inside when it contains outer");
+        free(bboxes);
+    }
+
+    TEST(cellBoundaryInsidePolygon_notInsideContainsHole) {
+        LatLng verts[] = {{0, 0}, {0, 1}, {1, 1}, {1, 0}};
+        GeoLoop geoloop = {.numVerts = 4, .verts = verts};
+
+        LatLng holeVerts[] = {{0.6, 0.6}, {0.6, 0.4}, {0.4, 0.4}, {0.4, 0.6}};
+
+        GeoLoop holeGeoLoop = {.numVerts = 4, .verts = holeVerts};
+
+        GeoPolygon polygon = {
+            .geoloop = geoloop, .numHoles = 1, .holes = &holeGeoLoop};
+
+        BBox *bboxes = calloc(sizeof(BBox), 2);
+        bboxesFromGeoPolygon(&polygon, bboxes);
+
+        CellBoundary boundary = {
+            .numVerts = 4,
+            .verts = {{0.9, 0.9}, {0.9, 0.1}, {0.1, 0.1}, {0.1, 0.9}}};
+        BBox boundaryBBox = {0.9, 0.1, 0.9, 0.1};
+
+        t_assert(!cellBoundaryInsidePolygon(&polygon, bboxes, &boundary,
+                                            &boundaryBBox),
+                 "not inside when it contains hole");
+        free(bboxes);
+    }
 }

src/apps/testapps/testPolygonToCellsExperimental.c

@@ -249,6 +249,46 @@ SUITE(polygonToCells) {
         free(hexagons);
     }
 
+    TEST(polygonToCellsHoleParentIssue) {
+        // This checks a specific issue where the bounding box of the parent
+        // cell fully contains the hole.
+        LatLng outer[] = {{0.7774570821346158, 0.19441847890170674},
+                          {0.7528853613617879, 0.19441847890170674},
+                          {0.7528853613617879, 0.23497118026107888},
+                          {0.7774570821346158, 0.23497118026107888}};
+        LatLng sanMarino[] = {{0.7662242554877188, 0.21790879024779208},
+                              {0.7660964275733029, 0.21688101821117023},
+                              {0.7668029019479251, 0.21636628570817204},
+                              {0.7676380769015895, 0.21713838446266925},
+                              {0.7677659048160054, 0.21823092566783267},
+                              {0.7671241996099247, 0.2184218123281233},
+                              {0.7662242554877188, 0.21790879024779208}};
+        GeoPolygon polygon = {
+            .geoloop = {.numVerts = ARRAY_SIZE(outer), .verts = outer},
+            .numHoles = 1,
+            .holes = (GeoLoop[]){
+                {.numVerts = ARRAY_SIZE(sanMarino), .verts = sanMarino}}};
+
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
+            &polygon, 6, CONTAINMENT_CENTER, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &polygon, 6, CONTAINMENT_CENTER, hexagons));
+
+        // This is the cell inside San Marino (i.e. inside the hole)
+        H3Index holeCell = 0x861ea3cefffffff;
+
+        int found = 0;
+        for (int64_t i = 0; i < numHexagons; i++) {
+            if (hexagons[i] == holeCell) found = 1;
+        }
+
+        t_assert(!found, "Did not include cell in hole");
+        free(hexagons);
+    }
+
     TEST(polygonToCellsEmpty) {
         int64_t numHexagons;
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(

src/h3lib/lib/polygon.c

@@ -119,9 +119,20 @@ bool cellBoundaryInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes,
                                    boundaryBBox)) {
         return false;
     }
-    // Check for line intersections with any hole
+
+    // Convert boundary to geoloop for point-inside check
+    const GeoLoop boundaryLoop = {.numVerts = boundary->numVerts,
+                                  // Without this cast, the compiler complains
+                                  // that using const LatLng[] here discards
+                                  // qualifiers. But this should be safe in
+                                  // context, all downstream usage expects const
+                                  .verts = (LatLng *)boundary->verts};
+
+    // Check for line intersections with, or containment of, any hole
     for (int i = 0; i < geoPolygon->numHoles; i++) {
-        if (cellBoundaryCrossesGeoLoop(&(geoPolygon->holes[i]), &bboxes[i + 1],
+        if (pointInsideGeoLoop(&boundaryLoop, boundaryBBox,
+                               &geoPolygon->holes[i].verts[0]) ||
+            cellBoundaryCrossesGeoLoop(&(geoPolygon->holes[i]), &bboxes[i + 1],
                                        boundary, boundaryBBox)) {
             return false;
         }