diff --git a/CHANGELOG.md b/CHANGELOG.md index 9312eb47f..f37236ec5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ The public API of this library consists of the functions declared in file [h3api.h](./src/h3lib/include/h3api.h). ## [Unreleased] +### Fixed +- Fixed bounding box bug for polygons crossing the antimeridian (#130) ### Changed - Longitude outputs are now guaranteed to be in the range [-Pi, Pi]. (#93) diff --git a/src/apps/testapps/testBBox.c b/src/apps/testapps/testBBox.c index 916a91bdf..7057467cd 100644 --- a/src/apps/testapps/testBBox.c +++ b/src/apps/testapps/testBBox.c @@ -88,9 +88,16 @@ TEST(transmeridian) { {-0.4, M_PI - 0.1}}; const Geofence geofence = {.numVerts = 4, .verts = verts}; const BBox expected = {0.4, -0.4, -M_PI + 0.1, M_PI - 0.1}; - const GeoCoord inside = {-0.1, M_PI}; + const GeoCoord insideOnMeridian = {-0.1, M_PI}; const GeoCoord outside = {1.0, M_PI - 0.5}; - assertBBox(&geofence, &expected, &inside, &outside); + assertBBox(&geofence, &expected, &insideOnMeridian, &outside); + + const GeoCoord westInside = {0.1, M_PI - 0.05}; + t_assert(bboxContains(&expected, &westInside), + "Contains expected west inside point"); + const GeoCoord eastInside = {0.1, -M_PI + 0.05}; + t_assert(bboxContains(&expected, &eastInside), + "Contains expected east outside point"); const GeoCoord westOutside = {0.1, M_PI - 0.5}; t_assert(!bboxContains(&expected, &westOutside), diff --git a/src/apps/testapps/testPolyfill.c b/src/apps/testapps/testPolyfill.c index 47408c898..74ba9b02c 100644 --- a/src/apps/testapps/testPolyfill.c +++ b/src/apps/testapps/testPolyfill.c @@ -26,76 +26,43 @@ GeoCoord sfVerts[] = { {0.659966917655, -2.1364398519396}, {0.6595011102219, -2.1359434279405}, {0.6583348114025, -2.1354884206045}, {0.6581220034068, -2.1382437718946}, {0.6594479998527, -2.1384597563896}, {0.6599990002976, -2.1376771158464}}; -Geofence sfGeofence; +Geofence sfGeofence = {.numVerts = 6, .verts = sfVerts}; GeoPolygon sfGeoPolygon; + GeoCoord holeVerts[] = {{0.6595072188743, -2.1371053983433}, {0.6591482046471, -2.1373141048153}, {0.6592295020837, -2.1365222838402}}; -Geofence holeGeofence; +Geofence holeGeofence = {.numVerts = 3, .verts = holeVerts}; GeoPolygon holeGeoPolygon; + GeoCoord emptyVerts[] = {{0.659966917655, -2.1364398519394}, {0.659966917655, -2.1364398519395}, {0.659966917655, -2.1364398519396}}; -Geofence emptyGeofence; +Geofence emptyGeofence = {.numVerts = 3, .verts = emptyVerts}; GeoPolygon emptyGeoPolygon; -GeoCoord primeMeridianVerts[] = { - {0.01, 0.01}, {0.01, -0.01}, {-0.01, -0.01}, {-0.01, 0.01}}; -Geofence primeMeridianGeofence; -GeoPolygon primeMeridianGeoPolygon; - -GeoCoord transMeridianVerts[] = {{0.01, -M_PI + 0.01}, - {0.01, M_PI - 0.01}, - {-0.01, M_PI - 0.01}, - {-0.01, -M_PI + 0.01}}; -Geofence transMeridianGeofence; -GeoPolygon transMeridianGeoPolygon; - -GeoCoord transMeridianHoleVerts[] = {{0.005, -M_PI + 0.005}, - {0.005, M_PI - 0.005}, - {-0.005, M_PI - 0.005}, - {-0.005, -M_PI + 0.005}}; -Geofence transMeridianHoleGeofence; -GeoPolygon transMeridianHoleGeoPolygon; -GeoPolygon transMeridianFilledHoleGeoPolygon; +static int countActualHexagons(H3Index* hexagons, int numHexagons) { + int actualNumHexagons = 0; + for (int i = 0; i < numHexagons; i++) { + if (hexagons[i] != 0) { + actualNumHexagons++; + } + } + return actualNumHexagons; +} BEGIN_TESTS(polyfill); -sfGeofence.numVerts = 6; -sfGeofence.verts = sfVerts; sfGeoPolygon.geofence = sfGeofence; sfGeoPolygon.numHoles = 0; -holeGeofence.numVerts = 3; -holeGeofence.verts = holeVerts; holeGeoPolygon.geofence = sfGeofence; holeGeoPolygon.numHoles = 1; holeGeoPolygon.holes = &holeGeofence; -emptyGeofence.numVerts = 3; -emptyGeofence.verts = emptyVerts; emptyGeoPolygon.geofence = emptyGeofence; emptyGeoPolygon.numHoles = 0; -primeMeridianGeofence.numVerts = 4; -primeMeridianGeofence.verts = primeMeridianVerts; -primeMeridianGeoPolygon.geofence = primeMeridianGeofence; -primeMeridianGeoPolygon.numHoles = 0; - -transMeridianGeofence.numVerts = 4; -transMeridianGeofence.verts = transMeridianVerts; -transMeridianGeoPolygon.geofence = transMeridianGeofence; -transMeridianGeoPolygon.numHoles = 0; - -transMeridianHoleGeofence.numVerts = 4; -transMeridianHoleGeofence.verts = transMeridianHoleVerts; -transMeridianHoleGeoPolygon.geofence = transMeridianGeofence; -transMeridianHoleGeoPolygon.numHoles = 1; -transMeridianHoleGeoPolygon.holes = &transMeridianHoleGeofence; - -transMeridianFilledHoleGeoPolygon.geofence = transMeridianHoleGeofence; -transMeridianFilledHoleGeoPolygon.numHoles = 0; - TEST(maxPolyfillSize) { int numHexagons = H3_EXPORT(maxPolyfillSize)(&sfGeoPolygon, 9); t_assert(numHexagons == 3169, "got expected max polyfill size"); @@ -112,12 +79,7 @@ TEST(polyfill) { H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&sfGeoPolygon, 9, hexagons); - int actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagons[i] != 0) { - actualNumHexagons++; - } - } + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); t_assert(actualNumHexagons == 1253, "got expected polyfill size"); free(hexagons); @@ -128,12 +90,7 @@ TEST(polyfillHole) { H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&holeGeoPolygon, 9, hexagons); - int actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagons[i] != 0) { - actualNumHexagons++; - } - } + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); t_assert(actualNumHexagons == 1214, "got expected polyfill size (hole)"); free(hexagons); @@ -144,12 +101,7 @@ TEST(polyfillEmpty) { H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&emptyGeoPolygon, 9, hexagons); - int actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagons[i] != 0) { - actualNumHexagons++; - } - } + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); t_assert(actualNumHexagons == 0, "got expected polyfill size (empty)"); free(hexagons); @@ -178,13 +130,7 @@ TEST(polyfillExact) { H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&someHexagon, 9, hexagons); - int actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagons[i] != 0) { - t_assert(hexagons[i] == origin, "Got origin back"); - actualNumHexagons++; - } - } + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); t_assert(actualNumHexagons == 1, "got expected polyfill size (1)"); free(hexagons); @@ -192,6 +138,35 @@ TEST(polyfillExact) { } TEST(polyfillTransmeridian) { + GeoCoord primeMeridianVerts[] = { + {0.01, 0.01}, {0.01, -0.01}, {-0.01, -0.01}, {-0.01, 0.01}}; + Geofence primeMeridianGeofence = {.numVerts = 4, + .verts = primeMeridianVerts}; + GeoPolygon primeMeridianGeoPolygon = {.geofence = primeMeridianGeofence, + .numHoles = 0}; + + GeoCoord transMeridianVerts[] = {{0.01, -M_PI + 0.01}, + {0.01, M_PI - 0.01}, + {-0.01, M_PI - 0.01}, + {-0.01, -M_PI + 0.01}}; + Geofence transMeridianGeofence = {.numVerts = 4, + .verts = transMeridianVerts}; + GeoPolygon transMeridianGeoPolygon = {.geofence = transMeridianGeofence, + .numHoles = 0}; + + GeoCoord transMeridianHoleVerts[] = {{0.005, -M_PI + 0.005}, + {0.005, M_PI - 0.005}, + {-0.005, M_PI - 0.005}, + {-0.005, -M_PI + 0.005}}; + Geofence transMeridianHoleGeofence = {.numVerts = 4, + .verts = transMeridianHoleVerts}; + GeoPolygon transMeridianHoleGeoPolygon = { + .geofence = transMeridianGeofence, + .numHoles = 1, + .holes = &transMeridianHoleGeofence}; + GeoPolygon transMeridianFilledHoleGeoPolygon = { + .geofence = transMeridianHoleGeofence, .numHoles = 0}; + int expectedSize; // Prime meridian case @@ -200,12 +175,7 @@ TEST(polyfillTransmeridian) { H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&primeMeridianGeoPolygon, 7, hexagons); - int actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagons[i] != 0) { - actualNumHexagons++; - } - } + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); t_assert(actualNumHexagons == expectedSize, "got expected polyfill size (prime meridian)"); @@ -218,12 +188,7 @@ TEST(polyfillTransmeridian) { H3Index* hexagonsTM = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&transMeridianGeoPolygon, 7, hexagonsTM); - actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagonsTM[i] != 0) { - actualNumHexagons++; - } - } + actualNumHexagons = countActualHexagons(hexagonsTM, numHexagons); t_assert(actualNumHexagons == expectedSize, "got expected polyfill size (transmeridian)"); @@ -234,24 +199,14 @@ TEST(polyfillTransmeridian) { H3Index* hexagonsTMFH = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&transMeridianFilledHoleGeoPolygon, 7, hexagonsTMFH); - int actualNumHoleHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagonsTMFH[i] != 0) { - actualNumHoleHexagons++; - } - } + int actualNumHoleHexagons = countActualHexagons(hexagonsTMFH, numHexagons); // Transmeridian hole case numHexagons = H3_EXPORT(maxPolyfillSize)(&transMeridianHoleGeoPolygon, 7); H3Index* hexagonsTMH = calloc(numHexagons, sizeof(H3Index)); H3_EXPORT(polyfill)(&transMeridianHoleGeoPolygon, 7, hexagonsTMH); - actualNumHexagons = 0; - for (int i = 0; i < numHexagons; i++) { - if (hexagonsTMH[i] != 0) { - actualNumHexagons++; - } - } + actualNumHexagons = countActualHexagons(hexagonsTMH, numHexagons); t_assert(actualNumHexagons == expectedSize - actualNumHoleHexagons, "got expected polyfill size (transmeridian hole)"); @@ -262,6 +217,29 @@ TEST(polyfillTransmeridian) { free(hexagonsTMH); } +TEST(polyfillTransmeridianComplex) { + // This polygon is "complex" in that it has > 4 vertices - this + // tests for a bug that was taking the max and min longitude as + // the bounds for transmeridian polygons + GeoCoord verts[] = {{0.1, -M_PI + 0.00001}, {0.1, M_PI - 0.00001}, + {0.05, M_PI - 0.2}, {-0.1, M_PI - 0.00001}, + {-0.1, -M_PI + 0.00001}, {-0.05, -M_PI + 0.2}}; + Geofence geofence = {.numVerts = 6, .verts = verts}; + GeoPolygon polygon = {.geofence = geofence, .numHoles = 0}; + + int numHexagons = H3_EXPORT(maxPolyfillSize)(&polygon, 4); + + H3Index* hexagons = calloc(numHexagons, sizeof(H3Index)); + H3_EXPORT(polyfill)(&polygon, 4, hexagons); + + int actualNumHexagons = countActualHexagons(hexagons, numHexagons); + + t_assert(actualNumHexagons == 1204, + "got expected polyfill size (complex transmeridian)"); + + free(hexagons); +} + TEST(polyfillPentagon) { H3Index pentagon; setH3Index(&pentagon, 9, 24, 0); diff --git a/src/apps/testapps/testPolygon.c b/src/apps/testapps/testPolygon.c index f7e8a0c6d..06361607c 100644 --- a/src/apps/testapps/testPolygon.c +++ b/src/apps/testapps/testPolygon.c @@ -127,10 +127,7 @@ TEST(pointInsideLinkedGeoLoop) { TEST(bboxFromGeofence) { GeoCoord verts[] = {{0.8, 0.3}, {0.7, 0.6}, {1.1, 0.7}, {1.0, 0.2}}; - - Geofence geofence; - geofence.verts = verts; - geofence.numVerts = 4; + Geofence geofence = {.numVerts = 4, .verts = verts}; const BBox expected = {1.1, 0.7, 0.7, 0.2}; @@ -139,6 +136,19 @@ TEST(bboxFromGeofence) { t_assert(bboxEquals(&result, &expected), "Got expected bbox"); } +TEST(bboxFromGeofenceTransmeridian) { + GeoCoord verts[] = {{0.1, -M_PI + 0.1}, {0.1, M_PI - 0.1}, + {0.05, M_PI - 0.2}, {-0.1, M_PI - 0.1}, + {-0.1, -M_PI + 0.1}, {-0.05, -M_PI + 0.2}}; + Geofence geofence = {.numVerts = 6, .verts = verts}; + + const BBox expected = {0.1, -0.1, -M_PI + 0.2, M_PI - 0.2}; + + BBox result; + bboxFromGeofence(&geofence, &result); + t_assert(bboxEquals(&result, &expected), "Got expected transmeridian bbox"); +} + TEST(bboxFromGeofenceNoVertices) { Geofence geofence; geofence.verts = NULL; @@ -154,14 +164,8 @@ TEST(bboxFromGeofenceNoVertices) { TEST(bboxesFromGeoPolygon) { GeoCoord verts[] = {{0.8, 0.3}, {0.7, 0.6}, {1.1, 0.7}, {1.0, 0.2}}; - - Geofence geofence; - geofence.verts = verts; - geofence.numVerts = 4; - - GeoPolygon polygon; - polygon.geofence = geofence; - polygon.numHoles = 0; + Geofence geofence = {.numVerts = 4, .verts = verts}; + GeoPolygon polygon = {.geofence = geofence, .numHoles = 0}; const BBox expected = {1.1, 0.7, 0.7, 0.2}; @@ -174,22 +178,14 @@ TEST(bboxesFromGeoPolygon) { TEST(bboxesFromGeoPolygonHole) { GeoCoord verts[] = {{0.8, 0.3}, {0.7, 0.6}, {1.1, 0.7}, {1.0, 0.2}}; - - Geofence geofence; - geofence.verts = verts; - geofence.numVerts = 4; + Geofence geofence = {.numVerts = 4, .verts = verts}; // not a real hole, but doesn't matter for the test GeoCoord holeVerts[] = {{0.9, 0.3}, {0.9, 0.5}, {1.0, 0.7}, {0.9, 0.3}}; + Geofence holeGeofence = {.numVerts = 4, .verts = holeVerts}; - Geofence holeGeofence; - holeGeofence.verts = holeVerts; - holeGeofence.numVerts = 4; - - GeoPolygon polygon; - polygon.geofence = geofence; - polygon.numHoles = 1; - polygon.holes = &holeGeofence; + GeoPolygon polygon = { + .geofence = geofence, .numHoles = 1, .holes = &holeGeofence}; const BBox expected = {1.1, 0.7, 0.7, 0.2}; const BBox expectedHole = {1.0, 0.9, 0.7, 0.3}; @@ -237,10 +233,7 @@ TEST(bboxFromLinkedGeoLoopNoVertices) { TEST(isClockwiseGeofence) { GeoCoord verts[] = {{0, 0}, {0.1, 0.1}, {0, 0.1}}; - - Geofence geofence; - geofence.verts = verts; - geofence.numVerts = 3; + Geofence geofence = {.numVerts = 3, .verts = verts}; t_assert(isClockwiseGeofence(&geofence), "Got true for clockwise geofence"); } diff --git a/src/h3lib/include/polygonAlgos.h b/src/h3lib/include/polygonAlgos.h index 6b2b2e55e..0681cc850 100644 --- a/src/h3lib/include/polygonAlgos.h +++ b/src/h3lib/include/polygonAlgos.h @@ -142,15 +142,17 @@ void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE* loop, BBox* bbox) { bbox->west = 0; return; } - double lat; - double lon; bbox->south = DBL_MAX; bbox->west = DBL_MAX; - bbox->north = -1.0 * DBL_MAX; - bbox->east = -1.0 * DBL_MAX; + bbox->north = -DBL_MAX; + bbox->east = -DBL_MAX; + double minPosLon = DBL_MAX; + double maxNegLon = -DBL_MAX; bool isTransmeridian = false; + double lat; + double lon; GeoCoord coord; GeoCoord next; @@ -165,6 +167,10 @@ void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE* loop, BBox* bbox) { if (lon < bbox->west) bbox->west = lon; if (lat > bbox->north) bbox->north = lat; if (lon > bbox->east) bbox->east = lon; + // Save the min positive and max negative longitude for + // use in the transmeridian case + if (lon > 0 && lon < minPosLon) minPosLon = lon; + if (lon < 0 && lon > maxNegLon) maxNegLon = lon; // check for arcs > 180 degrees longitude, flagging as transmeridian if (fabs(lon - next.lon) > M_PI) { isTransmeridian = true; @@ -172,9 +178,8 @@ void GENERIC_LOOP_ALGO(bboxFrom)(const TYPE* loop, BBox* bbox) { } // Swap east and west if transmeridian if (isTransmeridian) { - double tmp = bbox->east; - bbox->east = bbox->west; - bbox->west = tmp; + bbox->east = maxNegLon; + bbox->west = minPosLon; } }